Une énigme pour moi…

Une énigme pour moi…

Bonjour les gens ! En effectuant quelques recherches pour mon prochain livre à paraître uniquement sur mathweb.fr, je suis tombé sur le problème suivant:

On choisit au hasard un nombre entre 0 et 1, et on répète cela jusqu’à ce que la somme des nombres choisis dépasse 1.

Quelle est le nombre moyen de nombres choisis ?

Le problème peut s’implémenter en Python à l’aide de la fonction suivante:

from random import random

def moyenne(n):
    T = 0
    for k in range(n):
        s, i = 0, 0

        while s < 1:
            s += random()
            i += 1
            
        T += i
    
    return T / n

print( moyenne( 10**9 ) )

On obtient alors:

>>> moyenne(10**7)
2.7186881
>>> moyenne(10**8)
2.71829368
>>> moyenne(10**9)
2.718286834

La moyenne semble être égale à “e” (ce n’est qu’une conjecture).

Le problème est que je n’arrivais pas à établir la moindre formule…

La solution de Philippe Rackette

Fort heureusement, il y a parmi les visiteurs et abonné·e·s de ce site de bien meilleurs probabilistes que moi ! C’est ainsi que l’un d’eux m’a proposé la solution suivante (le lien vers le PDF a été mis en commentaire, mais il va sans doute arriver un jour où le fichier correspondant sera effacé, donc j’ai préféré le retranscrire ici.

Posons \(S_{k+1} = P(X>k+1)\) pour tout entier naturel \(k\). Alors,$$S_{k+1}=\int_{x_1+\cdots+x_{k+1}<1,\ x_i \geqslant0} 1\text{d}x_1\cdots \text{d}x_{k+1}$$En effet, la densité de probabilité est homogène et la probabilité qu’il faille plus de k+1 nombres nécessite que la somme des k+1 premiers nombres soit strictement inférieure à 1.

On peur alors écrire pour k non nul:$$S_{k+1}=\int_{x_{k+1}=0}^{x_{k+1}=1^-}\left( \int_{x_1+\ \cdots\ +x_k<1-x_{k+1},\ x_i\geqslant0}1\text{d}x_1\cdots\text{d}x_k \right)\text{d}x_{k+1}$$

En effectuant les changements de variables:$$u_i = \frac{1}{1-x_{k+1}}x_i$$pour i variant de 1 à k, on obtient:$$S_{k+1}=\int_{x_{k+1}=0}^{x_{k+1}=1^-}\left( \int_{u_1+\ \cdots\ +u_k<1,\ u_i\geqslant0}(1-x_{k+1})^k\text{d}u_1\cdots\text{d}u_k \right)\text{d}x_{k+1}$$que l’on peut aussi écrire: $$S_{k+1}=\int_{x_{k+1}=0}^{x_{k+1}=1^-}\left( \int_{u_1+\ \cdots\ +u_k<1,\ u_i\geqslant0}1\text{d}u_1\cdots\text{d}u_k \right)(1-x_{k+1})^k\text{d}x_{k+1}$$c’est-à-dire:$$S_{k+1}=\int_{x_{k+1}=0}^{x_{k+1}=1^-}S_k(1-x_{k+1})^k\text{d}x_{k+1}$$

Comme \(S_k\) est indépendant de \(x_{k+1}\), on a alors:$$S_{k+1}=S_k\int_{x_{k+1}=0}^{x_{k+1}=1^-}(1-x_{k+1})^k\text{d}x_{k+1}$$

Or,$$\int_0^1(1-x)^k\text{d}x=\left[-\frac{(1-x)^{k+1}}{k+1}\right]_0^1=\frac{1}{k+1}$$ d’où:$$S_{k+1}=\frac{1}{k+1}S_k.$$

De plus, \(S_1=1\) car l’événement “plus d’un nombre est nécessaire pour obtenir une somme strictement plus grande que 1” est certain. On a alors:$$\begin{align}S_k&=\frac{1}{k}S_{k-1}\\& = \frac{1}{k}\times\frac{1}{k-1}S_{k-2}\\&=\cdots\\&=\frac{1}{k}\times\frac{1}{k-1}\times\cdots\times\frac{1}{k-i}S_{k-i-1}\\&=\cdots\\&=\frac{1}{k}\times\frac{1}{k-1}\times\cdots\times\frac{1}{k-(k-2)}S_{k-(k-2)-1}\\&=\frac{1}{k}\times\frac{1}{k-1}\times\cdots\times\frac{1}{2}S_{1}\\& = \frac{1}{k!}\end {align}$$

Ainsi,$$P(X\leqslant k)=1-P(X>k)=1-S_k=1-\frac{1}{k!}$$

Par suite, on a:$$\begin{align}P(X=k) & =P(X\leqslant k)-P(X\leqslant k-1)\\&=\left(1-\frac{1}{k!}\right)-\left(1-\frac{1}{(k-1)!}\right)\\&=\frac{1}{(k-1)!}-\frac{1}{k!}\\&=\frac{k-1}{k!} \end{align}$$

On peut alors calculer le nombre moyen que l’on cherchait en déterminant l’espérance de X:$$\begin{align} E(X) & =\sum_{k\geqslant1} k \times P(X=k)\\& = \sum_{k\geqslant1}k\times\frac{k-1}{k!}\\&=\sum_{k\geqslant1}\frac{k-1}{(k-1)!}\\&=\sum_{p\geqslant0}\frac{1}{p!}\\&=\text{e}.\end{align}$$

L’énigme est ainsi résolue! Me voilà rassuré… La conjecture de départ était donc bonne, mais c’est bien mieux de voir pourquoi!

En mathématiques, tout est une question de point de vue: je n’avais pas le bon dès le départ et je ne pouvais donc pas avoir la bonne idée !

Stéphane Pasquet
Stéphane Pasquet

1 réflexion au sujet de « Une énigme pour moi… »

Philippe

Philippe RackettePublié le  9:39 - Avr 13, 2021

Est-ce que la valeur 1/6 pour P(X=2) est “confirmée” par la simulation ? Cela me semble peu. Je la verrais plutôt à 1/2 selon le “raisonnement” suivant : après un premier tirage aléatoire, qu’on sait d’avance ne pas être le dernier, la somme cumulée est “probablement autour de 1/2” donc il y a une chance sur deux pour qu’au tirage suivant la somme cumulée dépasse 1 et une chance sur deux qu’elle ne dépasse pas 1. En continuant… Si après 2 tirages, 1 n’a toujours pas était atteint, le somme cumulée est “probablement autour de 2/3” et il y a 1 chance sur 3 de ne pas dépasser la somme cumulée 1 au 3e tirage. Si après k tirages, la somme cumulée 1 n’a toujours pas était atteinte, elle est “probablement autour de k/(k+1)” et il y a une chance sur k+1 de ne pas dépasser la somme cumulée 1 au k+1 ème tirage. En admettant (ou expliquant) cela on peut en déduire que l’espérance est bien 1+1/2+1/3!+1/4!+1/5!+… = e

    Stéphane Pasquet

    Stéphane PasquetPublié le  10:58 - Avr 14, 2021

    Le raisonnement exposé dans l’article est loin d’être correct je pense… Donc 1/6, en effet, paraît peu. C’est d’ailleurs cette valeur qui indique que mon raisonnement est sans doute faux. Pour ce qui est de ce que vous expliquez, il faudrait un spécialiste des probabilités là… car tout ça est obscure pour moi (si Mickaël Launay venait à passer par ici…:-))

