Pour faire des simulations en proba, on a besoin de nombres pseudo-aléatoires.

Ensembles en Python

Simulation avec Python

Obtention

Recette

Versions jusqu'à la 2.7

Python 2.7 calculait les nombres pseudo-aléatoires avec l'algorithme de Wichmann-Hill, basé sur 3 générateurs congruentiels linéaires en parallèle (extrait du code source - on remarquera l'allusion à Fermat dans l'avant-dernière ligne):

        # This part is thread-unsafe:
        # BEGIN CRITICAL SECTION
        x, y, z = self._seed
        x = (171 * x) % 30269
        y = (172 * y) % 30307
        z = (170 * z) % 30323
        self._seed = x, y, z
        # END CRITICAL SECTION

        # Note:  on a platform using IEEE-754 double arithmetic, this can
        # never return 0.0 (asserted by Tim; proof too long for a comment).
        return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0

Ceci suggère quelques exercices d'arithmétique:

Vérifier que les modules 30269, 30307 et 30323 sont premiers.
Vérifier que les multiplicateurs 171, 172 et 170 sont primitifs modulo leurs modules respectifs.
Simuler l'un des générateurs congruentiels (x par exemple est une suite géométrique de raison 171 dans le corps de Galois $\mathbb {Z} /30269\mathbb {Z}$ ).
Il résulte du point 2 ci-dessus que la suite x est périodique de période 30268. L'addition de la dernière ligne donnerait une période de $30268\times 30306\times 30322$ . Calculer ce nombre et le décomposer en facteurs premiers...

Version 3.2

Python 3.2 utilise le Mersenne twister, implémenté en C. La période est annoncée comme égale à $2^{19\,937}-1$ . Démontrer la primalité de ce nombre de Mersenne est difficile, mais un exercice intéressant est de démontrer la primalité de son exposant 19 937.

Une fois qu'on sait construire un nombre entier pseudo-aléatoire compris entre 0 et un entier fixe N-1, il suffit de le diviser par N pour simuler un réel pseudo-aléatoire compris entre 0 (inclus en théorie) et 1 (exclu). Ce réel pseudo-aléatoire est celui appelé random et sa loi est uniforme sur l'intervalle allant de 0 à 1. Une transformation affine permet alors d'obtenir à partir de lui un réel pseudo-aléatoire uniforme entre deux réels donnés, puis une troncature permet de retrouver des entiers aléatoires, modulo un nombre entier donné.

Le module random de Python 3.2 ne sait pas simuler de variables aléatoires binomiales ou de Poisson, mais il fournit des variables aléatoires normales avec l'algorithme de Kinderman et Monahan :

NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0)



    def normalvariate(self, mu, sigma):
        random = self.random
        while 1:
            u1 = random()
            u2 = 1.0 - random()
            z = NV_MAGICCONST*(u1-0.5)/u2
            zz = z*z/4.0
            if zz <= -_log(u2):
                break
        return mu + z*sigma

Il est basé sur deux variables aléatoires u1 et u2 uniformes entre 0 et 1 (sauf que u2 ne peut être nulle). u1 est ensuite centrée, et divisée par u2. Le tout est normalisé avec la constante magique $4{\frac {e^{-{\frac {1}{2}}}}{\sqrt {2}}}$ , élevé au carré et on finit par lui appliquer une transformation affine pour que ses espérance et écart-type soient ceux désirés.

En plus, une variable aléatoire normale centrée réduite peut être obtenue avec gauss qui utilise l'algorithme de Box-Muller.

Python 3.2 possède aussi des simulateurs de tirage sans remise, et même de permutations aléatoires, où on mélange les cartes en multipliant les échanges de deux cartes (technique analogue à celle consistant à couper le jeu plusieurs fois de suite):

        for i in reversed(range(1, len(x))):
            # pick an element in x[:i+1] with which to exchange x[i]
            j = int(random() * (i+1))
            x[i], x[j] = x[j], x[i]

Cet algorithme illustre le fait que toute permutation est engendrée par des transpositions.

Syntaxe

Toutes ces merveilles se trouvant dans le module random, on doit charger celui-ci pour faire des simulations avec Python:

from random import *

Variables uniformes

Continue

Le plus simple à obtenir c'est le fameux nombre pseudo-aléatoire entre 0 et 1:

from random import *
print(random())

Pour avoir un nombre entre -2 et 3, on peut faire

from random import *
print(uniform(-2,3))

Entière

Pour lancer un dé, on peut utiliser l'une des méthodes suivantes:

from random import *
print(randint(1,6))
print(randrange(1,7))
print(choice(range(1,6)))
from math import *
print(ceil(6*random()))

Variables continues

La loi exponentielle de paramètre quelconque (ici 3) peut être simulée:

from random import *
print(expovariate(3))

Une variable gaussienne peut être simulée par deux méthodes (ici centrée et réduite):

from random import *
print(normalvariate(0,1))
print(gauss(0,1))

random simule aussi d'autres lois continues comme les lois Beta, Gamma, de Weibull etc.

Variables binomiales et de Poisson

