La “loi de Poisson” est une loi de probabilité associée à une variable aléatoire discrète (=ne prenant que des valeurs entières).

Voir ici pour la référence: http://fr.wikipedia.org/wiki/Loi_de_Poisson

Elle est souvent utilisée comme la loi des évènements rares. Par exemple:

• Une entreprise a 2 accidents du travail par mois. La loi de Poisson donne la probabilité qu'il y ait 0 , 1, 2, 3, etc… accident(s) dans un mois donné.

• La fabrication d'un tissus donne un défaut tous les 20 mètres. La loi de Poisson permet de calculer la probabilité qu'il y ait 0, 1, 2, 3, … défaut(s) dans une pièce de 20 mètres donnée

• Pendant un orage, il y a 3 éclairs par minutes en moyenne. La loi de Poisson permet de calculer la probabilité qu'il y ait 0, 1, 2, 3, 4, etc… éclair(s) dans une minute donnée.

La probabilité est donnée par la formule suivante:

<m>P(k)~=~{e^-m}*m_k_k</m>

Ou, en langage Python:

p=e**(-m)*m**k/fact(k)

Avec:

• e = 2.718281828459045

• m = paramètre caractéristique: c'est en même temps la moyenne et la variance de la loi.

• fact(k) = factorielle de k

1ère solution de codage:

On applique brutalement la formule de définition, mais il faut la fonction factorielle qui n'existe pas en Python de base.

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division
 
from math import *
 
def fact(n):
    """fact(n): calcule la factorielle de n (entier >= 0)"""
    x=1
    for i in xrange(2,n+1):
        x*=i
    return x
 
def poisson(k,m):
    """poisson(k,m): donne la probabilité d'avoir k évènements distribués selon une loi de Poisson de paramètre m"""
    return e**(-m)*m**k/fact(k)
 
# exemple d'utilisation
print poisson(2,3)  # affiche: 0.224041807655

2ème solution de codage:

On élimine en même temps le calcul de la puissance et le calcul de la factorielle pour éviter les dépassements de capacité.

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division
 
from math import *
 
def poisson(k,m):
    """poisson(k,m): donne la probabilité d'avoir k évènements distribués selon une loi de Poisson de paramètre m"""
    p=e**(-m)
    for i in xrange(0,k):
        p*=m/k
        k-=1
    return p
 
# exemple d'utilisation
print poisson(2,3)  # affiche: 0.224041807655

On va maintenant calculer la probabilité qu'il y ait 0 ou 1 ou 2 ou … k évènements dans une distribution de Poisson de paramètre m.

1ère solution de codage:

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division
 
from math import *
 
def poissoncum(k,m):
    """poissoncum(k,m): probabilité d'avoir k évènement(s) ou moins, dans une distribution de Poisson de paramètre m."""
    p=0
    for i in xrange(0,k+1):
        p+=poisson(i,m)
    return p
 
# exemple d'utilisation
print poissoncum(2,3)  # affiche: 0.423190081127 (idem à "poisson(0,3) + poisson(1,3) + poisson(2,3)")
print poissoncum(0,3)  # affiche: 0.0497870683679 (idem à "poisson(0,3)") 
print poissoncum(50,3)  # affiche: 1.0

2ème solution de codage:

Au lieu d'un calcul brutal, on s'appuie sur les 2 formules suivantes:

• <m>P(0)~=~e^{-m}</m>

• <m>P(k)~=~{P(k-1)}*{m/k}</m>

Ce qui permet de diminuer fortement la quantité de calcul:

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division
 
from math import *
 
def poissoncum(k,m):
    """poissoncum(k,m): probabilité d'avoir k évènement(s) ou moins, dans une distribution de Poisson de paramètre m."""
    pk=e**(-m)
    p=pk
    for i in xrange(1,k+1):
        pk*=m/i
        p+=pk
    return p
 
