Le plus grand de trois dés

Le programme.

import numpy as np
import numpy.random as rd

n = 100000
d1 = rd.randint(1, 7, n)
d2 = rd.randint(1, 7, n)
d3 = rd.randint(1, 7, n)
print(np.mean(np.maximum(np.maximum(d1, d2), d3)))

La borne haute de rd.randint est exclue. Et np.maximum compare deux tableaux terme à terme, d'où les deux appels emboîtés pour trois dés.

La valeur exacte. Pour $k$ entre $1$ et $6$, l'événement $[M\le k]$ signifie que les trois dés sont inférieurs ou égaux à $k$, donc

$$P(M\le k)={\left(\frac{k}{6}\right)}^3=\frac{k^3}{216},P(M=k)=\frac{k^3-(k-1{)}^3}{216}.$$

Les six probabilités sont donc $\frac{1}{216},\frac{7}{216},\frac{19}{216},\frac{37}{216},\frac{61}{216},\frac{91}{216}$, et

$$E(M)=\frac{1+14+57+148+305+546}{216}=\frac{1071}{216}\approx 4,958.$$

Contrôle. Les six numérateurs somment à $216$ : c'est bien une loi de probabilité.

La précision. Sur $100000$ parties, la simulation rend typiquement $4,95$ ou $4,96$ : les deux premières décimales sont fiables, la troisième ne l'est pas.

Le point à retenir. Annoncer « $E(M)=4,9583$ » d'après la simulation serait abusif. C'est le calcul, et lui seul, qui donne $\frac{1071}{216}$ ; porter un regard critique sur une simulation, c'est refuser d'en écrire plus de chiffres qu'elle n'en garantit.