Pour faire des simulations en proba, on a besoin de nombres pseudo-aléatoires. 
=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:
Version 3.2.
"Python" 3.2 utilise le Mersenne twister, implémenté en C. La période est annoncée comme égale à formula_3. 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" formula_4, é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. 
Dans cet excellent article (16 pages en Français), 
On constate qu'Euler utilisait le mot "espérance" là où aujourd'hui on écrit "probabilité", et sa conclusion s'écrit algébriquement formula_5 et formula_6, 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:
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:
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:
univers=[str(v)+' '+c for v in valeurs for c in couleurs]
from random import *
hand=sample(univers,5)
print(hand)

Pôles d'éco-citoyenneté/Les principaux Pôles d'éco-citoyenneté dans le Nord-Pas-de-Calais/Les CPIE/Chaîne des terrils