# exemple d'utilisation
print poissoncum(2,3)  # affiche: 0.423190081127 (idem à "poisson(0,3) + poisson(1,3) + poisson(2,3)")
print poissoncum(0,3)  # affiche: 0.0497870683679 (idem à "poisson(0,3)") 
print poissoncum(50,3)  # affiche: 1.0

Puisque dans la vie réelle, les prélèvements sont souvent faits pour faire des estimations et prendre des décisions, il est important de savoir quels sont les risques qu'on a de se tromper!

Problème: Sur un site industriel, on mesure qu'il y a eu k accidents en 1 semaine, mais ce nombre k est lié au hasard, et on voudrait évaluer la vrai moyenne des accidents correspondant aux risques habituels: que peut-on en dire?

Solution: En limitant le risque d'erreur à <m>alpha~=~5%</m>, la moyenne des accidents par semaine sera comprise entre m1 et m2 calculés de cette façon:

• m1 est tel que: <m>sum{i=k}{infty}{e^-m1}*m1_i_i~=~alpha/2</m> ou en utilisant notre fonction “loi de Poisson cumulée”: 1-poissoncum(k-1,m1) = <m>alpha/2</m>

• m2 est tel que: <m>sum{i=0}{k}{e^-m2}*m2_i_i~=~alpha/2</m> ou en utilisant notre fonction “loi de Poisson cumulée”: poissoncum(k,m2) = <m>alpha/2</m>

Le risque <m>alpha~=~5%</m> signifie que dans 2.5% la vraie moyenne sera inférieur à m1 et dans 2.5% des cas, elle sera supérieure à m2.

