En statistique inférentielle, on cherche à connaître l'inconnu. Pour ce faire, on émet des hypothèses ou on estime des grandeurs partiellement inconnues,
La simulation permet de prendre le problème à l'envers, en créant un modèle qui se comporte comme l'échantillon qu'on observe, et en estimant les probabilités difficiles à calculer par des fréquences. Il arrive que, de ces simulations, on puisse inférer (on est là pour ça) des conjectures, qui serviront à élaborer les algorithmes de test ou d'estimation. Certes, dire "on prend 2 parce que la simulation a suggéré que, pour 2, la probabilité est de 0,95", ce n'est pas très mathématique, mais la théorie qui disait elle aussi "on prend 1,96; concrètement 2" est très souvent basée sur des approximations normales de lois compliquées. On va donc expérimenter et observer sans expliquer ce qu'on observe, mais en se servant de ces observations pour expliquer pourquoi les statisticiens font comme ci et pas comme ça.
=Estimations=
Pour faire des statistiques, il faut un échantillon de données aléatoires ou non. Et pour avoir des données sous "Python", le plus simple est de les fabriquer sous "Python". Par exemple si on veut faire des statistiques sur les 100 premiers carrés d'entiers, on peut fabriquer une liste contenant ces 100 nombres:
donnees=[n**2 for n in range(100)]
print(len(donnees))
Moyenne.
Pour calculer la moyenne des nombres qui sont dans "donnees", on les additionne et on divise la somme par le nombre de nombres qu'il y a dans "donnees":
def moyenne(tableau):
 return sum(tableau, 0.0) / len(tableau)
print(moyenne(donnees))
L'algorithme est améliorable puisque si une donnée n'est pas numérique, il ne donne qu'un message d'erreur.
Variance.
La variance est définie comme la moyenne des carrés des écarts à la moyenne:
def variance(tableau):
 m=moyenne(tableau)
 return moyenne([(x-m)**2 for x in tableau])
print(variance(donnees))
Une variante pour la variance est donnée par la formule de Huyghens: Moyenne des carrés moins le carré de la moyenne.
Écart-type.
L'écart-type est défini comme la racine carrée de la variance:
def ecartype(tableau):
 return variance(tableau)**0.5
print(ecartype(donnees))
Échantillons.
On peut créer un échantillon de 100 nombres gaussiens d'espérance 16 et d'écart-type 2, puis calculer sa moyenne et son écart-type:
from random import *
echantillon=[gauss(16,2) for n in range(100)]
print(moyenne(echantillon))
print(ecartype(echantillon))
On voit que la moyenne est proche de 16 et l'écart-type proche de 2. C'est rassurant. Mais si on y regarde de plus près, on voit un problème: En prenant des échantillons plus petits, on s'attend à ce que leurs moyenne et écart-type fluctuent mais que la moyenne des moyennes (sur beaucoup de petits échantillons) soit 16 et que la moyenne des écarts-types soit proche de 2. C'est vrai pour la moyenne des moyennes mais visiblement pas pour la moyenne des écarts-types:
from random import *
m=[] #liste des moyennes des echantillons
s=[] #liste des ecarts-types des echantillons
for n in range(10000):
 echantillon=[gauss(16,2) for k in range(5)]
 m.append(moyenne(echantillon))
 s.append(ecartype(echantillon))
