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,

On se fixe des probabilités a priori de se tromper dans ses estimations ou son test;
On cherche (difficile) pour quels choix des paramètres (intervalle ou estimateur) ces probabilités sont atteintes;
On prend ses décisions à partir des choix faits ci-dessus.

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 ${\sqrt {\frac {5}{4}}}$ . 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 ${\sqrt {\frac {10}{7}}}$ pour avoir un estimateur sans biais de l'écart-type...

Intervalles de confiance

Pour des fréquences

Une situation fictive

Sur les 100 000 électeurs d'une ville, 43 000 s'apprêtent à voter pour le maire sortant, mais celui-ci ne le sait pas. Alors il commande un sondage, et l'institut de sondage constitue un échantillon de 100 habitants sur lesquels 52 disent vouloir voter pour le maire. Fou de joie, celui-ci achète le champagne pendant la campagne, avant les élections, et se fait battre lors des élections. Il accuse l'institut de sondage d'avoir triché. Qu'en penser?

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))

#On choisit 1000 echantillons de 100 et on compte combien sont favorables au maire:
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 ${\frac {2\sigma }{\sqrt {100}}}$ , 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:

Résultat du tirage	1	2	3	4	5	6
Effectifs	15	16	17	16	16	20

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: $d2=\sum _{k=1}^{6}\left(f_{k}-{\frac {1}{6}}\right)^{2}$ :

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:

valeurs={1,7,8,9,10,'Valet','Dame','Roi'}
couleurs={'carreau','cœur','pique','trèfle'}
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 ${\frac {12}{19}}$ , dont l'espérance est 12 et l'écart-type ${\sqrt {19\times {\frac {12}{19}}\times {7}{19}}}\simeq 2,103$ .