Vous ne trouvez pas de réponse à votre problème ? Alors posez la question dans le forum. Souvenez-vous qu'il n'y a jamais de question bête, mais rester dans l'ignorance parce que l'on n'ose pas poser une question, ça c'est une erreur !

MODULE DE FONCTIONS STATISTIQUES EN VB6


Information sur la source

Catégorie :Maths Niveau : Débutant Date de création : 18/05/2005 Date de mise à jour : 18/05/2005 10:28:45 Vu : 5 807

Note :
7 / 10 - par 3 personnes
7,00 / 10

  • 1

  • 2

  • 3

  • 4

  • 5

  • 6

  • 7

  • 8

  • 9

  • 10

Commentaire sur cette source (3)
Ajouter un commentaire et/ou une note

Description

Ce code contient les fonctions de Statistiques suivantes :
- Fonction de la table inverse de la loi normale centrée réduite
- Fonction renvoyant la probabilité d'une variable aléatoire z suivant la loi normale centrée réduite
- Fonction de la loi du Khi-deux
- Fonction de la loi inverse du Khi-deux
 

Source

  • ' L'algorithme du codage de la loi inverse du Chi-deux provient de la traduction d'un code source
  • ' en JavaScript sur le site internet suivant :
  • ' http://www.fourmilab.ch/rpkp/experiments/analysis/chiCalc.html
  • 'The following JavaScript functions for calculating normal and
  • 'chi-square probabilities and critical values were adapted by
  • 'John Walker from C implementations
  • 'written by Gary Perlman of Wang Institute, Tyngsboro, MA 01879.
  • 'Both the original C code and this JavaScript edition
  • 'are in the public domain.
  • 'Public Const Pi = 3.14159265359
  • 'densité de probabilité de la loi normale centrée réduite
  • Function PHI(ByVal u As Double) As Double
  • PHI = (1 / ((2 * Pi) ^ 0.5)) * Exp(-0.5 * (u ^ 2))
  • End Function
  • 'Fonction de la table inverse de la loi normale centrée réduite
  • Function NORMAL(ByVal p As Double) As Double
  • Dim xo, xx, y, h, f, f1, f2, fd, r, s As Double
  • Dim N, i, k As Integer
  • xo = 0.5
  • y = 0
  • xx = xo
  • N = 200
  • 'Suite convergente de Newton
  • Do While Abs(xx - y) > 0.00001
  • y = xx
  • h = 2 * xx / N
  • f1 = 0
  • f2 = 0
  • 'Résolution d'une intégrale numérique
  • For i = 0 To ((N - 2) / 2)
  • r = -xx + 2 * i * h
  • s = -xx + (2 * i + 1) * h
  • If i = 0 Then
  • f1 = f1 + 0
  • f2 = f2 + PHI(s)
  • Else
  • f1 = f1 + PHI(r)
  • f2 = f2 + PHI(s)
  • End If
  • Next
  • f = (h / 3) * (PHI(-xx) + PHI(xx) + 2 * f1 + 4 * f2) - p
  • fd = 2 * PHI(xx)
  • xx = xx - f / fd
  • Loop
  • NORMAL = xx
  • End Function
  • ' Renvoie la probabilité d'une variable aléatoire z suivant la loi normale centrée réduite
  • Function poz(ByVal z As Double) As Double
  • 'POZ -- probability of normal z value
  • 'Adapted from a polynomial approximation in:
  • 'Ibbetson D, Algorithm 209
  • 'Collected Algorithms of the CACM 1963 p. 616
  • 'Note:
  • 'This routine has six digit accuracy, so it is only useful for absolute
  • 'z values < 6. For z values >= to 6.0, poz() returns 0.0.
  • Dim y, x, w As Double
  • Dim z_max As Double
  • z_max = 6
  • If z = 0 Then
  • x = 0
  • Else
  • y = 0.5 * Abs(z)
  • If y >= z_max * 0.5 Then
  • x = 1
  • ElseIf y < 1 Then
  • w = y * y
  • x = ((((((((0.000124818987 * w - 0.001075204047) * w _
  • + 0.005198775019) * w - 0.019198292004) * w _
  • + 0.059054035642) * w - 0.151968751364) * w _
  • + 0.319152932694) * w - 0.5319230073) * w _
  • + 0.797884560593) * y * 2
  • Else
  • y = y - 2
  • x = (((((((((((((-0.000045255659 * y + 0.00015252929) * y _
  • - 0.000019538132) * y - 0.000676904986) * y _
  • + 0.001390604284) * y - 0.00079462082) * y _
  • - 0.002034254874) * y + 0.006549791214) * y _
  • - 0.010557625006) * y + 0.011630447319) * y _
  • - 0.009279453341) * y + 0.005353579108) * y _
  • - 0.002141268741) * y + 0.000535310849) * y _
  • + 0.999936657524
  • End If
  • End If
  • If z > 0 Then
  • poz = ((x + 1) * 0.5)
  • Else
  • poz = ((1 - x) * 0.5)
  • End If
  • Dim bigx As Double
  • bigx = 20
  • End Function
  • Function pochisq(ByVal x As Double, ByVal df As Integer) As Double
  • ' Adapted From:
  • 'Hill, I. D. and Pike, M. C. Algorithm 299
  • 'Collected Algorithms for the CACM 1967 p. 243
  • 'Updated for rounding errors based on remark in
  • 'ACM TOMS June 1985, page 185
  • Dim a, y, s, e, c, v As Double
  • Dim even As Boolean
  • 'even correspond à la parité de df, le degré de liberté
  • Dim lnpi, ipi As Double
  • lnpi = Log(Pi ^ 0.5)
  • ipi = 1 / Log(Pi)
  • If ((x <= 0) Or (df < 1)) Then
  • pochisq = 1
  • End If
  • a = 0.5 * x
  • If Fix(df / 2) = df / 2 Then
  • even = True
  • Else
  • even = False
  • End If
  • If df > 1 Then
  • y = Exp(-a)
  • End If
  • If even = True Then
  • s = y
  • Else
  • s = 2 * poz(-(x ^ 0.5))
  • End If
  • If df > 2 Then
  • x = 0.5 * (df - 1)
  • If even = True Then
  • z = 1
  • Else
  • z = 0.5
  • End If
  • If (a > bigx) Then
  • If even = True Then
  • e = 0
  • Else
  • e = lnpi
  • End If
  • c = Log(a)
  • Do While (z <= x)
  • e = Log(z) + e
  • s = s + Exp(c * z - a - e)
  • z = z + 1
  • Loop
  • pochisq = s
  • Else
  • If even = True Then
  • e = 1
  • Else
  • e = ipi / (a ^ 0.5)
  • End If
  • c = 0
  • Do While (z <= x)
  • e = e * (a / z)
  • c = c + e
  • z = z + 1
  • Loop
  • pochisq = c * y + s
  • End If
  • Else
  • pochisq = s
  • End If
  • End Function
  • Function critchi(ByVal p As Double, ByVal df As Integer) As Double
  • Dim epsilon, chimax, minchisq, maxchisq, chisqval As Double
  • epsilon = 0.000001
  • chimax = 99999
  • minchisq = 0
  • maxchisq = chimax
  • If p <= 0 Then
  • critchi = maxchisq
  • Else
  • If p >= 1 Then
  • critchi = 0
  • End If
  • End If
  • chisqval = df / (p ^ 0.5)
  • Do While ((maxchisq - minchisq) > epsilon)
  • If (pochisq(chisqval, df) < p) Then
  • maxchisq = chisqval
  • Else
  • minchisq = chisqval
  • End If
  • chisqval = (maxchisq + minchisq) * 0.5
  • Loop
  • critchi = chisqval
  • End Function
  • 'renvoie, pour la probabilité p et le degré de liberté df, la variable
  • Function invchi2(ByVal p As Double, ByVal df As Double) As Double
  • invchi2 = critchi(p, df)
  • End Function
  • 'renvoie, pour la variable x et le degré de liberté df, la probabilité
  • Function chi2(ByVal x As Double, ByVal df As Integer) As Double
  • chi2 = pochisq(x, df)
  • End Function
