{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TP 10 : P-valeurs\n",
    "Le but du TP est de réaliser des simulations de variable aléatoire et des calculs de P-valeurs.\n",
    "\n",
    "| Numéro        | Compétence           | \n",
    "|:------------- |:--------------------:|\n",
    "|TS001|\tSavoir définir une hypothèse nulle|\n",
    "|TS002|\tConnaître les différences entre les deux types d’erreurs|\n",
    "|TS003|\tSavoir interpréter une p-valeur|\n",
    "|PY601|\tSavoir calculer une p-valeur|\n",
    "|PY602|\tÊtre capable de calculer une statistique avec scipy|\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.io as sio\n",
    "import scipy.stats as st"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ex. 0 : Échauffement\n",
    "1. Qu'est-ce que le module scipy.stats ?\n",
    "2. Que font les méthodes `pdf`, `cdf`, `sf`, `ppf` et `isf` de `scipy.stats.norm` ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "source": [
    "https://docs.scipy.org/doc/scipy/reference/stats.html     \n",
    "\n",
    "scipy.stats.norm correspond à une loi normale\n",
    "\n",
    " - pdf  $$ f(x) = \\frac{1}{\\sqrt{2 \\pi \\sigma^2}} \\exp^{-\\frac{(x-m)^2}{2 \\sigma^2}}$$\n",
    " - cdf  $$ F(x) = P(X \\leq x) = \\int_{-\\infty}^x f(t) dt $$\n",
    " - sf  $$sf(x) = 1 - F(x)$$\n",
    " - ppf\n",
    " - isf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "$$p = \\mathbb{P}(X < t)$$\n",
    "$$p = cdf(t)$$\n",
    "$$t = ppf(p)$$\n",
    "$$s = \\mathbb{P}(X \\geq t)$$\n",
    "$$ s = sf(t) $$\n",
    "$$ t = isf(s) $$\n",
    "\n",
    "et on a $$sf(t) = 1 - cdf(t)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3\\. Dessinez les densités de probabilité d'une loi normale centrée réduite et une autre loi Normale $\\mathcal{N}(\\frac{1}{2}, 1.5^2)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import norm\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "x = np.linspace(-4, 4, 100)\n",
    "plt.plot(x, norm.?(?))\n",
    "plt.plot(x, norm.?(?, loc=?, scale=?))\n",
    "plt.legend(['N(0,1)','N(1/2,1.5^2)'])\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "from scipy.stats import norm\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "x = np.linspace(-4, 4, 100)\n",
    "plt.plot(x, norm.pdf(x))\n",
    "plt.plot(x, norm.pdf(x, loc=0.5, scale=1.5))\n",
    "plt.legend(['N(0,1)','N(1/2,1.5^2)'])\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4\\. Calculez, pour $X$ une variable aléatoire suivant une loi normale $\\mathcal{N}(0,1)$ et $t \\in R^+$, $\\mathbb{P}(X < t)$ et $\\mathbb{P}(|X| < t)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = 1.96\n",
    "#P(X <t)\n",
    "pr1 = norm.?(?)\n",
    "pr2 = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "t = 1.96\n",
    "#P(X <t)\n",
    "pr1 = norm.cdf(t) #norm.cdf calcule toute l'aire \"à gauche\" de t. Faire un dessin de la gaussienne pour illustrer\n",
    "#P(|X| < t), 3 manieres de faire\n",
    "pr2 = norm.cdf(t) - norm.cdf(-t)\n",
    "pr3 = 1 - 2*(1 - norm.cdf(t))\n",
    "pr4 = 1 - 2*norm.sf(t)\n",
    "print(pr1, pr2, pr3, pr4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ex. 1: Pile ou face"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Générez $n=1000$ répétitions d’un échantillon de $p=10$ tirages pile ou face avec une pièce « honnête » (telle que: $pr=\\mathbb{P}(X=pile) =.5)$\n",
    "\n",
    ">**aide** : vous pouvez utiliser la méthode `choice` de `numpy.random`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n, p = 1000, 10\n",
    "rand = np.random.choice(?, size=[?, ?])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "n, p = 1000, 10\n",
    "rand = np.random.choice([0,1], size=[n, p]) #face est codé par 0, pile par 1\n",
    "print(rand.shape)\n",
    "print(rand[0:50,0:10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Comptez le nombre de pile de chaque tirage et affichez l'histogramme correspondant. Sachant qu'une loi binomiale tend vers une loi normale $\\mathcal{N}(p*pr,p*pr*(1-pr))$, tracez la densité de probabilité correspondante."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nb_pile =  ?\n",
    "plt.hist(nb_pile,bins=11)\n",
    "\n",
    "pr = 1/2\n",
    "\n",
    "x = np.linspace(-1, 11, 100)\n",
    "plt.plot(x, ?*norm.?(x, loc=?, scale=?))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "nb_pile = np.sum(rand, axis=1)\n",
    "print(nb_pile.shape)\n",
    "print(nb_pile[:10])\n",
    " \n",
    "plt.hist(nb_pile,bins=11)\n",
    "\n",
    "pr = 1/2\n",
    "x = np.linspace(-1, 11, 100)\n",
    "\n",
    "plt.plot(x, n*norm.pdf(x, loc=p*pr, scale=np.sqrt(p*pr*(1-pr))))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "source": [
    "Car la loi binômiale $B(n,p)$ tend vers une loi normale d'espérence $np$ et de variance $npq$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Comptez le nombre d’échantillons dans lesquels vous avez tiré de $k=0 $ à $p= 10$ fois piles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mods,freqs = ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "mods, freqs = np.unique(nb_pile, return_counts=True)\n",
    "mods, freqs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Comptez le nombre théorique moyen d’échantillon dans lesquels vous auriez du avoir de $k=0$ à $p=10$ fois piles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "import scipy.stats as st\n",
    "print(p,n)\n",
    "#loi binomiale de param pr pour p répetitions\n",
    "# pmf : probability mass function\n",
    "pr_observer_i_pile = np.array([st.binom.pmf(i, p, pr) for i in range(11)])\n",
    "pr_observer_i_pile*n # transformation des probas en effectif"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. Approchez ce nombre en utilisant la loi normale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "# param loi normale (TCL)\n",
    "normal = np.array([st.norm.pdf(i, p*pr, np.sqrt(p*pr*(1-pr))) for i in range(11)])\n",
    "normal*n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "6. Affichez vos résultats. À savoir les effectifs observés, les effectifs prédits par la loi binomiale ainsi que les effectifs prédits par la loi normale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "plt.plot(np.arange(len(freqs)), freqs ,'o-', label='Echantillons')\n",
    "plt.plot(np.arange(len(pr_observer_i_pile)), pr_observer_i_pile*n ,'x-r', label='Binomial')\n",
    "plt.plot(np.arange(len(normal)), normal*n ,'+:g', label='Normale')\n",
    "plt.legend()\n",
    "#Les échantillons sont légeement différents"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "7. que se passe t’il lorsque $n=1000000$ et $p= 20$ ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n, p = 1000000, 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "n, p = 1000000, 20\n",
    "rand = np.random.choice([0,1], size=[n, p])\n",
    "sums = np.sum(rand, axis=1)\n",
    "mods, freqs = np.unique(sums, return_counts=True)\n",
    "\n",
    "# Il y a une chance de tirer 0 ou 20 piles !\n",
    "if 0 not in mods:\n",
    "    mods = np.concatenate([[0], mods])\n",
    "    freqs = np.concatenate([[0], freqs])\n",
    "if 20 not in mods:\n",
    "    mods = np.concatenate([mods, [20]])\n",
    "    freqs = np.concatenate([freqs, [20]])\n",
    "\n",
    "\n",
    "binom = n * np.array([st.binom.pmf(i, p, 0.5) for i in range(p+1)])\n",
    "normal = n * np.array([st.norm.pdf(i, p/2, np.sqrt(p/4)) for i in range(p+1)])\n",
    "plt.plot(np.arange(p+1), freqs ,'o-', label='Echantillons')\n",
    "plt.plot(np.arange(p+1), binom ,'x-r', label='Binomial')\n",
    "plt.plot(np.arange(p+1), normal ,'+:g', label='Normale')\n",
    "plt.legend()\n",
    "# on est bien plus proche de l'équivalence loi binomiale/normale/echantillons"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ex. 2: La binômiale"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On propose un questionnaire comprenant 10 questions avec chacune deux réponses possibles,l’une vraie et l’autre fausse. Pour tester si une personne interrogée répond correctement, on adopte la règle suivante : si 7 réponses ou plus sont correctes on admet que la personne n'a pas répondu au hasard.\n",
    "\n",
    "1. Donner les deux hypothèses du tests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "source": [
    "l'élève à travaillé : p > 1/2\n",
    "l'élève a répondu au hasard : p = 1/2\n",
    "\n",
    "$\\begin{cases}\n",
    "        \\mathcal{H}_0: p  = 1/2\\\\\n",
    "        \\mathcal{H}_1: p > 1/2\n",
    "       \\end{cases}\n",
    "       $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "source": [
    "2. Quelle est la p-valeur du test d’une personne qui a donné 8 réponses correctes ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "source": [
    "$$\\mbox{p-val} = P(X \\geq 8) = 1 - P(X \\leq 8)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "n = 10 # 10 questions\n",
    "p = 1/2 # Hypothese H_0\n",
    "p_val = (1 - st.binom.cdf(8,n,p)) #evenement plus rare que au moins 8 réponses correctes\n",
    "p_val2 = st.binom.sf(8,n,p) # deux manieres de calculer\n",
    "print(p_val,p_val2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Quelle est la probabilité de décider qu’une personne à répondu correctement alors qu’elle a répondu au hasard ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si $X$ est une v.a. binômiale $B(n,p)$ avec $p=1/2$, \n",
    "$$\n",
    "P(X \\geq 7)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = (1 - st.binom.cdf(7,n,p)) # idem qu'avant\n",
    "print(alpha)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4\\. Que devient cette probabilité lorsque chacune des questions posées comporte trois réponses dont une seule est vraie ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = 1/3\n",
    "alpha = (1 - st.binom.cdf(6,n,p))\n",
    "alpha"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ex. 3: La grande échelle FDR (Bonus mais intéressant !!)\n",
    "\n",
    "Les industries pharmaceutiques ont toujours eu pour objectif de développer des méthodes permettant la production de nouvelles molécules thérapeutiques. L’approche privilégiée consiste à générer un grand nombre de molécules de les tester afin de décider si elles sont actives ou non. On parle alors de criblage ou screening. Nous allons supposer dans cet exercice que nous savons que parmi 1000 molécules testées, seulement 100 ont un effet avéré (et donc sont actives). Nous allons supposer que le risque de première espèce du test utilisé est de 5 % et que la puissance de ce test est de 0,7 et utiliser la théorie des tests statistiques vue en cours.\n",
    "\n",
    "1. Quelle est, en moyenne, le nombre de faux positifs (le nombre de molécules réellement inactives que le test va déclarer actives) ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "source": [
    "on pose les hypothèses :\n",
    "* H_0 : molécule inactive\n",
    "* H_1 : molecule active\n",
    "\n",
    "$\\alpha$ :P(decider H1 | H0 est vraie) <- les faux positifs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "# parmi toutes les réellement inactives, je commets 5 % d 'erreur \n",
    "alpha = 0.05 # mon risque alpha\n",
    "nb_H0 = 900 # les inactives\n",
    "alpha*nb_H0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Quelle est, en moyenne, le nombre de faux négatifs (le nombre de molécules actives que le test va déclarer inactives) ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "#beta : Decider H0 alors que H1 est vrai\n",
    "puissance = .7\n",
    "beta= 1-puissance \n",
    "nb_H1 = 1000-nb_H0\n",
    "beta*nb_H1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Donnez une estimation de la probabilité de déclarer une molécule active alors qu’elle ne l’est pas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "# L'ensemble des molecules décidées active : P(D1|H0) + P(D1|H1)\n",
    "decid_H1 = alpha*nb_H0 + (1-beta)*nb_H1\n",
    "\n",
    "print(f\"Nb molecules décidées active : {decid_H1}\")\n",
    "\n",
    "# parmi toutes les molécules décidées actives, quelle proportion de molécules ne le sont pas ?\n",
    "p_err = (alpha*nb_H0)/decid_H1 # part des réellement H0 parmi les décidées H1\n",
    "print(f\"proba : {p_err}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Qu’en concluez vous ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "source": [
    "Si je décide une molécule active, j'ai encore beaucoup d'incertitudes... Pourtant ma p valeur est standard à .05 !"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. Donnez la matrice de confusion « moyenne » du test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "|/|molécule active|molecule inactive|totaux|\n",
    "|---|---|---|---|\n",
    "|test positif|?|?|?|\n",
    "|test negatif|?|?|?|\n",
    "|totaux|100|900|1000|\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "source": [
    "\n",
    "|/|molcule active|molecule inactive|totaux|\n",
    "|---|---|---|---|\n",
    "|test positif|70|45|115|\n",
    "|test negatif|30|855|885|\n",
    "|totaux|100|900|1000|\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "mlenv",
   "language": "python",
   "name": "mlenv"
  },
  "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}