Les scripts de la ruine

Les programmes du laboratoire Python du DM3, complétés, exécutés et confrontés aux valeurs exactes. Comme dans le DM, tout est écrit avec des boucles et des accumulateurs, sans listes. Attention aux scripts qui contiennent random() : ce sont de vrais tirages, les sorties affichées sont celles de nos exécutions, et les vôtres différeront. C’est précisément le sujet de la dernière section.

1. Une partie sous vos yeux

Le cœur de la simulation du DM, extrait et rendu bavard : une marche aléatoire sur $\llbracket 0\,;N \rrbracket$ qui part de $k$ et affiche la fortune coup par coup, jusqu’à l’absorption. Ici $k = 5$ et $N = 10$, comme sur la figure des préliminaires.

from random import random

k, N, p = 5, 10, 0.5

fortune = k
coups = 0
while 0 < fortune < N:
    if random() < p:
        fortune = fortune + 1
    else:
        fortune = fortune - 1
    coups = coups + 1
    print(coups, ":", fortune)

Sortie (une exécution parmi d’autres) :

1 : 6
2 : 5
3 : 4
4 : 5
5 : 6
6 : 7
7 : 8
8 : 7
9 : 8
10 : 7
11 : 6
12 : 7
13 : 8
14 : 7
15 : 8
16 : 7
17 : 8
18 : 9
19 : 10

Dix-neuf coups, des hésitations autour de $7$, et Alice rafle tout. Relancez : vous obtiendrez une autre partie, parfois expédiée en cinq coups, parfois interminable.

2. Le calcul exact (questions Q10 et Q11)

Le premier programme du laboratoire, complété : la formule close de Q8.b), valable pour $p \neq \frac{1}{2}$ (dans le cas équitable, le DM a démontré en Q4 que $\ \forall k \in \llbracket 0\,;N \rrbracket,\ r_k = 1 - \frac{k}{N}$). La fonction seuil cherche la plus petite fortune commune $a$ pour laquelle la ruine dépasse $0{,}99$.

def ruine(k, N, p):
    rho = (1 - p)/p
    return (rho**N - rho**k)/(rho**N - 1)

def seuil(p):
    a = 1
    while ruine(a, 2*a, p) <= 0.99:
        a = a + 1
    return a

print(ruine(10, 20, 18/37))
print(seuil(18/37))
print(ruine(85, 170, 18/37))

Sortie :

0.6319687776828178
85
0.9900058784715651

La première valeur est le $r_{10}$ de la question Q10.a) : environ $0{,}632$, au lieu de $0{,}5$ pour un jeu équitable. La deuxième répond à Q11.b) : à partir de $85$ jetons chacun, la roulette ruine Alice avec une probabilité qui dépasse $0{,}99$, et la troisième ligne le vérifie. Aucun hasard dans ce script : ces trois nombres sont exacts (aux arrondis de la machine près) et retomberont à l’identique chez vous.

3. La simulation (questions Q26.a et Q26.b)

Le second programme du laboratoire, complété : il ne connaît aucune formule, il joue parties parties et compte. Il renvoie la fréquence des ruines et le nombre moyen de coups observés.

from random import random

def simule(k, N, p, parties):
    ruines = 0
    total_coups = 0
    for i in range(parties):
        fortune = k
        coups = 0
        while 0 < fortune < N:
            if random() < p:
                fortune = fortune + 1
            else:
                fortune = fortune - 1
            coups = coups + 1
        if fortune == 0:
            ruines = ruines + 1
        total_coups = total_coups + coups
    return ruines/parties, total_coups/parties

print(simule(10, 20, 18/37, 10000))
print(simule(1, 100, 0.5, 10000))

Sortie (une exécution) :

(0.6328, 97.1094)
(0.9908, 92.9378)

