Calcul de la loi normale et de la fonction d'erreur de Gauss

Commenter

Contenu du snippet 

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

Compatibilité : VB 2005, VB 2008, VB.NET 1.x

A voir également

S'inscrire maintenant

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.

Du même auteur (circular)