' L'algorithme du codage de la loi inverse du Chi-deux provient de la traduction d'un code source
' en JavaScript sur le site internet suivant :
'   http://www.fourmilab.ch/rpkp/experiments/analysis/chiCalc.html

	'The following JavaScript functions for calculating normal and
	'chi-square probabilities and critical values were adapted by
	'John Walker from C implementations
	'written by Gary Perlman of Wang Institute, Tyngsboro, MA 01879.
	'Both the original C code and this JavaScript edition
	'are in the public domain.



'Public Const Pi = 3.14159265359




'densité de probabilité de la loi normale centrée réduite
Function PHI(ByVal u As Double) As Double
    PHI = (1 / ((2 * Pi) ^ 0.5)) * Exp(-0.5 * (u ^ 2))
End Function



'Fonction de la table inverse de la loi normale centrée réduite
Function NORMAL(ByVal p As Double) As Double

    Dim xo, xx, y, h, f, f1, f2, fd, r, s As Double
    Dim N, i, k As Integer
    xo = 0.5
    y = 0
    xx = xo
    N = 200
    'Suite convergente de Newton
    Do While Abs(xx - y) > 0.00001
        y = xx
        h = 2 * xx / N
        f1 = 0
        f2 = 0
        'Résolution d'une intégrale numérique
        For i = 0 To ((N - 2) / 2)
            r = -xx + 2 * i * h
            s = -xx + (2 * i + 1) * h
            If i = 0 Then
                f1 = f1 + 0
                f2 = f2 + PHI(s)
            Else
                f1 = f1 + PHI(r)
                f2 = f2 + PHI(s)
            End If
        Next
        f = (h / 3) * (PHI(-xx) + PHI(xx) + 2 * f1 + 4 * f2) - p
        fd = 2 * PHI(xx)
        xx = xx - f / fd
    Loop
    NORMAL = xx
    