print(moyenne(m)) # Voisin de 16, c'est rassurant!
print(moyenne(s)) # Largement plus petit que 2!
print(moyenne(s)*2**0.25)
print(ecartype(m)) # Le moyennage resserre les ecarts-types
print(2/5**0.5) # en les divisant par la racine de la taille de l'echantillon
En théorie, le nombre par lequel on doit multiplier la moyenne des écarts-types pour estimer l'écart-type de la population est formula_1. Ce n'est pas le cas ici: Il semble que l'algorithme de "Python" pour effectuer des tirages avec remise introduise un biais. Par contre on découvre expérimentalement ici que la moyenne des écarts-types doit être multipliée par formula_2 pour avoir un estimateur sans biais de l'écart-type...
=Intervalles de confiance=
Pour des fréquences.
Comme "Python" sait faire des tirages sans remise, on peut constituer une liste ordonnée de "pour" et de "contre", et y puiser des échantillons au hasard. On peut estimer la proportion d'entre eux qui donne de faux espoirs au maire (au moins 50 "pour" parmi les 100).
from random import *
population=['contre' for n in range(57000)]+['pour' for n in range(43000)]
shuffle(population)
print(len(population))
p=len([n for n in range(1000) if len([v for v in sample(population,100) if v=='pour'])>=50])
print(p)
print(p/1000)
L'échantillon choisi par l'institut de sondage, s'il a réellement été choisi au hasard, était favorable au maire avec une probabilité d'environ 0,1. Cette probabilité n'est pas si ridicule que ça, et l'institut de sondage aurait pu répondre au maire "c'est la faute à pas de chance": Il est tombé sur les 10 % d'échantillons favorables... Un épisode analogue s'est déroulé lors des élections présidentielles de 1995, le ministre du budget de l'époque ayant un peu trop vite cru aux sondages !
Plus sérieux (et plus prudent, pour éviter la vindicte de l'ancien maire, désormais dans l'opposition, et qui a maintenant le temps de mener une croisade contre les instituts de sondage) eût été la publication par l'institut de sondage, d'un intervalle de confiance, par exemple à 95% (c'est-à-dire un intervalle qui contient en moyenne 95% des échantillons). Expérimentalement, on peut s'inventer un intervalle et compter la fréquence des échantillons de 100 personnes qui sont dedans. Ce sera un estimateur de la probabilité que l'échantillon soit représentatif de l'ensemble de la population:
from random import *
h=0.1 
p=0
for n in range(1000):
 pourcentage=len([v for v in sample(population,100) if v=='pour'])/100
 if pourcentage>0.43-h and pourcentage<0.43+h:
 p+=1
print(p/1000)
On voit que l'intervalle "[0,33 ; 0,53]" obtenu avec "h=0,1" est un intervalle à 95 %. En modifiant la valeur de "h" on constate que si "h" diminue (l'intervalle rétrécit), on perd de la confiance (la probabilité qu'il soit bon diminue aussi). On trouve par tâtonnements la valeur de "h" pour laquelle la confiance de l'intervalle vaut 95 %, puis par changement de la taille de l'échantillon, on peut conjecturer le lien entre "h" et la taille de l'échantillon.
Pour des moyennes.
Là encore, on peut facilement tester des intervalles de confiance, pour des moyennes de variables aléatoires normales, par exemple d'espérance 16, d'écart-type 2 et indépendantes entre elles:
from random import *
h=0.1 
p=0
for n in range(1000):
 echantillon=[]
 for k in range(100):
 echantillon.append(gauss(16,2))
 m=moyenne(echantillon)
 if m>16-h and m<16+h:
 p+=1
print(p/1000)
On découvre que l'intervalle de confiance "[15,9 ; 16,1]" donné ci-dessus (pour "h=0,1") est à environ 40% de confiance. En modifiant la valeur de "h", on retrouve expérimentalement que pour celle-ci égale à environ formula_3, l'intervalle est à 95 % de confiance.
=Test d'équirépartition=
En lançant un dé 100 fois, on constate que le 6 est sorti un peu souvent par rapport aux autres nombres:
On se demande si le dé est équilibré. Pour cela, on se choisit comme critère de test la somme des carrés des écarts aux fréquences théoriques: formula_4:
d2=(0.15-1/6)**2+(0.16-1/6)**2+(0.17-1/6)**2+(0.16-1/6)**2+(0.16-1/6)**2+(0.2-1/6)**2
Soit mais encore, que faire avec "d2": Est-ce qu'on doit dire que 0,0015 est anormalement élevé et que le dé est truqué, ou que 0,0015 est suffisamment petit pour attribuer ces résultats à la fluctuation d'échantillonnage? Pour le savoir, on va simuler 10000 lancers de dés et calculer l'équivalent de "d2" pour chacun d'entre eux, puis faire une étude statistique sur le "d2" observé. La réponse statistique à la question "Qu'est-ce qui est normal?" est en général fournie par les déciles: On dira que le dé est vraisemblablement truqué si le "d2" observé est supérieur au neuvième décile de la série. Pour calculer ce décile, on va devoir trier les données.
from random import *
khi2=[]
for n in range(10000):
 effectifs=[0]*6
 for k in range(100):
 effectifs[randint(0,5)]+=1
 dsquare=0
 for e in effectifs:
 dsquare+=(e/100-1/6)**2
 khi2.append(dsquare)
khi2.sort()
decile9=khi2[int(0.9*len(khi2))]
print(decile9)
print(d2)
"d2" est environ 10 fois plus petit que le neuvième décile de la série, donc on se trompe de plus de 10 % en considérant que le dé est truqué: Il est parfaitement normal pour autant qu'on sache.
=Problème des rencontres d'Euler=
Euler affirme que, si on joue 19 fois au jeu de treize, le joueur A gagnera "probablement 12 fois". On peut détailler cette affirmation, en jouant 19 fois au jeu pour chaque échantillon, et en répétant l'expérience 1000 fois (donc 1000 échantillons de 19 parties). Voici le script:
jeu1=[str(v)+' '+c for v in valeurs for c in couleurs]
jeu2=[str(v)+' '+c for v in valeurs for c in couleurs]
from random import *
pg=[0]*20
for n in range(1000):
 rencontres=0
 for p in range(19):
 shuffle(jeu2)
 for i in range(32):
 if jeu1[i]==jeu2[i]:
 rencontres+=1
 break
 pg[rencontres]+=1
print(pg)
Le tableau affiché par ce script est celui des parties gagnées par A parmi les 19; il s'agit d'un tableau d'effectifs (le total des nombres entiers affichés est d'ailleurs 1000). Voici un exemple de diagramme en bâtons de ce tableau, traîtreusement déguisé en histogramme:
On voit que le mode de cette série est 12, c'est peut-être ce que voulait dire Euler. Mais peut-être aussi voulait-il dire que l'espérance de la série (en fait, sa moyenne puisqu'on fait de la statistique) est proche de 12. C'est bien le cas, cette moyenne étant proche de 11,9 (avec un écart-type proche de 2). En fait le phénomène aléatoire simulé ci-dessus est assez bien modélisé par une variable aléatoire binomiale de paramètres 19 et formula_5, dont l'espérance est 12 et l'écart-type formula_6.