On l'a vu plus haut, random ne simule pas de variables binomiales. Mais une variable binomiale de paramètres n et p pouvant être définie comme somme de n variables de Bernoulli de probabilité p indépendantes entre elles, est simulable par boucle:

from random import *

def Bernoulli(p):
    if random()<p:
        return 1
    else:
        return 0


def binomial(n,p):
    somme=0
    for k in range(n):
        somme+=Bernoulli(p)
    return somme


print(binomial(25,0.2))

Cet algorithme nécessite d'être amélioré, en effet rien n'empêche de l'appeler avec des valeurs de p supérieures à 1 ou des valeurs de n non entières.

Pour simuler une loi de Poisson, on peut utiliser une variable binomiale qui l'approche, comme par exemple

def Poisson(l):
    return binomial(1000000,l/1000000)

mais c'est très long à calculer. Une méthode plus rapide utilise la loi exponentielle de même paramètre:

from random import *
def Poisson(l):
    t=1.0/expovariate(l)
    return int(t)

Permutations et arrangements

Permutations

En 1741, dans les Mémoires de l'Académie des Sciences de Berlin, Leonhard Euler publiait un article titré Calcul de la probabilité du jeu de rencontre. Le nom initial du jeu de rencontre était jeu de treize parce qu'il se jouait à 13 cartes. Mais Euler généralise le jeu à n cartes.

Il le définit ainsi:

Le jeu de rencontre est un jeu de hasard, où deux personnes ayant chacune un entier jeu de cartes, en tirent à la fois une carte après l'autre, jusqu'à ce qu'il arrive, qu'elles rencontrent la même carte: et alors l'une des deux personnes gagne. Or, lorsqu'une telle rencontre n'arrive point du tout, alors c'est l'autre des deux personnes qui gagne. Cela posé, on demande la probabilité, que l'une et l'autre de ces deux personnes aura de gagner.

Dans cet excellent article (16 pages en Français),

Euler montre que

Pourvu donc que le nombre de cartes ne soit pas moindre que 12, l'espérance de A sera toujours à celle de B à peu près comme 12 à 7 ... Ou bien parmi 19 jeux qu'on joue, il y en aura probablement 12 qui font gagner A, et 7 qui feront gagner B.

On constate qu'Euler utilisait le mot espérance là où aujourd'hui on écrit probabilité, et sa conclusion s'écrit algébriquement $P\left(A\right)\simeq {\frac {12}{19}}$ et $P\left(B\right)\simeq {\frac {7}{19}}$ , expliquant pourquoi dans les diligences de l'époque, les gens pariaient à 12 contre 7 sur une rencontre dans ce jeu.

Pour vérifier cela avec 2 jeux de 32 cartes, ça peut prendre du temps (bien que le jeu s'arrête dès la première rencontre). Alors la simulation offerte par Python permet de rapidement simuler un grand nombre de parties (Python est une sorte de diligence supersonique).

Pour mélanger un jeu de 32 cartes, on doit d'abord le construire, avec la technique vue précédemment, sauf que le mélange ne peut être fait que sur une liste et pas sur un ensemble:

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 *

shuffle(jeu2)


diagnostic='pas de rencontre'
for i in range(32):
    if jeu1[i]==jeu2[i]:
        diagnostic='rencontre sur la '+str(i)+'ème carte, qui est le '+jeu1[i]
        break

print(diagnostic)

On a créé deux jeux de 32 cartes, mélangé l'un des deux (jeu2) puis comparé carte à carte les deux jeux. Par défaut le texte de sortie est pas de rencontre. En effet c'est seulement s'il y a une rencontre qu'il est modifié, et remplacé par le nom et le numéro de la carte commune aux deux jeux. L'instruction break sert à éviter de continuer à jouer après la fin du jeu.

Tirage sans remise

Un magicien demande à une personne de l'assistance de tirer une carte d'un jeu de 32. Quelle est la probabilité que ce soit l'as de pique?

Pour simuler l'expérience, on va utiliser la fonction choice qui permet de tirer une carte. On va alors construire le jeu de cartes comme dans l'article précédent, à ceci près que choice attend une liste et non un ensemble:

valeurs={1,7,8,9,10,'Valet','Dame','Roi'}
couleurs={'carreau','cœur','pique','trèfle'}
univers=[str(v)+' '+c for v in valeurs for c in couleurs]

from random import *
for n in range(100):
    if choice(univers)=="1 pique":
        print('victoire !')

Une estimation de la probabilité de l'évènement peut alors se faire en comptant le nombre de fois que le mot victoire a été écrit: Il est, en pourcents, la fréquence de l'évènement.

Pour jouer au poker, on peut simuler le choix d'une main par un tirage de 5 éléments (sans répétition) parmi les 32:

valeurs={1,7,8,9,10,'Valet','Dame','Roi'}
couleurs={'carreau','cœur','pique','trèfle'}
univers=[str(v)+' '+c for v in valeurs for c in couleurs]

from random import *

hand=sample(univers,5)
print(hand)

Mathématiques avec Python et Ruby/Nombres pseudoaléatoires en Python

Sections