.. _tdnote2021rst: ============================ 1A - Enoncé 24 novembre 2020 ============================ .. only:: html **Links:** :download:`notebook `, :downloadlink:`html `, :download:`python `, :downloadlink:`slides `, :githublink:`GitHub|_doc/notebooks/exams/td_note_2021.ipynb|*` Correction de l’examen du 24 novembre 2020. .. code:: ipython3 from jyquickhelper import add_notebook_menu add_notebook_menu() .. contents:: :local: Exercice 1 : guérison --------------------- On commence par générer des données artificielles à partir de véritables données. .. code:: ipython3 import pandas df = pandas.read_csv("https://www.data.gouv.fr/en/datasets/r/63352e38-d353-4b54-bfd1-f1b3ee1cabd7", sep=";") gr = df[["jour", "rad", "dc"]].groupby(["jour"]).sum() gr.head() .. raw:: html
rad dc
jour
2020-03-18 1627 435
2020-03-19 2322 642
2020-03-20 3128 890
2020-03-21 3580 1041
2020-03-22 4188 1251
.. code:: ipython3 diff = gr.diff().reset_index(drop=False) diff.head() .. raw:: html
jour rad dc
0 2020-03-18 NaN NaN
1 2020-03-19 695.0 207.0
2 2020-03-20 806.0 248.0
3 2020-03-21 452.0 151.0
4 2020-03-22 608.0 210.0
On convertit la date en jour de l’année puis on simule un loi exponentielle de paramètre 14 pour avoir la date de sortie. .. code:: ipython3 import numpy.random as rnd def donnees_artificielles(hosp, mu=14, nu=21): dt = pandas.to_datetime(hosp['jour']) res = [] for i in range(hosp.shape[0]): date = dt[i].dayofyear h = hosp.iloc[i, 1] delay = rnd.exponential(mu, int(h)) for j in range(delay.shape[0]): res.append([date - int(delay[j]), date, 1]) h = hosp.iloc[i, 2] delay = rnd.exponential(nu, int(h)) for j in range(delay.shape[0]): res.append([date - int(delay[j]), date , 0]) return pandas.DataFrame(res, columns=["entree", "sortie", "issue"]) data = donnees_artificielles(diff[1:].reset_index(drop=True)) data.head() .. raw:: html
entree sortie issue
0 39 79 1
1 71 79 1
2 58 79 1
3 79 79 1
4 53 79 1
.. code:: ipython3 data.to_csv("examen2021.csv", index=False) Q1 ~~ On récupère les données. .. code:: ipython3 import pandas df = pandas.read_csv("http://www.xavierdupre.fr/enseignement/complements/examen2021.zip") df.head() .. raw:: html
entree sortie issue
0 49 79 1
1 27 79 1
2 73 79 1
3 74 79 1
4 48 79 1
Q2 : durée de guérison ~~~~~~~~~~~~~~~~~~~~~~ .. code:: ipython3 duree = df['sortie'] - df['entree'] duree = duree.values # conversion en numpy issue = df['issue'].values duree[:5], issue[:5] .. parsed-literal:: (array([30, 52, 6, 5, 31], dtype=int64), array([1, 1, 1, 1, 1], dtype=int64)) Q3 : estimateur Kaplan-Meier (1) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. code:: ipython3 t = 10 nt = duree[(duree >= t)].shape[0] dt = duree[(duree == t) & (issue == 0)].shape[0] st = 1. - dt / nt st .. parsed-literal:: 0.98965342710248 Q4 : courbe de Kaplan-Meier ~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. code:: ipython3 T = [0] St = [1.] for t in range(0, 150): nt = duree[(duree >= t)].shape[0] dt = duree[(duree == t) & (issue == 0)].shape[0] st = 1. - dt / nt T.append(t) St.append(st * St[-1]) Q5 : graphe ~~~~~~~~~~~ .. code:: ipython3 import matplotlib.pyplot as plt fig, ax = plt.subplots(1, 1, figsize=(4, 4)) ax.plot(T, St); .. image:: td_note_2021_17_0.png .. code:: ipython3 import numpy from lifelines import KaplanMeierFitter fig, ax = plt.subplots(1, 1, figsize=(10, 4)) ax.plot(T, St, label="custom", lw=10) kmf = KaplanMeierFitter() kmf.fit(duree, (issue == 0).astype(numpy.int32)) kmf.plot(ax=ax) ax.legend(); .. image:: td_note_2021_18_0.png Q6 : application aux données publiques ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Les données accessibles librement sur le portail `data.gouv.fr `__ recensent les entrées et les sorties des personnes sans relier une entrée et une sortie spécifique. Si N est personnes sont sorties guéries, on ne sait pas quand elles sont entrées. Donc le calcul ci-dessus n’est pas possible. Exercice 2 ---------- Q1 : t + 1 ~~~~~~~~~~ .. code:: ipython3 import numpy N = 10 M = numpy.zeros((N, N)) M[4, 5] = 1 M .. parsed-literal:: array([[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]]) .. code:: ipython3 def propagation(M): M2 = M.copy() M2[1:, :] = numpy.maximum(M2[1:, :], M[:-1, :]) M2[:-1, :] = numpy.maximum(M2[:-1, :], M[1:, :]) M2[:, 1:] = numpy.maximum(M2[:, 1:], M[:, :-1]) M2[:, :-1] = numpy.maximum(M2[:, :-1], M[:, 1:]) return M2 propagation(M.copy()) .. parsed-literal:: array([[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.], [0., 0., 0., 0., 1., 1., 1., 0., 0., 0.], [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]]) Q2 : après T itération ~~~~~~~~~~~~~~~~~~~~~~ .. code:: ipython3 def propagation_n(M, t): for i in range(t): M = propagation(M) return M propagation_n(M, 3) .. parsed-literal:: array([[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.], [0., 0., 0., 0., 1., 1., 1., 0., 0., 0.], [0., 0., 0., 1., 1., 1., 1., 1., 0., 0.], [0., 0., 1., 1., 1., 1., 1., 1., 1., 0.], [0., 0., 0., 1., 1., 1., 1., 1., 0., 0.], [0., 0., 0., 0., 1., 1., 1., 0., 0., 0.], [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]]) Q3 : vaccin ~~~~~~~~~~~ .. code:: ipython3 p = 0.3 vaccine = (numpy.random.rand(N, N) <= p).astype(numpy.int32) vaccine .. parsed-literal:: array([[0, 0, 0, 0, 0, 0, 1, 0, 1, 0], [0, 0, 0, 0, 0, 0, 1, 0, 0, 1], [1, 0, 0, 0, 0, 0, 1, 0, 1, 1], [0, 1, 0, 1, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 1, 1, 0, 0, 1, 0, 1, 0], [0, 0, 1, 0, 0, 0, 0, 1, 0, 0], [0, 0, 1, 0, 0, 0, 1, 0, 0, 1], [0, 0, 1, 0, 0, 0, 0, 1, 0, 0], [0, 0, 0, 1, 0, 0, 0, 1, 0, 0]]) .. code:: ipython3 def propagation_vaccine(M, vaccine): M2 = M.copy() M2[1:, :] = numpy.maximum(M2[1:, :], M[:-1, :]) M2[:-1, :] = numpy.maximum(M2[:-1, :], M[1:, :]) M2[:, 1:] = numpy.maximum(M2[:, 1:], M[:, :-1]) M2[:, :-1] = numpy.maximum(M2[:, :-1], M[:, 1:]) M2 = numpy.minimum(M2, 1 - vaccine) return M2 vaccine[4, 5] = 0 propagation_vaccine(M, vaccine) .. parsed-literal:: array([[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.], [0., 0., 0., 0., 1., 1., 1., 0., 0., 0.], [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]]) Q4 : après T heures ~~~~~~~~~~~~~~~~~~~ .. code:: ipython3 def propagation_n_vaccine(M, t, vaccine): for i in range(t): M = propagation_vaccine(M, vaccine) return M propagation_n_vaccine(M, 3, vaccine) .. parsed-literal:: array([[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.], [0., 0., 0., 0., 1., 1., 0., 0., 0., 0.], [0., 0., 0., 0., 1., 1., 1., 1., 0., 0.], [0., 0., 1., 1., 1., 1., 1., 1., 1., 0.], [0., 0., 0., 0., 1., 1., 0., 1., 0., 0.], [0., 0., 0., 0., 1., 1., 1., 0., 0., 0.], [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]]) Q5 : variation ~~~~~~~~~~~~~~ .. code:: ipython3 import pandas res = [] for p in [0.7, 0.8, 0.85, 0.86, 0.87, 0.88, 0.89, 0.9, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.0]: cont = [] for test in range(0, 10): vaccine = (numpy.random.rand(N, N) <= p).astype(numpy.int32) M[4, 5] = 1 vaccine[4, 5] = 0 M = propagation_n_vaccine(M, 3, vaccine) contamine = M.ravel().sum() cont.append(contamine) cont = numpy.array(cont) res.append(dict(p=p, prop=cont.mean())) df = pandas.DataFrame(res) df .. raw:: html
p prop
0 0.70 9.1
1 0.80 9.6
2 0.85 3.8
3 0.86 2.4
4 0.87 2.4
5 0.88 1.4
6 0.89 1.6
7 0.90 1.2
8 0.94 1.2
9 0.95 1.5
10 0.96 1.0
11 0.97 1.2
12 0.98 1.2
13 0.99 1.1
14 1.00 1.0
.. code:: ipython3 ax = df.plot(x="p", y="prop") ax.plot([0.7, 1.], [2, 2], "--"); .. image:: td_note_2021_33_0.png Dans cette configuration, en supposant qu’un seul élève est contaminé, il faudrait vacciner à plus de 85% pour avoir une chance de ne pas avoir une nouvelle contamination. Ce résultat ne change pas si la taille de la matrice change. Il change en revanche en fonction du nombre de tirages, ici 10. Pour n’avoir aucune propagation, il faut entourer la personne contaminé de 4 personnes vaccinées. :math:`p` est la probabilité pour une personné d’être vaccinée (ou immunisée si la vaccination est faite sur 100% des personnes). :math:`p^4` est la probabilité d’avoir 4 personnes vaccinées. :math:`q=1-p^4` est la probabilité qu’une personne ne soient pas immunisées parmi les 4. :math:`1 - (1-q)^{10}` est la probabilité que 10 cours de 3h se passent bien sans contamination. .. code:: ipython3 p = 0.9 p4 = 1 - (1 - p) **4 p4 .. parsed-literal:: 0.9999 .. code:: ipython3 df2 = df.copy() df2["P4**100"] = 1 - df2["p"] ** 40 .. code:: ipython3 import matplotlib.pyplot as plt fig, ax = plt.subplots(1, 2, figsize=(10, 4)) df2.plot(x="p", y="prop", ax=ax[0]) ax[0].plot([0.7, 1.], [2, 2], "--") df2.plot(x="p", y="P4**100", ax=ax[1]); .. image:: td_note_2021_37_0.png