End Function



' Renvoie la probabilité d'une variable aléatoire z suivant la loi normale centrée réduite
Function poz(ByVal z As Double) As Double
'POZ  --  probability of normal z value

'Adapted from a polynomial approximation in:
'Ibbetson D, Algorithm 209
'Collected Algorithms of the CACM 1963 p. 616

'Note:
    'This routine has six digit accuracy, so it is only useful for absolute
    'z values < 6.  For z values >= to 6.0, poz() returns 0.0.
    
    
    Dim y, x, w As Double
    Dim z_max As Double
    z_max = 6
    If z = 0 Then
        x = 0
    Else
        y = 0.5 * Abs(z)
        If y >= z_max * 0.5 Then
            x = 1
        ElseIf y < 1 Then
            w = y * y
            x = ((((((((0.000124818987 * w - 0.001075204047) * w _
                    + 0.005198775019) * w - 0.019198292004) * w _
                    + 0.059054035642) * w - 0.151968751364) * w _
                    + 0.319152932694) * w - 0.5319230073) * w _
                    + 0.797884560593) * y * 2
        Else
            y = y - 2
            x = (((((((((((((-0.000045255659 * y + 0.00015252929) * y _
                    - 0.000019538132) * y - 0.000676904986) * y _
                    + 0.001390604284) * y - 0.00079462082) * y _
                    - 0.002034254874) * y + 0.006549791214) * y _
                    - 0.010557625006) * y + 0.011630447319) * y _
                    - 0.009279453341) * y + 0.005353579108) * y _
                    - 0.002141268741) * y + 0.000535310849) * y _
                    + 0.999936657524
        End If
    End If
    If z > 0 Then
        poz = ((x + 1) * 0.5)
    Else
        poz = ((1 - x) * 0.5)
    End If
    Dim bigx As Double
    bigx = 20
    
