''' <summary>
''' Loi normale et fonction d'erreur
''' </summary>
Public Class NormalDistribution
Shared rel_error As Double = 0.000000000001
''' <summary>
''' Fonction d'erreur de Gauss
''' </summary>
Shared Function erf(ByVal x As Double) As Double
'erf(x) = 2/sqrt(pi)*integral(exp(-t^2),t,0,x)
' = 2/sqrt(pi)*[x - x^3/3 + x^5/5*2! - x^7/7*3! + ...]
' = 1-erfc(x)
Const two_sqrtpi As Double = 1.1283791670955126
If (Math.Abs(x) > 2.2) Then
Return 1.0 - erfc(x) 'use continued fraction when |x| > 2.2
End If
Dim sum As Double = x, term As Double = x, xsqr As Double = x * x
Dim j As Integer = 1
Do
term *= xsqr / j
sum -= term / (2 * j + 1)
j += 1
term *= xsqr / j
sum += term / (2 * j + 1)
j += 1
Loop While Math.Abs(term / sum) > rel_error
Return two_sqrtpi * sum
End Function
''' <summary>
''' Complémentaire à 1 de la fonction d'erreur
''' </summary>
Shared Function erfc(ByVal x As Double) As Double
'erfc(x) = 2/sqrt(pi)*integral(exp(-t^2),t,x,inf)
' = exp(-x^2)/sqrt(pi) * [1/x+ (1/2)/x+ (2/2)/x+ (3/2)/x+ (4/2)/x+ ...]
' = 1-erf(x)
'expression inside [] is a continued fraction so '+' means add to denominator only
Dim one_sqrtpi As Double = 0.56418958354775628
If Math.Abs(x) < 2.2 Then
Return 1.0 - erf(x) 'use series when |x| < 2.2
End If
If x <= 0 Then 'continued fraction only valid for x>0
Return 2.0 - erfc(-x)
End If
Dim a As Double = 1, b As Double = x 'last two convergent numerators
Dim c As Double = x, d As Double = x * x + 0.5 'last two convergent denominators
Dim q1 As Double, q2 As Double = b / d 'last two convergents (a/c and b/d)
Dim n As Double = 1.0, t As Double
Do
t = a * n + b * x
a = b
b = t
t = c * n + d * x
c = d
d = t
n += 0.5
q1 = q2
q2 = b / d
Loop While Math.Abs(q1 - q2) / q2 > rel_error
Return one_sqrtpi * Math.Exp(-x * x) * q2
End Function
''' <summary>
''' Logarithme néperien de la fonction complémentaire de la fonction d'erreur
''' </summary>
Shared Function logerfc(ByVal x As Double) As Double
'erfc(x) = 2/sqrt(pi)*integral(exp(-t^2),t,x,inf)
' = exp(-x^2)/sqrt(pi) * [1/x+ (1/2)/x+ (2/2)/x+ (3/2)/x+ (4/2)/x+ ...]
' = 1-erf(x)
'expression inside [] is a continued fraction so '+' means add to denominator only
Dim one_sqrtpi As Double = 0.56418958354775628
If Math.Abs(x) < 2.2 Then
Return Log(1.0 - erf(x)) 'use series when |x| < 2.2
End If
If x <= 0 Then 'continued fraction only valid for x>0
Return Log(2.0 - erfc(-x))
End If
Dim a As Double = 1, b As Double = x 'last two convergent numerators
Dim c As Double = x, d As Double = x * x + 0.5 'last two convergent denominators
Dim q1 As Double, q2 As Double = b / d 'last two convergents (a/c and b/d)
Dim n As Double = 1.0, t As Double
Do
t = a * n + b * x
a = b
b = t
t = c * n + d * x
c = d
d = t
n += 0.5
q1 = q2
q2 = b / d
Loop While Math.Abs(q1 - q2) / q2 > rel_error
Return Log(one_sqrtpi * q2) - x * x
End Function
''' <summary>
''' Fonction mathématique : produit de erfc et de exp
''' </summary>
''' <param name="y">Argument pour la fonction erfc</param>
''' <param name="argExp">Argument pour la fonction exponentielle</param>
Shared Function erfexp(ByVal y As Double, ByVal argExp As Double) As Double
If y <= 0 Or argExp <= 0 Then
Return NormalDistribution.erfc(y) * Exp(argExp)
End If
Dim ArgAux As Double = logerfc(y) + argExp
Return Exp(argAux)
End Function
''' <summary>
''' Fonction de répartition de la loi normale, c'est-à-dire
''' la probabilité P(X <= z)
''' </summary>
Shared Function Phi(ByVal z As Double)
Return 0.5 * (1 + erf(z / Sqrt(2)))
End Function
End Class