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 / Exemple :
' 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é !!!
Vous n'êtes pas encore membre ?
inscrivez-vous, c'est gratuit et ça prend moins d'une minute !
Les membres obtiennent plus de réponses que les utilisateurs anonymes.
Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.
Le fait d'être membre vous permet d'avoir des options supplémentaires.