End Function




Function pochisq(ByVal x As Double, ByVal df As Integer) As Double
' Adapted From:
'Hill, I. D. and Pike, M. C.  Algorithm 299
'Collected Algorithms for the CACM 1967 p. 243
'Updated for rounding errors based on remark in
'ACM TOMS June 1985, page 185
    
    Dim a, y, s, e, c, v As Double
    Dim even As Boolean
    'even correspond à la parité de df, le degré de liberté
    Dim lnpi, ipi As Double
    lnpi = Log(Pi ^ 0.5)
    ipi = 1 / Log(Pi)
    If ((x <= 0) Or (df < 1)) Then
        pochisq = 1
    End If
    a = 0.5 * x
    If Fix(df / 2) = df / 2 Then
        even = True
    Else
        even = False
    End If
    If df > 1 Then
        y = Exp(-a)
    End If
    If even = True Then
        s = y
    Else
        s = 2 * poz(-(x ^ 0.5))
    End If
    If df > 2 Then
        x = 0.5 * (df - 1)
        If even = True Then
            z = 1
        Else
            z = 0.5
        End If
        If (a > bigx) Then
            If even = True Then
                e = 0
            Else
                e = lnpi
            End If
            
            c = Log(a)
            
            Do While (z <= x)
                e = Log(z) + e
                s = s + Exp(c * z - a - e)
                z = z + 1
            Loop
            pochisq = s
        Else
            If even = True Then
                e = 1
            Else
                e = ipi / (a ^ 0.5)
            End If
            c = 0
            Do While (z <= x)
                e = e * (a / z)
                c = c + e
                z = z + 1
            Loop
            pochisq = c * y + s
        End If
    Else
        pochisq = s
    End If
    
End Function


Function critchi(ByVal p As Double, ByVal df As Integer) As Double

    Dim epsilon, chimax, minchisq, maxchisq, chisqval As Double
    epsilon = 0.000001
    chimax = 99999
    minchisq = 0
    maxchisq = chimax
    If p <= 0 Then
        critchi = maxchisq
    Else
        If p >= 1 Then
            critchi = 0
        End If
    End If
    chisqval = df / (p ^ 0.5)
    Do While ((maxchisq - minchisq) > epsilon)
        If (pochisq(chisqval, df) < p) Then
            maxchisq = chisqval
        Else
            minchisq = chisqval
        End If
        chisqval = (maxchisq + minchisq) * 0.5
    Loop
    critchi = chisqval
    
End Function

'renvoie, pour la probabilité p et le degré de liberté df, la variable
Function invchi2(ByVal p As Double, ByVal df As Double) As Double

    invchi2 = critchi(p, df)
    
End Function

'renvoie, pour la variable x et le degré de liberté df, la probabilité
Function chi2(ByVal x As Double, ByVal df As Integer) As Double

    chi2 = pochisq(x, df)
    
End Function

Conclusion

Pour la loi inverse de la loi normale centrée réduite, il s'agit d'un calcul de résolution d'une suite numérique convergente (méthode des tangentes de Newton). Je pense qu'il peut être optimisé.

Pour la loi du Khi-deux, j'avoue, ce n'est pas de moi (les info sur le code source sont en commentaires...), d'ailleurs, il est vraiment efficace. Il fonctionne à partir d'approximation polynomiales. Il fonctionne jusqu'à plus de 10000 degrés de liberté !!!

 

Historique

18 mai 2005 10:28:46 :
problème de mise en page des commentaires

Commentaires et avis

signaler à un administrateur
Commentaire de jack le 18/05/2005 18:50:06 administrateur CS

Salut
J'y comprends pas grand choses aux stats, mais côté programmation, plusieurs questions :
- A quoi peuvent donc servir les deux dernières lignes de ta fonction POZ ?
    Dim bigx As Double
    bigx = 20

