Simuler une loi discrète donnée par un tableau

Le principe. On empile les probabilités : les réels de $[0;0,2[$ donneront $0$, ceux de $[0,2;0,6[$ donneront $1$, ceux de $[0,6;0,9[$ donneront $2$, ceux de $[0,9;1[$ donneront $3$. Comme rd.random() rend un réel uniforme de $[0,1[$, chaque valeur sort avec la probabilité voulue, qui est exactement la longueur de son intervalle.

Les seuils. Ce sont les probabilités cumulées, que np.cumsum calcule : $0,2;0,6;0,9;1,0$.

Le programme.

import numpy as np
import numpy.random as rd

val = np.array([0, 1, 2, 3])
pro = np.array([0.2, 0.4, 0.3, 0.1])
cum = np.cumsum(pro)

def tirage():
    u = rd.random()
    k = 0
    while u >= cum[k]:
        k = k + 1
    return val[k]

rd.seed(3)
ech = np.array([tirage() for i in range(100000)])
print(np.mean(ech))
# 1.30469

La valeur exacte. $$E(X)=0\times 0,2+1\times 0,4+2\times 0,3+3\times 0,1=0,4+0,6+0,3=1,3.$$

Ce que la simulation prouve, et ce qu'elle ne prouve pas. Elle affiche $1,305$ contre $1,300$ : elle corrobore le calcul, elle ne le démontre pas. Et le troisième chiffre après la virgule n'est pas garanti — l'annoncer serait la faute que la compétence de regard critique vise justement.

Le point délicat. La boucle while s'arrête forcément, puisque le dernier seuil vaut $1$ et que rd.random() rend un réel strictement inférieur à $1$. Si les probabilités ne sommaient pas à $1$, l'indice k sortirait du tableau : c'est le contrôle à faire avant d'exécuter.