(Modifié le 24/4/2008 pour élargir la plage d'application.)

Référence pour la définition: voir http://fr.wikipedia.org/wiki/Loi_binomiale

Prenons un sac contenant beaucoup de billes (ce sac est appelé la “population”). La plupart des billes sont blanches, mais il y en a qui sont vertes et qui représentent p% de l'ensemble du sac.

On tire du sac un échantillon de n billes. On se demande quelle est la probabilité d'avoir k billes vertes dans cet échantillon.

On considère ici que l'échantillon est petit par rapport à la quantité de billes du sac (disons <10%), ou alors qu'on fait le prélèvement de l'échantillon bille par bille et qu'on remet chaque bille prélevée dans le sac avant le tirage suivant. Dans tous les cas, il faut que chaque prélèvement d'une bille ne change pas de façon significative la probabilité pour les billes suivantes.

Pour les fans, cette expérience est attribuée à Bernouilli (tirage dans une urne de Bernouilli)

La formule donnant la probabilité cherchée est la suivante:

<m>p(k)~=~{{C_n}^k}p^k(1-p)^{n-k}</m> Avec <m>{{C_n}^k}</m> = combinaison de n objets pris k à k (voir la page du site sur l'analyse combinatoire) **1ère solution de codage** Pour le code, nous utiliserons la fonction "combin(n,k)" calculant la combinaison de n objets pris k à k déjà présentée sur ce site ([[http://python.jpvweb.com/mesrecettespython/combinatoire]]). <code python> #!/usr/bin/python # -*- coding: utf-8 -*- from __future__ import division def binom(k,n,p): """binom(k,n,p): probabilité d'avoir k réussite(s) dans n évènements indépendants, chaque évènement ayant une probabilité p% de réussite""" x=combin(n,k)*pow(p,k)*pow(1-p,n-k) return x # Exemple d'utilisation: print binom(2,10,0.25) # => affiche 0.281567573547363 </code> **2ème solution de codage** Dans certains cas extrêmes, le code précédent génère une erreur de dépassement de capacité dans les calculs intermédiaires. Ce qu'il n'est pas facile ici, c'est de rendre le code capable de calculer avec des valeurs importantes comme k=1000, n=3000 sans plantage et sans que ça dure trop longtemps. Vous voyez dans le code suivant qu'au fur et à mesure que le calcul se déroule (boucle for), la multiplication par n/k fait augmenter, et la multiplication par p ou par q fait diminuer. On utilise donc une astuce (les 2 boucles while) pour que la valeur calculée à chaque boucle for ne dépasse jamais les capacités de calcul en nombres réels, et ceci, quelque soit le k et le n donné. Code proposé: <code python> #!/usr/bin/python # -*- coding: utf-8 -*- from __future__ import division def binom(k,n,p): """binom(k,n,p): probabilité d'avoir k réussite(s) dans n évènements indépendants, chaque évènement ayant une probabilité p% de réussite""" if k==n: return p**n if p==0: return 0.0 if p==1.0: # NB: le cas k==n a déjà été traité plus haut return 0.0 q=1-p if k==0: return q**n j1=k j2=n-k xmin=1.0e-200 x=1.0 for i in xrange(0,k): x*=n/k while (x>xmin) and (j1>0): x*=p j1-=1 while (x>xmin) and (j2>0): x*=q j2-=1 n-=1 k-=1 return x*(p**j1)*(q**j2) # Exemple d'utilisation: print binom(2,10,0.25) # => affiche 0.281567573547363 </code> L'exemple donné signifie: dans un échantillon de 10 billes, j'ai une probabilité de 28% d'avoir 2 billes vertes, sachant qu'il y en a 25% dans la population (=le sac). \\ Mais à quoi ça sert? Les occasions de prélever des billes dans un sac sont plutôt limitées dans la vie courante... En fait, cette expérience va servir de "modèle" pour d'autres expériences plus intéressantes qui lui ressemblent. Deux exemples très concrets: * Contrôle statistique: dans l'industrie, on veut prendre des décisions qui soient en même temps efficaces et peu couteuses, d'où les décisions sur prélèvement d'échantillon. On va donc prélever périodiquement des pièces (=l'échantillon), les mesurer et les qualifier "mauvaises" (=billes vertes) ou "bonnes" (=billes blanches). On va en déduire une décision: refuser le lot, arrêter la production, régler les machines, etc... * Sondage avant votes: Imaginez un sondage avant un référendum: on va consulter 1000 personnes (=l'échantillon) sur leur prévision de vote "oui" (=billes vertes) ou "non" (=billes blanche). On va en déduire une estimation du résultat. \\ Concernant la loi binomiale, on peut ajouter deux caractéristiques: * sa moyenne est égale à <m>M~=~np</m> ou, en Python, m=n*p * son écart-type est égal à <m>sigma~=~sqrt{np(1-p)}</m> ou, en Python, s=sqrt(n*p*(1-p)) ===== Loi binomiale cumulée ===== La loi binomiale cumulée est la probabilité d'obtenir 0 ou 1 ou 2 ... ou k billes vertes dans un échantillon de n billes, sachant qu'il en a p% dans la population. La formule est donc en notation mathématique: <m>sum{i=0}{k}{{C_n}^i}{{p}^i}(1-p)^{n-i}</m> **1ère solution de codage:** Le code Python est simple à en déduire en utilisant la fonction binomiale définie plus haut: binom(k,n,p) <code python> #!/usr/bin/python # -*- coding: utf-8 -*- from __future__ import division def binomcum(k,n,p): """binomcum(k,n,p): probabilité d'avoir 0,1,..,k réussite(s) dans n évènements indépendants, chaque évènement ayant une probabilité p% de réussite""" x=0 for i in xrange(0,k+1): x+=binom(i,n,p) return x # Exemple d'utilisation: print binomcum(2,10,0.25) # => affiche 0.525592803955078 </code> **2ème solution de codage:** La solution ci-dessus n'est guère astucieuse, parce qu'à chaque nouvelle valeur de binom, on recalcule tout à partir de zéro. Mais, comme pour la fonction binom(), ce qu'il n'est pas facile ici, c'est de rendre le code capable de calculer des valeurs importantes comme k=1000, n=3000 sans plantage et sans que ça dure trop longtemps. L'astuce utilisée ici est la suivante: les probabilités binomiales seront réparties autour de int(n*p). Mais pour des valeurs importantes de k et de n, ces probabilités laisseront avant (jusqu'à 0) et après (jusqu'à n) une plage importante avec des valeurs faibles ou nulles. Il s'agit ici, pour gagner du temps, de s'arrêter de calculer lorsqu'on a atteint ces valeurs. Code proposé: <code python> #!/usr/bin/python # -*- coding: utf-8 -*- from __future__ import division def binomcum(k,n,p): """binomcum(k,n,p): probabilité d'avoir 0,1,..,k réussite(s) dans n évènements indépendants, chaque évènement ayant une probabilité p% de réussite""" if k==n: return 1.0 if p==0: return 0.0 if p==1.0: # NB: le cas k==n a déjà été traité plus haut return 0.0 j=int(p*n) b=0 xmin=1.0e-20 if k<=j: for i in xrange(k,-1,-1): x=binom(i,n,p) b+=x if x<xmin: break else: for i in xrange(k+1,n+1): x=binom(i,n,p) b+=x if x<xmin: break b=1.0-b return b # Exemple d'utilisation: print binomcum(2,10,0.25) # => affiche 0.525592803955078 </code> ===== Intervalles de confiance sur la loi binomiale ===== Puisqu'on a dit que, dans la vie réelle, les prélèvements étaient faits pour estimer des proportions et prendre des décisions, il est important de savoir quels sont les risques qu'on a de se tromper! Problème: On prélève un échantillon de n billes et on en trouve k billes verte. Mais on ne sait pas quel pourcentage de billes vertes il y a dans le sac: que peut-on en dire? Solution: En limitant le risque d'erreur à <m>alpha~=~5%</m>, le pourcentage p de billes vertes dans la population sera compris entre pinf et psup calculés de cette façon: * pinf est tel que: <m>sum{i=k}{n}{{C_n}^i}{{p_inf}^i}(1-p_inf)^{n-i}~=~alpha/2</m> Ou en utilisant notre fonction "loi binomiale cumulée": 1-binomcum(k-1,n,pinf) = <m>alpha</m>/2 * psup est tel que: <m>sum{i=0}{k}{{C_n}^i}{{p_sup}^i}(1-p_sup)^{n-i}~=~alpha/2</m> Ou en utilisant notre fonction "loi binomiale cumulée": binomcum(k,n,psup) = <m>alpha</m>/2 Le risque <m>alpha~=~5%</m> signifie que dans 2.5% le vrai pourcentage sera inférieur à pinf et dans 2.5% des cas, il sera supérieur à psup. Donc, dans (<m>1-alpha</m>)=95% des cas (c'est à dire 19 fois sur 20), nous aurons raison en disant que le vrai pourcentage est entre pinf et psup. Compte tenu du type de formules, on ne peut pas les calculer simplement comme p=f(), mais comme f(p)=0. On fait donc une recherche de la valeur de pinf et de psup qui annule leur fonction, et on le fait par approche successive. On rencontre ce type de calcul (trouver x tel que f(x)=0) assez souvent, et il est intéressant d'en conserver un "modèle": - on fixe une valeur de départ p1, l'incrément de départ dp et on calcule une première valeur f1(p1) - on calcule p2=p1+dp et f2(p2) - si f2 est très proche de zéro (à eps près), on a trouvé! et on renvoie p2 - si f1 et f2 sont de même signe et f2 est plus grand que f1 => on s'éloigne du point cherché => revenir en arrière avec un incrément moitié - si f1 et f2 ne sont pas de même signe, on a dépassé le point zéro => revenir en arrière avec un incrément moitié - les nouvelles valeurs de départ sont p1=p2 et f1=f2 et on recommence! La fonction binomconfsup(k,n,conf) calcule psup, la borne supérieure, et binomconfinf(k,n,conf) calcule pinf, la borne inférieure. L'argument conf est ici la probabilité que la valeur de p 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 à pinf et 2.5% qu'on soit supérieur à psup, ça fait bien une confiance de 95% d'être entre pinf et psup, et un risque de 5% d'être en dehors. La fonction binomconf(k,n,conf) renvoie la liste des 2 bornes pinf et psup avec la confiance conf d'être entre les 2. \\ <code python> #!/usr/bin/python # -*- coding: utf-8 -*- from __future__ import division def binomconfsup(k,n,conf): """binomconfsup(k,n,conf): borne supérieure de confiance sur la proportion donnée par un échantillon (k,n) avec une confiance conf""" eps=1.0e-8 r=1.0-conf p1=k/n dp=0.1 f1=binomcum(k,n,p1)-r while True: p2=p1+dp f2=binomcum(k,n,p2)-r if abs(f2)<eps: break if ((f1>0)and(f2>0))or((f1<0)and(f2<0)): if abs(f2)>abs(f1): dp*=-0.5 else: dp*=-0.5 p1=p2 f1=f2 return p2 def binomconfinf(k,n,conf): """binomconfinf(k,n,conf): borne inférieure de confiance sur la proportion donnée par un échantillon (k,n) avec une confiance conf""" eps=1.0e-8 r=1-conf dp=-0.1 p1=k/n f1=1-binomcum(k-1,n,p1)-r while True: p2=p1+dp f2=1-binomcum(k-1,n,p2)-r if abs(f2)<eps: break if ((f1>0)and(f2>0))or((f1<0)and(f2<0)): if abs(f2)>abs(f1): dp*=-0.5 else: dp*=-0.5 p1=p2 f1=f2 return p2 def binomconf(k,n,conf): """binomconf(k,n,conf): intervalle de confiance sur le pourcentage trouvé avec un échantillon (k,n) et une confiance conf """ bi=binomconfinf(k,n,conf+(1-conf)/2) bs=binomconfsup(k,n,conf+(1-conf)/2) return [bi,bs] # Exemple d'utilisation: print binomconfinf(4,20,0.975) # affiche: 0.057333993911743181 print binomconfsup(4,20,0.975) # affiche: 0.43661398887634278 print binomconf(4,20,0.95) # affiche: [0.057333993911743181, 0.43661398887634278] print binomconf(510,1000,0.95) # affiche: [0.478525838851929,0.541415128707886] print binomconf(1000,2000,0.95) # affiche: [0.477850553393364,0.522149446606636] print binomconf(1020,2000,0.95) # affiche: [0.487840253710747,0.532130388021469] </code> Exemple pratique: pour évaluer le résultat des votes à un référendum, on questionne 1000 personnes, et on trouve: 510 réponses "oui" et 490 réponses "non". On a envie de dire que le résultat estimé est 51%-49%. Mais calculons l'intervalle de confiance sur la proportion. binomconf(510,1000,0.95) => [0.478525838851929,0.541415128707886] Ce qui veut dire que la vrai estimation sera entre 48% et 54%, et que donc on ne peut rien en déduire. Avec un échantillon plus grand, disons 2000 personnes, avec le même résultat 51%-49%: binomconf(1020,2000,0.95) => [0.487840253710747,0.532130388021469] On a encore une estimation entre 49% et 53% avec la même conséquence: la fourchette est plus petite, mais on ne sait toujours pas quel réponse va gagner. Vous commencez à comprendre pourquoi, lorsqu'on est aussi près de 50%, les journaux peuvent annoncer des résultats contradictoires... Si vous essayez ces derniers calculs, vous verrez que binomconf(1020,2000,0.95) demande déjà beaucoup de temps (près de 4 minutes chez moi). Au delà, il faut calculer les intervalles de confiance par d'autres méthodes qui ne sont pas (encore) décrites ici. Par exemple,on peut remplacer la loi binomiale par une loi normale de moyenne np et de variance np(1-p) dès que np(1-p)>18, l'erreur étant d'autant plus faible que n est grand et que p est proche de 0.5. Par contre, si p est faible (<10%) et n est grand, la loi binomiale peut être remplacée par une loi de Poisson de paramètre np, l'erreur étant d'autant plus faible que n est grand et p est faible. ===== Génération de valeurs tirées au hasard selon une distribution binomiale ===== 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é binomiale cumulée. * grâce au calcul de la probabilité binomiale cumulée, on cherche à quel variable aléatoire cela correspond (ici = k) Pour alléger les calculs, on utilise ici le calcul de la probabilité binomiale et pas celle cumulée. En fait, on la recalcule dans la recherche. <code python> #!/usr/bin/python # -*- coding: utf-8 -*- from __future__ import division import random def hbinom(n,p,nb=0): """hbinom(n,p,nb=0): Génération de valeurs tirées au hasard selon une distribution binomiale""" def _hbinom(n,p): ph=random.random() i=0 pc=binom(i,n,p) while ph>pc: i+=1 pc+=binom(i,n,p) return i if nb==0: return _hbinom(n,p) else: R=[] for j in xrange(0,nb): R.append(_hbinom(n,p)) return R # exemple d'utilisation print hbinom(20,0.3) # affiche par exemple: 7 = une valeur qui peut aller de k=0 à k=20 print hbinom(20,0.3,10) # affiche par exemple: [5, 7, 6, 10, 5, 8, 5, 7, 6, 9] </code> Avec cette fonction, on peut faire des simulations de tirage au hasard. <html> <head> <style type="text/css"> <!-- body {background-image:url(fondcorps.jpg);} --> </style> </head> <body> </body> </html>