Première ligne, la roulette de Q26.a) : la fréquence $0{,}6328$ colle à la valeur exacte $0{,}632$ du script précédent. Seconde ligne, le cas équitable de Q26.b) avec $k = 1$ et $N = 100$ : les valeurs exactes du DM sont $r_1 = 0{,}99$ et $D_1 = 1 \times 99 = 99$, et la simulation répond $0{,}9908$ et $92{,}9$. Le calcul et la simulation tombent d’accord : c’est la meilleure validation possible de l’un comme de l’autre.

4. L’expérience de la question Q26.c)

On reprend le même appel avec $N = 10\,000$, cinq fois de suite. Les valeurs exactes sont $r_1 = 1 - \frac{1}{10\,000} = 0{,}9999$ et $D_1 = 1 \times 9\,999 = 9\,999$ coups. Nécessite la fonction simule du script 3.

for essai in range(5):
    print(simule(1, 10000, 0.5, 10000))

Sortie (nos cinq exécutions, environ une minute de calcul au total) :

(1.0, 12869.5328)
(0.9999, 11137.1658)
(0.9999, 10142.9102)
(1.0, 3341.4784)
(0.9998, 16761.8576)

Le contraste annoncé par l’énoncé est là. La première valeur renvoyée est d’une fidélité totale : cinq estimations entre $0{,}9998$ et $1$, pour une valeur exacte de $0{,}9999$. La seconde est déboussolée : de $3\,341$ à $16\,762$ coups, du simple au quintuple, autour d’une valeur exacte de $9\,999$ qu’aucune des cinq exécutions n’approche vraiment.

5. Pourquoi la moyenne ne se stabilise pas

Les deux estimateurs n’ont pas le même ennemi. La fréquence des ruines moyenne des valeurs qui valent $0$ ou $1$ : rien ne peut y peser lourd, et sur $10\,000$ parties son écart-type est inférieur à $\frac{1}{2\sqrt{10\,000}} = 0{,}005$. La durée moyenne, elle, additionne des nombres sans plafond, et le paradoxe de Q21 travaille contre elle : presque toutes les parties sont courtes (partie de $1$ jeton, la ruine arrive le plus souvent en quelques coups), mais de très rares parties voient la fortune s’élever loin des barrières et durent alors des millions de coups. Ce sont ces parties presque introuvables qui portent l’essentiel de la moyenne : selon qu’une exécution en attrape zéro, une ou deux, son estimation s’effondre ou s’envole.

On peut le chiffrer. La méthode de conditionnement de R1, appliquée au moment d’ordre deux de la durée (on admet encore son existence), donne dans le cas équitable, pour un départ à $k = 1$, en notant $T_1$ la durée, aléatoire, d’une telle partie :

$$\forall N \geqslant 2, \qquad V(T_1) = \frac{N\,(N-1)\,(N-2)}{3},$$

formule que nous avons vérifiée par résolution exacte du système linéaire pour plusieurs petits $N$. Pour $N = 10\,000$, l’écart-type d’une seule partie vaut environ $577\,000$ coups, cinquante-huit fois la moyenne $D_1 = 9\,999$ ; et l’écart-type de la moyenne de $10\,000$ parties vaut encore environ $5\,800$ coups, soit $58\ \%$ de la valeur à estimer. Nos cinq exécutions, dispersées de $3\,341$ à $16\,762$, sont exactement dans cette fourchette : le programme est juste, le calcul de Q21.a) est juste, et c’est la moyenne empirique qui n’est pas un bon instrument ici. Vous tenez la version numérique du sommet du DM : un événement peut être certain, et son attente moyenne hors de portée de toute expérience.

Sources

• DM3 du chapitre A1, annexe « Le laboratoire Python » et question Q26 : les deux programmes ci-dessus en sont les complétions exactes.
• DM3, questions Q10, Q11, Q21 : les valeurs exactes utilisées comme références ($r_{10} \approx 0{,}632$, seuil $a = 85$, $r_1 = 0{,}99$, $D_1 = k\,(N-k)$).
• La formule de la variance de la durée : vérifiée par résolution exacte (fractions) du système de conditionnement pour $N = 4$, $6$, $8$, $12$, puis appliquée à $N = 10\,000$.