''' <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
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.