Philippe

Philippe RackettePublié le  11:55 - Avr 14, 2021

On pourrait commencer par essayer de montrer que pour tout entier k supérieur à 1, \(S_{k+1}=\frac{1}{k+1}S_k\) avec \(S_k=\int_{x_1+…x_k \leq, 1 x_i \geq 0} 1 dx_1 … dx_k\) (la densité de probabilité est homogène). Au rang 1, cela signifie que la moitié inférieure gauche du carré unitaire a une aire de 1/2, au rang 2, que la pyramide (O,I,J,K) du repère canonique de l’espace a un volume de 1/6. \(S_k\) correspond à l’espérance de \(X_1+…X_k\) sachant que \(X_1+…X_k\) est inférieur à 1 où les \(X_k\) suivent une loi uniforme sur [0,1[… et \(S_k=\int_{x_1+…+x_k \leq 1, x_i \geq 0} 1 dx_1 … dx_k\) représente \(P(X > k)=1-P(X \leq k)\)

Christophe BalPublié le  11:07 - Avr 16, 2021

Le code suivant confirme l’intuition de P. Rackette.

from random import random

moy = 0
rep = 10**6

for _ in range(rep):
    N1 = random()
    N2 = random()

    if N1 + N2 > 1:
        moy += 1
        
moy /= rep

print(moy)
Stéphane Pasquet

Stéphane PasquetPublié le  11:14 - Avr 16, 2021

Ce dernier code donne :


0.499667
0.500387
0.499749

qui “montre” que la probabilité que la somme des deux premiers nombres pris au hasard semble être égale à \(\frac{1}{2}\) (précision apportée pour les visiteurs).

Christophe BalPublié le  11:22 - Avr 16, 2021

Une idée géométrique… À creuser je pense.

Pour \(X_1\) et \(X_2\), l’ensemble des points \((X_1 ; X_2)\) est dans un carré de côté 1, et les points tels que \(X_1 + X_2 \leqslant 1\) correspondent à la moitié basse du carré d’aire \(\frac{1}{2}\) d’où le 1/2 de proba.

Pour \(X_1,\ X_2, X_3\), l’ensemble des points \((X_1 ; X_2 ; X_3)\) est dans un cube de côté 1 et les points tels que \(X_1 + X_2 + X_3 \leqslant 1\) correspondent à une pyramide de volume \(\frac{1}{2}\times\frac{1}{3}\) d’où le 1/6 de proba.

Il faudrait continuer avec de n simplexes droits.

Stéphane Pasquet

Stéphane PasquetPublié le  11:34 - Avr 16, 2021

Cette façon de voir est en effet intéressante… mais pour la généraliser à la dimension \(n\), c’est une autre histoire… 🙂

Christophe BalPublié le  12:25 - Avr 16, 2021

Voici comment généraliser.

Prenons \(n\) vecteurs unitaires \(e_1,\ , e_2 ,\ \ldots ,\ e_n\) formant un repère orthonormé.

Il est connu que le volume du polytope formé par ces vecteurs est \( | det(e_1 , e_2 , … , e_n) |\) soit ici 1 .

Il y a n! façons d’ordonner les vecteurs. Pour un ordre choisi \(e_{k_1} ,\ e_{k_2} ,\ \ldots ,\ e_{k_n}\) , on a un repère affine \( (e_{k_1} ,\ e_{k_2} ,\ \ldots ,\ e_{k_n}) \) d’origine \( e_{k_1} \) avec des abus évidents de notations.

Dans ce repère, nous avons un volume représentant de façon naturelle les points \( (X_{k_1} ; X_{k_2} ; \ldots ; X_{k_n}) \) tels que \( X_{k_1} + X_{k_2} + \ldots + X_{k_n} \leqslant 1\).

On a donc qu’un tel volume aura \( \frac{| det(e_1 , e_2 , … , e_n) |}{n!} = \frac{1}{n!} \) pour volume.

CHRISTOPHE BALPublié le  4:05 - Avr 16, 2021

Ma réponse n’est pas claire dans le choix des repères (ce n’est pas bon en fait). J’ai du mal à exprimer ce choix. L’idée est d’aller chercher les simplexes dans chaque “coin” afin de couvrir le polytope.

Christophe BalPublié le  5:35 - Avr 16, 2021

Je vais essayer de dégager du temps pour rédiger un truc au propre.

Je vous tiendrais au courant.

    Stéphane Pasquet

    Stéphane PasquetPublié le  6:03 - Avr 16, 2021

    Si vous rédigez un document au format PDF (c’est sans doute plus simple que sur un site internet), n’hésitez pas à partager le lien: peut-être que d’autres personnes seraient intéressées 🙂

LAU BOULANGERPublié le  6:18 - Avr 18, 2021

Bonjour Stéphane,
en ce qui concerne ta première interrogation c’est effectivement un problème d’indépendance.Si on pose X ta variable aléatoire, il faut calculer l’espérance de cette variable aléatoire.
P(X=1)=0
P(X=2)=0.5, il faut calculer une intégrale double sur le domaine 0<x<1 et 1-x<y2) c’est à dire P(X1+X2+X3>1|X1+Y<1)
Donc calculer l'aire du domaine 0<x<1, 0<y<1-x, 1-x-y<z2).
On obtient 1/6.
On continue ainsi, j’imagine que l’on peut faire une récurrence, mais j’ai la flemme de rédiger .
Ensuite l’espérance n’est rien d’autre que la somme des inverse des factorielles, du coup on trouve le bon résultat.
Bravo pour ton travail sur LaTeX, en tout cas.
Par contre si tu veux intégrer des résultats statistiques dans LaTeX, tu pourrais regarder Sweave pour R en utilisant l’IDE Rstudio.
Bon dimanche
amicalement
Laurent

LAU BOULANGERPublié le  6:21 - Avr 18, 2021

Désolé je n’avais pas vu la réponse de Philippe !

CHRISTOPHE BALPublié le  9:38 - Avr 18, 2021

Merci Philippe.

Il resterait à passer via un argument géométrique au lieu de calculs d’intégrales pour que ce soit plus joli (je n’ai toujours pas le temps de taper ceci proprement mais c’est promis je le ferais d’ici mi juin et je posterais ceci ici).

Laissez votre message