Distribution binomiale
1. Introduction
On considère une variable aléatoire x pouvant prendre la valeur 1 avec une probabilité p, ou la valeur 0 avec une probabilité 1-p. L'espérance de cette variable est E(x)=p et sa variance var(x)=p(1-p).
Dans un premier temps, on fera une simulation de Monte-Carlo de cette variable aléatoire, pour calculer la moyenne et la variance empiriques. On étudiera aussi la distribution de la moyenne, que l'on comparera à une distribution normale.
La distribution binomiale donne la probabilité d'obtenir n fois la valeur 1 lorsqu'on fait N tirages. On fera des simulations de Monte-Carlo pour obtenir cette distribution, pour différentes valeurs de N et p.
Pour générer des nombres pseudo-aléatoires, on utilisera la fonction random.random().
2. Moyenne et variance
2.a. Fonctions de calcul
[1] Écrire une fonction qui génère des échantillons de la variable aléatoire x, en utilisant la fonction random.random().
[2] Écrire une fonction qui calcule, pour N tirages, la moyenne, l'écart-type, et une estimation de l'écart-type de la moyenne. On rappelle que la variance de la moyenne est égale à la variance de la variable x divisée par N.
[3] Tester cette fonction pour différentes valeurs de p, et pour des valeurs de N croissantes.
2.b. Distribution des valeurs moyennes
Considérons Nt calculs de la valeur moyenne, chacun étant fait avec N tirages. Notons mi ces valeurs moyennes. Pour obtenir la distribution de ces valeurs moyennes, on construit un histogramme. L'intervalle de définition de la variable aléatoire, ici [0,1], est subdivisé en M classes, c'est-à-dire M sous-intervalles (de largeur 1/M). L'histogramme est un tableau H dont l'élément H[i] contient le nombre de moyennes contenues dans l'intervalle correspondant. Pour représenter graphiquement l'histogramme, on utilisera simplement la fonction matplotlib.pyplot.plot avec l'option "o" pour tracer des ronds.
On remarquera que la moyenne mi est toujours un multiple de 1/N, puisqu'elle est égale à un nombre entier inférieur à N divisé par N. En conséquence, il est judicieux de prendre M=N, afin que toutes les valeurs possibles de mi soient représentées dans l'histogramme et que chaque classe corresponde bien à une valeur possible.
Pour que l'histogramme représente une densité de probabilité, il faudra diviser toutes les valeurs du tableau H, d'une part par le nombre Nt de moyennes calculées, d'autre part par la largeur de chaque intervalle (ici 1/M).
[4] Écrire une fonction qui génére l'histogramme, pour N, Nt et M donnés.
[5] Que se passe-t-il si M>N ? Est-il intéressant de choisir M< N ?
[6] Tracer l'histogramme, pour des valeurs de N et Nt croissantes. Choisir une valeur de N qui donne un écart-type bien visible.
[7] Comparer la distribution obtenue à la distribution normale rappelée ci-dessous, où μ est l'espérance et σ l'écart-type de la moyenne :
avec μ=E(x)=p et σ2=var(x)/N=p(1-p)/N. D'après le théorème de la limite centrale, les valeurs moyennes suivent effectivement une loi normale pour N assez grand.
[8] À partir de quelle valeur de N peut-on considérer que la loi normale est vérifiée pour la distribution des valeurs moyennes ?
3. Distribution binomiale
Pour étudier la distribution binomiale par une simulation de Monte-Carlo, on effectue Nt fois le tirage de N valeurs de x. On génère un histogramme H, où H[n] est le nombre de tirages donnant n fois la valeur 1. Cet histogramme représente des probabilités; il faudra donc diviser H par le nombre total de tirages.
[9] Écrire la fonction générant l'histogramme de la distribution binomiale.
[10] La tester pour N=10, et pour des valeurs croissantes de Nt. Tester pour différentes valeurs de la probabilité p. Vérifier que la moyenne est bien Np.
[11] Écrire une fonction calculant la moyenne et l'écart-type de la distribution binomiale. Comparer à l'écart-type théorique .
4. Solution
import numpy
import random
from matplotlib.pyplot import *
class Binomiale:
def __init__(self,p):
self.p = p
def tirage(self):
x = random.random()
if x<self.p:
return 1
else:
return 0
def moyenne(self,N):
somme = 0
somme2 = 0
for k in range(N):
t = self.tirage()
somme += t
somme2 += t*t
moyenne = somme*1.0/N
variance = somme2*1.0/N-moyenne*moyenne
ecart = numpy.sqrt(variance*1.0/N)
return (moyenne,variance,ecart)
def distribution_moyennes(self,N,N_tirages,M):
hist = numpy.zeros(M)
for t in range(N_tirages):
(m,v,e) = self.moyenne(N)
i = int(m*M)
hist[i] += 1
return (numpy.arange(M)*1.0/M,hist/N_tirages*M)
def binomiale(self,N,N_tirages):
hist = numpy.zeros(N+1)
for t in range(N_tirages):
s = 0
for k in range(N):
s += self.tirage()
hist[int(s)] += 1
return hist/N_tirages
Voici tout d'abord un calcul de moyenne dans le cas p=0,5, pour N=1000 échantillons :
p=0.5
binom = Binomiale(p)
N=1000
(m,v,e) = binom.moyenne(N)
print((m,v,e)) --> (0.495, 0.249975, 0.015810597711661632)
Le dernier nombre est l'estimation de l'écart-type de la moyenne, qui indique quelle confiance on doit donner à l'estimation de la moyenne.
On génère l'histogramme des valeurs moyennes pour Nt=10000 calculs de la moyenne, puis on le compare à la distribution normale, avec μ=p et . Le nombre de calculs de la moyenne doit être grand par rapport à M. La moyenne est égale à un nombre entier inférieur à N (le nombre de tirages de 1) divisé par N. Il y a donc N valeurs possibles de la moyenne et on a intérêt à choisir M=N.
M=N
Nt=100000
(moy,hist) = binom.distribution_moyennes(N,Nt,M)
x=numpy.linspace(0,1,1000)
sigma=numpy.sqrt(p*(1-p)/N)
normale=1.0/(numpy.sqrt(2*numpy.pi)*sigma)*numpy.exp(-(x-p)**2/(2*sigma**2))
figure()
plot(moy,hist,"ko")
plot(x,normale,"r")
xlabel("m")
ylabel("P")
grid()
figA.pdf
Voici un histogramme de la distribution binomiale, pour N=10 :
p=0.5
binom = Binomiale(p)
hist = binom.binomiale(10,100000)
figure()
plot(hist,"o-")
xlabel("n")
ylabel("P")
grid()
title("Binomiale")
figB.pdf
et pour une probabilité différente :
p=0.3
binom = Binomiale(p)
hist = binom.binomiale(10,100000)
figure()
plot(hist,"o-")
xlabel("n")
ylabel("P")
grid()
title("Binomiale")
figC.pdf