- Le dimensionnement de tes variables (de toutes tes fonctions) est mal écrit : Quand on fait un Dim suivi d'une liste de variables, celles ci ne prennent pas le type précisé à la fin, mais le vulgaire Variant. Attention aux précisions des valeurs qui peuvent s'en retrouver altérées.
    Dim y, x, w As Double  --->  Dim y As Double, x As Double, w As Double

- Ce qui serait bien, serait de nous donner un petit Zip qui contiendrait un exemple concret de ces fonctions. Il doit bien y avoir une banalité en stat qui pourrait servir d'exemple, non ?

Sujet pas évident à traiter. (pas fréquent sur ce site, bravo)

signaler à un administrateur
Commentaire de marinmarais le 19/05/2005 09:20:13

Salut Jack,

Merci d'avoir pris le temps de regarder ce code.
Il est vrai que c'est un peu spécialisé, plus précisément, on en a pas souvent besoin. Mais quand celui-ci se ressent, on est bien embêté. Du coup, j'ai programmé une partie de ce module tout seul, pour l'autre, (n'étant pas assez compétent en mathématiques, ainsi que mes professeurs...) j'ai traduit un code source en JavaScript un peu à l'arrachée, d'où certaines imperfections :
- Les deux dernières lignes de la fonction POZ sont la définition d'un critère de rejet dans la fonction POCHISQ (cf ligne 154).
- Pour le dimensionnement de mes variables, je ne savais pas que je risquais de perdre de la précision. Je vais donc mettre ça à jour prochainement, j'en profiterai pour ajouter un exemple zipé.
- En attendant cet exemple plus parlant, je vais expliquer mes besoins :

J'ai programmé un soft servant à calculer les erreurs systématiques d'un appareil de topographie. En clair, la relation entre la valeur vraie et l'observation est une droite. Donc, le soft calcule une droite de régression par la méthode des moindres carrés.
Ces fonctions statistiques servent à vérifier si il n'y a pas de problèmes.
Ainsi, le test du Khi-deux consiste à tester le facteur unitaire de variance (en gros les écarts-types sur mes observations, donc leur précision) entre deux bornes calculées par la fonction INVCHI2. S'il est entre les deux bornes, c'est OK, si non, tu peux supposer qu'il y a un problème.
Dans le cas où le test précédent est un échec, le soft teste chaque résidus d'observation (l'écart entre la valeur mesurée et la valeur ajustée) par la loi normale centrée réduite. Là, par exemple, tu peux dire que le premier test a échoué parce que telle ou telles observations ne sont pas bonnes.

Ces bornes du Khi-deux sont très dure à calculer, en général, on les trouve dans des tableaux. Dans mon cas, c'est impossible car je suis obligé d'ajuster globalement plusieurs milliers d'observations : il n'a pas de tableaux assez grand ;-)

Voila, à bientôt pour l'exemple zipé...

Tom.

signaler à un administrateur
Commentaire de Alain Proviste le 19/05/2005 12:21:12 administrateur CS

jack fallait à l'école ptetre que t'aurais compris \o/

Ajouter un commentaire



Nos sponsors

Sondage...

CalendriCode

Janvier 2009
LMMJVSD
   1234
567891011
12131415161718
19202122232425
262728293031 

Consulter la suite du CalendriCode



Développement réalisé par Nicolas SOREL (Nix) avec l'aide de : Cyril DURAND et Emmanuel BAÏSE, Merci à Vincent pour ses précieux conseils
CodeS-SourceS.com© Toute reproduction même partielle est interdite sauf accord écrit du Webmaster
CodeS-SourceS.com© est une marque déposée tous droits réservés
Temps d'éxécution de la page : 0,218 sec

Google Coop CodeS-SourceS Google Coop CodeS-SourceS


Certaines images présentes sur le site (notament certains avatars) sont issues des collections IconShock, donc si vous souhaitez utiliser ces icons vous devez les acheter, ne les copiez pas et ne utilisez pas dans vos sites et applications sans les avoir commandé.