{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pile ou face"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Marc - Lorenzi - 16 avril 2018"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Succès, échec, quelques remarques philosophiques sur les générateurs de nombres aléatoires, attente d'un succès ... "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Succès, échec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Variables de Bernoulli"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $(\\Omega, P)$ un espace probabilisé. Soit $p\\in[0,1]$. Soit $A\\subset \\Omega$ un événement de probabilité $p$. Les éléments de $\\Omega$ sont les issues d'une expérience. Pour tout $\\omega\\in\\Omega$ nous dirons que l'expérience réussit (relativement à l'ensemble $A$) lorsque $\\omega\\in A$, sinon l'expérience échoue. Le réel $p$ est donc la probabilité de réussite de l'expérience. Par exemple, $p=\\frac 1 2$ pour un jeu de pile ou face avec une pièce parfaite, $p=\\frac 1 5$ lorsque $A$ est l'événement \"je suis gaucher\", $p\\simeq 5\\ 10^{-8}$ pour l'événement \"je gagne au loto\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parlons en termes de variables aléatoires, ce sera plus efficace. Soit $X:\\Omega\\to\\{0,1\\}$ une variable aléatoire ne pouvant prendre que les valeurs 0 et 1. Une telle v.a. est appelée variable aléatoire de Bernoulli. On pose $p=P(X=1)$. Le réel $p$ est appelé le paramètre de Bernoulli de la v.a. et on dit que $X$ suit une loi de Bernoulli de paramètre $p$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Simulation "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On peut simuler un variable de Bernoulli par la fonction ci-dessous. La fonction calcule un réel aléatoire appartenant à $[0,1]$ (fonction `random.random`). Si ce réel est dans $[0,p]$ elle renvoie 1, sinon elle renvoie 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rand(p):\n",
    "    if random.random() <= p: return 1\n",
    "    else: return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print([rand(0.3) for k in range(100)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ainsi, la fonction `rand` est une variable de Bernoulli de paramètre $p$. Euh mais non, ça ne va pas du tout !!! On voit bien que l'ensemble d'arrivée est $\\{0, 1\\}$ mais il n'y a pas d'ensemble de départ ! En fait `rand` n'est même pas une fonction au sens mathématique puisque deux appels successifs à cette fonction renvoient des valeurs qui peuvent être différentes. On peut contourner ce problème en imaginant qu'il y a un paramètre $\\omega$ implicite qui appartient à un univers $\\Omega$ qui est ... implicite aussi. Lors de $n$ évaluations successives de `rand(p)`, on obtient $n$ valeurs $X_1(\\omega),\\ldots,X_n(\\omega)$ où $X_1,\\ldots,X_n$ sont des v.a. qui suivent une loi de Bernoulli de paramètre $p$. \n",
    "\n",
    "Un générateur de nombres aléatoires bien écrit nous permettra d'affirmer avec confiance que ces v.a. sont indépendantes (enfin autant qu'on peut l'espérer)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Remarque__ : les générateurs de nombres aléatoires ne sont pas des objets magiques, ils sont pour la plupart bâtis autour de suites récurrentes du type $u_{n+1}=f(u_n)$, le paramètre $u_n$ appartenant à un ensemble plus ou moins compliqué et caché à l'utilisateur. Pour information, à l'heure où est écrit ce notebook, Python utilise comme générateur le \"Mersenne Twister\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 Moyennes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tout le monde le sait, le nombre moyen de fois que `rand(p)` renvoie 1 tend vers $p$ lorsque le nombre de tirages tend vers l'infini. Ceci soulève quelques questions :\n",
    "- Est-ce vrai ?\n",
    "- Si oui, pourquoi ?\n",
    "- Si non, pourquoi pas ? Et c'est non non ou un peu non ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici une fonction effectuant $n$ tirages aléatoires avec une probabilité de réussite $p$. La fonction renvoie la liste des tirages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def tirages(n, p):\n",
    "    return [rand(p) for k in range(n)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(tirages(50, 0.3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Que contient réellement cette liste ? Une fois encore les variables aléatoires sont nos amies. Sans justifier complètement ce que je vais dire, je peux considérer que le $k$ième élément de cette liste est $X_k(\\omega)$ où $\\omega\\in\\Omega$ est un élément inconnu d'un univers inconnu, et $X_k$ est une variable aléatoire suivant une loi de Bernoulli de paramètre $p$. Mieux, les v.a. $X_1,\\ldots,X_n$ sont indépendantes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `moyennes` prend en paramètres une liste $[x_1,\\ldots,x_n]$. Elle renvoie la liste $[y_1,\\ldots,y_n]$ où $y_k=\\frac 1 k\\sum_{i=1}^k x_i$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def moyennes(s):\n",
    "    somme = 0.\n",
    "    s1 = []\n",
    "    for k in range(len(s)):\n",
    "        somme += s[k]\n",
    "        s1.append(somme / (k + 1))\n",
    "    return s1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On regroupe en une seule fonction les appels aux deux fonctions précédentes : la fonction `pile_face` renvoie les moyennes des tirages. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pile_face(n, p):\n",
    "    s = tirages(n, p)\n",
    "    return moyennes(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Avec les notations ci-dessus, `pile_face(n, p)` renvoie la liste $[Y_1(\\omega),\\ldots,Y_n(\\omega)]$ où $\\omega\\in\\Omega$ et $Y_k=\\frac 1 k \\sum_{i=1}^k X_i$ pour $k=1,\\ldots,n$. Il n'y a plus qu'à tester. Plusieurs fois, parce que à chaque fois c'est différent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n, p = 1000, 0.6\n",
    "s = pile_face(n, p)\n",
    "plt.plot(s, 'k-')\n",
    "plt.plot(n* [p], 'r-')\n",
    "#plt.axis([1, n, p - 0.3, p + 0.3])\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Oui, bon, ben oui. Les moyennes ont bien l'air de converger vers $p$. Mais qu'entend-on par là ? La réponse à laquelle on rêve est que cela signifie que pour tout $\\omega\\in\\Omega$, $Y_n(\\omega)\\to p$ lorsque $n\\to\\infty$. Mais cela NE PEUT PAS être vrai. Imaginons que j'aie beaucoup, beaucoup de chance, je réussis à tous les coups. On aura alors $Y_n(\\omega)=1$ pour tout $n$, la limite sera 1 et pas $p$. En fait, en continuant à imaginer, on voit que n'importe quel réel entre 0 et 1 est une limite possible pour la suite $(Y_n)$. Pire, il se pourrait très bien que pour certains $\\omega$, $Y_n(\\omega)$ diverge ... Toute la question est : qu'entend-on par beaucoup, beaucoup de chance ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prenons-nous à rêver que $Y_n(\\omega)$ tend vers $p$ pour __presque tous__ les $\\omega\\in\\Omega$ ? Mathématiquement, cela revient à considérer l'événement $\\{Y_n\\to p\\}$. A-t-on $P(Y_n\\to p) = 1$ ? Affirmatif, mais c'est difficile à montrer. C'est la __loi forte des grands nombres__. Que signifie-t-elle exactement ? Soit $E=\\{\\omega\\in\\Omega, Y_n(\\omega)\\not\\to p\\}$. La loi forte des grands nombres nous dit que $P(E)=0$. Cela ne signifie pas que $E$ est vide, mais que si on tombe sur un élément de $E$ on n'a vraiment, vraiment pas de chance (ou beaucoup, beaucoup :-).\n",
    "\n",
    "Nous allons prouver le résultat plus faible suivant :"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Loi faible des grands nombres__ : Pour tout $\\varepsilon >0$, $P(|Y_n-p|> \\varepsilon)\\to 0$ lorsque $n$ tend vers l'infini."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4 Un peu de théorie"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $(X_n)_{n\\ge 1}$ une suite de v.a. indépendantes suivant une même loi de Bernoulli de paramètre $p\\in[0,1]$. Pour tout $n\\ge 1$, posons $S_n=\\sum_{k=1}^n X_k$. Soit $Y_n=\\frac 1 n S_n$. Soit $\\varepsilon> 0$. On utilise l'inégalité de Bienaymé-Tchebychev, démontrée en cours :\n",
    "$$P(|S_n-E(S_n)|> n\\varepsilon)\\le \\frac{V(S_n)}{n^2\\varepsilon^2}$$\n",
    "Maintenant, par linéarité de l'espérance, $E(S_n)=\\sum_{k=1}^n E(X_k)=np$. Et comme les $X_k$ sont indépendantes, $V(S_n)=\\sum_{k=1}^n V(X_k)=np(1-p)$. On réinjecte pour obtenir\n",
    "$$P(|S_n-np|> n\\varepsilon)\\le \\frac{np(1-p)}{n^2\\varepsilon^2}$$\n",
    "\n",
    "Bien évidemment, $|S_n-np|> n\\varepsilon$ si et seulement si $|nY_n-np|> n\\varepsilon$, ou encore $|Y_n-p|>\\varepsilon$. Ainsi,\n",
    "$$P(|Y_n-p|> \\varepsilon)\\le \\frac C n$$\n",
    "où $C=\\frac{p(1-p)}{\\varepsilon^2}$. D'où le résultat cherché par le théorème des gendarmes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Premier succès dans un tirage à pile ou face"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici deux affirmations maintes fois entendues dans les bars, autour d'un apéro pastis-cacahuètes :\n",
    "\n",
    "1- Si je joue vingt fois au loto j'ai plus de chances de gagner que si je joue 10 fois ...\n",
    "\n",
    "2- Ça fait dix fois de suite que je perds au loto ! Donc les 10 prochaines fois j'aurai plus de chances de gagner ... \n",
    "\n",
    "Ne remettons pas en doute la première affirmation (elle est vraie, mais cela ne m'intéresse pas ici). Intéressons nous plutôt à le seconde, elle est la raison d'être du joueur addict :-)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $(X_1,X_2,\\ldots)$ une suite de variables aléatoires indépendantes suivant une loi de Bernoulli de paramètre $p\\ne 0$, définies sur un espace probabilisé $(\\Omega, P)$. Soit $N:\\Omega\\to\\mathbb N$ définie par $N(\\omega)=\\min\\{n\\in\\mathbb N^*, X_n(\\omega)=1\\}$. On voit tout de suite un petit problème : et si pour tout entier $n$, $X_n(\\omega)=0$ ? On peut dans ce cas convenir que $N(\\omega)=\\infty$ où $\\infty$ est un symbole sans signification particulière.\n",
    "\n",
    "__Question__ : quelle est la loi de $N$ ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Commençons par une petite simulation. La fonction `premier_succes` est en fait la fonction $N$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def premier_succes(p):\n",
    "    count = 0\n",
    "    while rand(p) != 1: count += 1\n",
    "    return count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = [premier_succes(0.05) for k in range(10000)]\n",
    "plt.plot(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bon, c'est le fouillis. Mais c'est la loi de $N$ qui nous intéresse, pas les valeurs prises individuellement par $N$. Déterminons donc les probabilités $P(N=k)$ pour $k\\in\\mathbb N$. La fonction `histogramme` fait le travail."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def histogramme(s):\n",
    "    m = max(s)\n",
    "    h = (m + 1) * [0]\n",
    "    c = 0\n",
    "    for k in s:\n",
    "        h[k] = h[k] + 1\n",
    "        c = c + 1\n",
    "    for i in range(len(h)): h[i] = h[i] / c\n",
    "    return h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = [premier_succes(0.05) for k in range(10000)]\n",
    "plt.plot(histogramme(s), 'k-')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "C'est plus joli. La fonction $k\\mapsto P(N=k)$ est assez décroissante et tend super vite vers 0 lorsque $k$ tend vers l'infini. Apparemment. Vrai ou faux ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Un peu de théorie"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Avant de démarrer, voici un petit résultat sur les séries."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Proposition__ : soit $q\\in]-1,1[$. La série $\\sum q^n$ est une série convergente, et $\\sum_{n=0}^\\infty q^n=\\frac 1 {1-q}$.\n",
    "\n",
    "__Démonstration__ : Soit $N\\in\\mathbb N$. On a $\\sum_{n=0}^N q^n=\\frac{1-q^{N+1}}{1-q}$ et cette quantité tend bien vers $\\frac 1 {1-q}$ lorsque $N$ tend vers l'infini. Une telle série est appelée série géométrique de raison $q$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Après cette petite digression, revenons à notre fonction de premier succès. Nous reprenons les notations de l'introduction. Soit $\\omega\\in\\Omega$. On a $N(\\omega)=k$ si et seulement si pour tout $i<k$, $X_i(\\omega)=0$, et $X_k(\\omega)=1$. Ainsi, $\\{N=k\\}=\\left(\\cap_{i=1}^{k-1}\\{X_i=0\\}\\right)\\cap\\{X_k=1\\}$. Mais les $X_i$ sont indépendantes, la probabilité de l'intersection est donc le produit des probabilités :\n",
    "$$P(N=k)=P(X_k=1)\\prod_{i=1}^{k-1}P(X_i=0)=p(1-p)^{k-1}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On a fait un petit oubli ... et $P(N=\\infty)$ ? Eh bien $\\{N=\\infty\\}$ est le complémentaire de la réunion de tous les événements $\\{N=k\\}$, événements deux à deux incompatibles. Donc, $P(N=\\infty)=1-\\sum_{k=0}^\\infty P(N=k)$. Poursuivons :\n",
    "$$P(N=\\infty)=1-\\sum_{k=0}^\\infty p(1-p)^k=1-p\\sum_{k=0}^\\infty (1-p)^k$$\n",
    "On reconnaît là une série géométrique de raison $0\\le 1-p< 1$. Donc\n",
    "$$P(N=\\infty)=1-p\\frac{1}{1-(1-p)}=\\ldots 0$$\n",
    "L'événement $\\{N=0\\}$ n'est pas impossible, mais il est de probabilité nulle."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tout ceci suggère la définition suivante :\n",
    "\n",
    "__Définition__ : On dit qu'une variable aléatoire $X$ à valeurs dans $\\mathbb N$ suit une loi géométrique de paramètre $p\\in]0,1]$ lorsque, pour tout $k\\in\\mathbb N$, $P(X=k)=p(1-p)^{k-1}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici une illustration de la loi géométrique ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def geom(p, k):\n",
    "    return p * (1 - p) ** (k - 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = [geom(0.05, k) for k in range(1, 10000)]\n",
    "plt.grid()\n",
    "plt.plot(s[:200], 'r-')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Superposons la simulation et la théorie."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n, p = 10000, 0.05\n",
    "s = [premier_succes(p) for k in range(n)]\n",
    "s1 = [geom(p, k) for k in range(1, n)]\n",
    "plt.plot(histogramme(s), 'k-')\n",
    "plt.plot(s1[:max(s) + 1], 'r-')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ça colle. Et même plutôt bien. Et le loto dans tout ça ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 La propriété d'oubli de la loi géométrique"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $X$ une v.a. suivant une loi géométrique de paramètre $p$. Soit $n\\in\\mathbb N$. On a $P(X>n)=\\sum_{k=n+1}^\\infty p(1-p)^{k-1}=p(1-p)^n\\sum_{k=n+1}^\\infty(1-p)^{k-1-n}$. Le changement d'indice $k'=k-n-1$ dans la série nous montre que sa somme est $\\frac 1 {1-(1-p)}=\\frac 1 p$. Ainsi, $P(X>n)=(1-p)^n$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prenons maintenant deux entiers $m$ et $n$. On a $P(X>m+n| X>n)=\\frac{P(X>m+n,X>n)}{P(X>n)}=\\frac{P(X>m+n)}{P(X>n)}=\\frac{(1-p)^{m+n}}{(1-p)^n}=(1-p)^m=P(X>m)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On a utilisé dans notre calcul le fait que $\\{X>m+n\\}\\subset \\{X>n\\}$. Donc $\\{X>m+n\\}\\cap \\{X>n\\}=\\{X>m+n\\}$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Revenons à notre loto gagnant. Quelle est la probabilité de devoir jouer au moins $m+n$ fois avant de gagner, sachant qu'on a perdu $n$ fois ? Eh bien c'est la probabilité de devoir jouer au moins $m$ fois avant de gagner, sachant rien du tout. La morale de l'histoire c'est que ce n'est pas parce qu'on sait qu'on a déjà joué avant qu'on a plus de chances de gagner après. Et donc celui qui n'a jamais joué a autant de chances de gagner que celui qui a passé sa vie à perdre :-)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}