Dans ce TP on se propose de simuler une loi exponentielle avec Python.

La loi exponentielle est utilisée dans le cadre de durée de vie. On fabrique une ampoule et on se demande combien de temps elle va tenir. La durée de vie de l'ampoule est aléatoire, c'est donc une variable aléatoire et on peut s'interroger sur sa densité de probabilité.

Pour bien comprendre le fonctionnement de la loi exponentielle, on va passer par un jeu.

Dans ce jeu, à chaque tour, on jette un dé (ou une pièce, ou un osselet, ou ce qu'on veut…)

Par exemple on pourrait dire qu'à chaque tour on jette le dé et que si on fait 6, on a perdu. On dirait alors qu'on la probabilité $p = \frac{1}{6}$ de perdre à chaque tour.

De façon plus générale, on supposera qu'on a la probabilité $p$ de perdre à chaque tour – c'est à dire $p$ pas forcément égal à $\frac{1}{6}$.

On s'intéresse au nombre de tours où on a survécu.

On réalise d'abord une fonction qui simule le dé (ou tirage de cartes, ou autre chose…)

On veut écrire une fonction perdu(p) qui renvoie vrai avec la probabilité p.

from random import random

def perdu(p):
    """
    renvoie True avec la probabilité p
    """
    if p > random():
        return True
    else:
        return False

Entrez ce code et vérifiez que perdu(0.5) renvoie True à peu près une fois sur deux.

Une partie est une succession de lancers de dés. Chaque fois qu'on gagne, on compte +1. Quand c'est perdu, c'est fini. On s'intéresse au nombre de tours gagnés, donc à la durée de vie de la partie. On l'appelle T.

Ajoutez à la suite :

def T(p):
    """
    p: probabilité de perdre à chaque tour
    simule une partie et renvoie le nombre de lancer gagnants
    """
    n = 0
    while ...
        n = ...
    return n

Vous pouvez tester.

Si on simule, cela donne quelque chose comme :

>>> simulation(0.2)
10

Dans cet exemple, on a survécu 10 tours. Cela ne dit pas grand chose. Il faut répéter l'expérience et compter les résultats pour avoir une idée de la répartition des scores et de leur probabilité.

import matplotlib.pyplot as plt  # pour le graphique
import numpy as np               # contient une fonction de comptage

p = 0.2   # par exemple...
N = 10000 # nombre d'essais

# simulation de N parties     
essais = [simulation(p) for i in range(N)]
     
# comptage des scores obtenus et des effectifs correspondants
valeurs, effectifs = np.unique(essais , return_counts=True)
# on préfère des fréquences :
freqs = [e/N for e in effectifs]
     
# figure
plt.figure()
plt.bar(valeurs, freqs)
plt.title('Répartition des scores')
plt.xticks(valeurs)
plt.xlabel('durée de vie')
plt.ylabel('fréquence')
plt.show()

On reconnaît bien une forme exponentielle : $\lambda \cdot \text{e}^{-\lambda\cdot t}$.

Pour trouver la valeur de $\lambda$ correspondant à notre simulation, on peut exploiter la formule $E(T) = \frac{1}{\lambda}$.

Donc en prenant la moyenne de notre simulation et en inversant, on aura $\lambda$.

# à mettre avant le tracer du graphique
moyenne = np.mean(essais)
L = 1/moyenne
print(1/moyenne)

# on peut ajouter un graphique
# juste avant le plt.bar :
plt.plot(valeurs, [L*np.exp(-L*v) for v in valeurs], 'r')

Vous constatez que la courbe suit à peu près les barres.

On peut montrer que la bonne valeur pour $\lambda$ est : L = -np.log(1 - p). Vous pouvez essayer.