Donc, dans <m>(1-alpha)=95%</m> des cas (c'est à dire 19 fois sur 20), nous aurons raison en disant que la vraie moyenne est entre m1 et m2.

Compte tenu du type de formules, on ne peut pas les calculer simplement comme m=f(k,<m>alpha</m>), mais comme f(m,k,<m>alpha</m>)=0. On fait donc une recherche de la valeur de m1 et de m2 qui annule leur fonction, et on le fait par approche successive.

La fonction poissonconfsup(k,conf) calcule m2, la borne supérieure, et poissonconfinf(k,conf) calcule m1, la borne inférieure. Dans chacun des 2 cas, l'argument conf est ici la probabilité que la valeur de m ne dépasse pas la borne. Dans l'exemple ci-dessus, ce serait 0.975, puisque un risque de 2.5% qu'on soit inférieur à m1 et 2.5% qu'on soit supérieur à m2, ça fait bien une confiance de 95% d'être entre m1 et m2, et un risque de 5% d'être en dehors.

La fonction poissonconf(k,n,conf) renvoie la liste des 2 bornes de l'intervalle de confiance m1 et m2 avec la confiance conf (ici 0.95) d'être entre les 2.

Voilà le codage proposé:

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division
 
from math import *
 
def poissonconfsup(k,conf):
    """poissonconfsup(k,conf): borne supérieure sur la moyenne trouvée avec un échantillon (k) et une confiance conf """
    dp=0.1
    eps=1.0e-8
    m1=0.5
    r=1-conf
    f1=poissoncum(k,m1)-r
    while True:
        m2=m1+dp
        f2=poissoncum(k,m2)-r
        if abs(f2)<eps:
            break
        if (abs(f2)>abs(f1)) or (f1*f2<0):
            dp=-dp/2
        m1=m2
        f1=f2
    return m2
 
# exemple d'utilisation
print poissonconfinf(1,0.975) # affiche: 0.0253178119659
print poissonconfinf(2,0.975) # affiche: 0.242209243774
print poissonconfinf(3,0.975) # affiche: 0.618672180176

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division
 
from math import *
 
def poissonconfinf(k,conf):
    """poissonconfinf(k,conf): borne inférieure sur la moyenne trouvée avec un échantillon (k) et une confiance conf """
    dp=0.1
    eps=1.0e-8
    m1=0.5
    r=1-conf
    f1=1-poissoncum(k-1,m1)-r
    while True:
        m2=m1+dp
        f2=1-poissoncum(k-1,m2)-r
        if abs(f2)<eps:
            break
        if (abs(f2)>abs(f1)) or (f1*f2<0):
            dp=-dp/2
        m1=m2
        f1=f2
    return m2
 
# exemple d'utilisation
print poissonconfsup(1,0.975) # affiche: 5.57164306641
print poissonconfsup(2,0.975) # affiche: 7.22468719482
print poissonconfsup(3,0.975) # affiche: 8.76727294922

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division
 
from math import *
 
def poissonconf(k,conf):
    """poissonconf(k,n,conf): intervalle de confiance sur la moyenne trouvée avec un échantillon (k) et une confiance conf """
    bi=poissonconfinf(k,conf+(1-conf)/2)
    bs=poissonconfsup(k,conf+(1-conf)/2)
    return [bi,bs]
 
# exemple d'utilisation
print poissonconf(1,0.95) # affiche: [0.025317811965942395, 5.5716430664062466]
print poissonconf(2,0.95) # affiche: [0.24220924377441408, 7.2246871948242113]
print poissonconf(3,0.95) # affiche: [0.61867218017578118, 8.7672729492187322]

Application pratique:

Sur un site industriel et sur une assez longue période stable, j'ai calculé une moyenne de 2 accidents du travail par semaine. Alors je mène des actions pour les réduire, et après la mise en place de ces actions, je mesure 1 seul accident pendant une semaine donnée: est-ce que ça démontre une baisse effective des risques?

Calculons la borne supérieure m2 avec un risque de 95%:

m2=poissonconfsup(1,0.95) ce qui affiche 4.74386444092

Donc, je peux dire qu'avec 1 accident pour une semaine donnée, la vraie moyenne d'accident est inférieur à 4.7 accidents par semaine: je ne sais donc pas si les actions que j'ai menées ont eu un effet.

En fait, cette période d'observation est trop courte. Imaginons que sur les 8 semaines qui suivent les actions, je n'ai observé que 8 accidents, donc toujours une moyenne de 1 accident par semaine:

m2=poissonconfsup(8,0.95)/8 ce qui affiche 1.80433120728

Là, j'en suis sûr (à 95%): les actions que j'ai menées ont effectivement fait baissé les accidents!

Cette fonction peut fournir une seule valeur ou une liste de valeurs (longueur de la liste donnée par le paramètre optionnel nb).

On applique le principe général de ce genre de calcul:

• on tire au hasard un nombre décimal compris entre 0 et 1 avec random(), et on le considère comme une probabilité cumulée.

• grâce au calcul de la probabilité cumulée de la loi de Poisson, on cherche à quel variable aléatoire k cela correspond

Pour alléger les calculs, on utilise ici le calcul de la probabilité normale et pas celle cumulée. En fait, on la recalcule dans la recherche.

Codage proposé:

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division
 
from math import *
import random
 
def hpoisson(m,nb=0):
    """hpoisson(m,nb=0): Génération de valeurs tirées au hasard selon une distribution de Poisson"""
    def _hpoisson(m):
        ph=random.random()
        k=0
        pc=poisson(k,m)
        while ph>pc:
            k+=1
            pc+=poisson(k,m)
        return k
 
    if nb==0:
        return _hpoisson(m)
    else:
        R=[]
        for j in xrange(0,nb):
            R.append(_hpoisson(m))
        return R
 
# exemple d'utilisation:
print hpoisson(2)   # affiche par exemple: 3 = une valeur qui peut aller de k=0 à l'infini  
print hpoisson(2)   # idem
print hpoisson(2)   # idem
print hpoisson(2,10)  # affiche par exemple: [1, 1, 1, 0, 3, 2, 3, 2, 0, 2] 

Avec cette fonction, on peut faire des simulations de tirage au hasard.