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 considérer aussi que pour u>=7.56 (7.56 x écarts-types), la probabilité est égale à 1. On peut en effet vérifier que les limites du calcul en double précision sont atteintes au delà. De même pour u<=-7.56, valeur pour laquelle la fonction renverra une probabilité nulle.

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division
 
from math import *
 
def pgaussred(x):
    """pgaussred(x): probabilité qu'une variable aléatoire distribuée selon une loi normale réduite soit inférieure à x"""
    if x==0:
        return 0.5
    if x>=7.56:
        return 1.0
    if x<=-7.56:
        return 0.0
    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 xrange(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:
        return p
    else:
        return 1.0-p
 
# exemples d'utilisation
print pgaussred(-8) # affiche: 0.0 
print pgaussred(-7) # affiche: 1.26010313295e-012
print pgaussred(-6) # affiche: 9.86570047878e-010
print pgaussred(-5) # affiche: 2.86651707593e-007
print pgaussred(-4) # affiche: 3.1671252966e-005
print pgaussred(-3) # affiche: 0.0013498983086
print pgaussred(-2) # affiche: 0.0227501341978
print pgaussred(-1) # affiche: 0.158655258973
print
print pgaussred(0) # affiche: 0.5
print
print pgaussred(1) # affiche: 0.841344741027
print pgaussred(2) # affiche: 0.977249865802
print pgaussred(3) # affiche: 0.998650101691
print pgaussred(4) # affiche: 0.999968328747
print pgaussred(5) # affiche: 0.999999713348
print pgaussred(6) # affiche: 0.999999999013
print pgaussred(7) # affiche: 0.999999999999
print pgaussred(8) # affiche: 1.0

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]