[modification le 17/3/2014: fonction de répartition de la loi normale centrée réduite]

Article de référence ici: http://fr.wikipedia.org/wiki/Loi_normale

La loi normale est souvent appelée “courbe en cloche”:

Son équation dans le cas général est:

<m>f(x)~=~1/{sigma~sqrt{2pi}}e^{-1/2lbracem_f_u_1_sqrt_2pie^{-1/2{u}^2}</m>

La loi normale est souvent appelée la “loi des grands nombres”, parce qu'une variable qui dépend de très nombreux phénomènes aléatoires converge souvent vers une telle distribution .

Il faut calculer la surface qui se trouve à gauche de la valeur u donnée, et on aura trouvé la probabilité pour que la variable aléatoire se trouve inférieure à u.

Cette surface a pour formule mathématique:

<m>Prob(<u)~=~int{-infty}{u}{1/{sqrt{2pi}}e^{-1/2{u}^2}}du</m>

Pour calculer cette surface, on va faire quelque chose qui ressemble à l'intégrale: on va décomposer la surface en “tranches” verticales toute petites (2000 tranches par écart-type!). Pour chaque tranche, on peut en calculer la hauteur puisqu'on a l'équation de la courbe, et on peut donc en calculer la surface puisqu'on a la largeur. On fera la somme de tous ces petits rectangles et on aura la surface recherchée, donc la probabilité recherchée.

Comme la courbe est symétrique, on va aussi considérer que la surface correspondante à u<=0 est égale à 0.5, ce qui nous évitera de calculer jusqu'à - l'infini.

Ce qui fait que nous calculerons toujours la surface pour abs(u), et nous corrigerons si le u initial était négatif.

On va arrondir le résultat à 7 chiffres après la virgule: il est rare que les tables donnent plus de 5 chiffres après la virgule, et, compte tenu de la grande quantité de calculs, l'accumulation des petites erreurs de calcul liées à la précision des flottants, nous obligent à limiter la précision obtenue. Pour une plus grande précision, il faudrait faire les calculs avec le module decimal de Python.

from math import *
 
def pgaussred(x):
    """fonction de répartition de la loi normale centrée réduite
       (= probabilité qu'une variable aléatoire distribuée selon
       cette loi soit inférieure à x)
       calcul par intégration numérique: 2000 tranches par écart-type
       résultat arrondi à 7 chiffres après la virgule
    """
    if x==0:
        return 0.5
    u=abs(x)
    n=int(u*2000)
    du=u/n
    k=1/sqrt(2*pi)
    u1=0
    f1=k
    p=0.5
    for i in range(0,n):
        u2=u1+du
        f2=k*exp(-0.5*u2*u2)
        p=p+(f1+f2)*du*0.5
        u1=u2
        f1=f2
    if x<0:
        p = 1.0-p
    return round(p, 7)

Exemples d'utilisation

for x in range(-6, 7):
    P = pgaussred(x)
    print(x, P)

Ce qui donne:

-6 0.0
-5 3e-07
-4 3.17e-05
-3 0.0013499
-2 0.0227501
-1 0.1586553
0 0.5
1 0.8413447
2 0.9772499
3 0.9986501
4 0.9999683
5 0.9999997
6 1.0

Certains livres de mathématiques propose des formules approchées qui donnent de bons résultats plus rapidement. J'ai retenu une des formules proposées par l'excellent livre “Handbook of Mathematical Functions” de Abramovitz & Stegun. Voilà une proposition de codage Python:

from math import *
 
def pgaussred(x):
    """fonction de répartition de la loi normale centrée réduite
       (= probabilité qu'une variable aléatoire distribuée selon
       cette loi soit inférieure à x)
       formule simplifiée proposée par Abramovitz & Stegun dans le livre
       "Handbook of Mathematical Functions" (erreur < 7.5e-8)
    """
    u = abs(x) # car la formule n'est valable que pour x>=0
 
    Z = 1/(sqrt(2*pi))*exp(-u*u/2) # ordonnée de la LNCR pour l'absisse u
 
    b1 = 0.319381530
    b2 = -0.356563782
    b3 = 1.781477937
    b4 = -1.821255978
    b5 = 1.330274429
 
    t = 1/(1+0.2316419*u)
    t2 = t*t
    t4 = t2*t2
 
    P = 1-Z*(b1*t + b2*t2 + b3*t2*t + b4*t4 + b5*t4*t)
 
    if x<0:
        P = 1.0-P # traitement des valeurs x<0
 
    return round(P, 7) # retourne une valeur arrondie à 7 chiffres

Cette formule donne de bons résultats avec une erreur < 7.5e-8. Il est d'ailleurs très rare qu'une table donne plus que 5 chiffres après la virgule.

Exemple d'utilisation:

for x in range(-6, 7):
    P = pgaussred(x)
    print(x, P)

Ce qui donne:

-6 0.0
-5 3e-07
-4 3.17e-05
-3 0.00135
-2 0.0227501
-1 0.1586553
0 0.5
1 0.8413447
2 0.9772499
3 0.99865
4 0.9999683
5 0.9999997
6 1.0

Résultats quasi identiques au code précédent

Comme nous l'avons dit plus haut, avec le changement de variable <m>u~=~(x-m)/sigma</m>, on utilise les probabilités calculées sur la loi centrée réduite.

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division
 
from math import *
 
def pgauss(x,m,s):
    """pgauss(x): probabilité qu'une variable aléatoire distribuée selon une loi normale(m,s) soit inférieure à x"""
    return pgaussred((x-m)/s)
 
# exemples d'utilisation
print pgauss(2,1,2) # affiche: 0.691462457607
print pgauss(-2,1,2) # affiche: 0.0668072053163
print pgauss(0,0,1) # affiche: 0.5
print pgauss(7,1,2) # affiche: 0.998650101691
print pgauss(3,0,1) # affiche: 0.998650101691

C'est le contraire du cas précédent, et c'est souvent appelé dans les tables “papier”: “fractiles de la loi normale”: calculer le u correspondant à une probabilité P donnée que la variable soit inférieure à u.

On va aussi utiliser des petites tranches verticales, à part qu'on ne sait pas d'avance combien il y en aura. En fait, on va calculer par approche successive.

• on commence par calculer la surface par petites tranches du=1/2000 (=2000 tranches par écart-type).

• à un moment, la surface dépasse la valeur cherchée, et donc on revient en arrière avec des tranches plus petites (moitié).

• ceci jusqu'à ce qu'on soit assez près de la surface donnée, et on a trouvé le “u” cherché.

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division
 
from math import *
 
def xgaussred(p):
    """xgaussred(p): renvoie la variable x distribuée selon une LNR, et correspondant à une probabilité donnée p"""
    if p==0.5:
        return 0.0
    if p==0.0:
        return -7.56
    if p==1.0:
        return +7.56
    if p>0.5:
        pc=p
    else:
        pc=1-p
    du=1/2000
    eps=1.0E-15
    k=1/sqrt(2*pi)
    u1=0.0
    f1=k
    p1=0.5
    while True:
        u2=u1+du
        f2=k*exp(-0.5*u2*u2)
        s=(f1+f2)*0.5*du
        p2=p1+s
        if abs(p2-pc)<eps:
            break
        if ((p1<pc) and (p2<p1)) or  ((p1>pc) and (p2>p1)):
            # =on est en train de s'éloigner du point recherché
            du=-du/2
        u1=u2
        f1=f2
        p1=p2
    if p>0.5:
        return u2
    else:
        return -u2
 
# exemples d'utilisation:
xgaussred(0.97725)  => affiche (arrondi à 5 chiffres après la virgule): 2.00000
xgaussred(0.998650101691) => affiche (arrondi à 5 chiffres après la virgule): 3.00000

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division
 
from math import *
 
def xgauss(p,m,s):
    """xgauss(p,m,s): renvoie la variable x distribuée selon une LN(m,s), et correspondant à une probabilité donnée p"""
    return m+s*xgaussred(p)
 
# exemples d'utilisation:
xgauss(0.97725,0,1)  => affiche (arrondi à 5 chiffres après la virgule): 2.00000
xgauss(0.998650101691,0,1) => affiche (arrondi à 5 chiffres après la virgule): 3.00000

On applique un principe simple: on tire au hasard un nombre compris entre 0 et 1 avec random() du module random, et on considère cette valeur comme une probabilité gaussienne dont on veut avoir le u correspondant. On va donc naturellement utiliser la fonction précédente.

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division
 
from math import *
 
def hgaussred(n=0):
    """hgaussred(): renvoie une valeur au hasard de la variable x distribuée selon une LN réduite"""
    if n==0:
        return xgaussred(random())
    else:
        L=[]
        for i in range(0,n):
            L.append(xgaussred(random()))
        return L
 
# exemples d'utilisation:
hgaussred()  # affiche une valeur au hasard entre -infini et +infini distribuée selon une loi normale réduite comme 0.851116409419517
hgaussred(5)  # affiche une liste de 5 valeurs au hasard entre - l'infini et + l'infini distribuée selon une loi normale réduite comme [1.60380862369612,1.47044806011366,-1.33672988492718,-1.07610674009864,-0.771293171367761]

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division
 
from math import *
 
def hgauss(m,s,n=1):
    """hgauss(m,s): renvoie une valeur au hasard de la variable x distribuée selon une LN(m,s)"""
    if n==1:
        return xgauss(random(),m,s)
    else:
        L=[]
        for i in range(0,n):
            L.append(xgauss(random(),m,s))
        return L
 
# exemples d'utilisation:
hgauss(2,5)  # affiche une valeur au hasard entre -infini et +infini distribuée selon une loi normale m = 2, s = 5 comme 4.5505296649258
hgauss(2,5,5)  # affiche une liste de 5 valeurs quelconques entre - l'infini et + l'infini distribuée selon une loi normale m = 2, s = 5 comme [-2.94898945875669,-0.706714534003603,6.31447343047467,5.77588811971859,1.20454578795292]