diff --git a/module2/exo1/montecarlo.png b/module2/exo1/montecarlo.png new file mode 100644 index 0000000000000000000000000000000000000000..56a96283d90230655c9da5fad73bba733bd7b947 Binary files /dev/null and b/module2/exo1/montecarlo.png differ diff --git a/module2/exo1/toy_document_orgmode_python_fr.org b/module2/exo1/toy_document_orgmode_python_fr.org index c7157ba42216cf2e1d291112bb351ce48811115c..6b548bff4501f7b2da682d6c498223c6f84cf194 100644 --- a/module2/exo1/toy_document_orgmode_python_fr.org +++ b/module2/exo1/toy_document_orgmode_python_fr.org @@ -1,6 +1,6 @@ -#+TITLE: Votre titre -#+AUTHOR: Votre nom -#+DATE: La date du jour +#+TITLE: À propos du calcul de $\pi$ +#+AUTHOR: François Févotte +#+DATE: 2020-04-06 #+LANGUAGE: fr # #+PROPERTY: header-args :eval never-export @@ -11,83 +11,67 @@ #+HTML_HEAD: #+HTML_HEAD: -* Quelques explications +* En demandant à la lib maths -Ceci est un document org-mode avec quelques exemples de code -python. Une fois ouvert dans emacs, ce document peut aisément être -exporté au format HTML, PDF, et Office. Pour plus de détails sur -org-mode vous pouvez consulter https://orgmode.org/guide/. +Mon ordinateur m'indique que $\pi$ vaut /approximativement/ : +#+begin_src python :results value :exports both :session + from math import * + pi +#+end_src -Lorsque vous utiliserez le raccourci =C-c C-e h o=, ce document sera -compilé en html. Tout le code contenu sera ré-exécuté, les résultats -récupérés et inclus dans un document final. Si vous ne souhaitez pas -ré-exécuter tout le code à chaque fois, il vous suffit de supprimer -le # et l'espace qui sont devant le ~#+PROPERTY:~ au début de ce -document. +#+RESULTS: +: 3.141592653589793 -Comme nous vous l'avons montré dans la vidéo, on inclue du code -python de la façon suivante (et on l'exécute en faisant ~C-c C-c~): +* En utilisant la méthode des aiguilles de Buffon -#+begin_src python :results output :exports both -print("Hello world!") -#+end_src +Mais calculé avec la *méthode* des [[https://fr.wikipedia.org/wiki/Aiguille_de_Buffon][aiguilles de Buffon]], on obtiendrait comme +*approximation* : +#+BEGIN_SRC python :results value :exports both :session + import numpy as np + np.random.seed(seed=42) + N = 10000 + x = np.random.uniform(size=N, low=0, high=1) + theta = np.random.uniform(size=N, low=0, high=pi/2) + 2/(sum((x+np.sin(theta))>1)/N) +#+END_SRC #+RESULTS: -: Hello world! - -Voici la même chose, mais avec une session python, donc une -persistance d'un bloc à l'autre (et on l'exécute toujours en faisant -~C-c C-c~). -#+begin_src python :results output :session :exports both -import numpy -x=numpy.linspace(-15,15) -print(x) -#+end_src +: 3.128911138923655 + +* Avec un argument "fréquentiel" de surface + +Sinon, une méthode plus simple à comprendre et ne faisant pas intervenir d'appel +à la fonction sinus se base sur le fait que si $X \sim U(0,1)$ et $Y \sim +U(0,1)$ alors $P\left[X^2 + Y^2 \leq 1\right] = \pi/4$ (voir [[https://fr.wikipedia.org/wiki/M%C3%A9thode_de_Monte-Carlo#D%C3%A9termination_de_la_valeur_de_%CF%80][méthode de Monte +Carlo sur Wikipedia]]). Le code suivant illustre ce fait : +#+BEGIN_SRC python :results output file :exports both :session :var matplot_lib_filename="montecarlo.png" + import matplotlib.pyplot as plt + + np.random.seed(seed=42) + N = 1000 + x = np.random.uniform(size=N, low=0, high=1) + y = np.random.uniform(size=N, low=0, high=1) + + accept = (x*x+y*y) <= 1 + reject = np.logical_not(accept) + + fig, ax = plt.subplots(1) + ax.scatter(x[accept], y[accept], c='b', alpha=0.2, edgecolor=None) + ax.scatter(x[reject], y[reject], c='r', alpha=0.2, edgecolor=None) + ax.set_aspect('equal') + + plt.savefig(matplot_lib_filename) + print(matplot_lib_filename) +#+END_SRC #+RESULTS: -#+begin_example -[-15. -14.3877551 -13.7755102 -13.16326531 -12.55102041 - -11.93877551 -11.32653061 -10.71428571 -10.10204082 -9.48979592 - -8.87755102 -8.26530612 -7.65306122 -7.04081633 -6.42857143 - -5.81632653 -5.20408163 -4.59183673 -3.97959184 -3.36734694 - -2.75510204 -2.14285714 -1.53061224 -0.91836735 -0.30612245 - 0.30612245 0.91836735 1.53061224 2.14285714 2.75510204 - 3.36734694 3.97959184 4.59183673 5.20408163 5.81632653 - 6.42857143 7.04081633 7.65306122 8.26530612 8.87755102 - 9.48979592 10.10204082 10.71428571 11.32653061 11.93877551 - 12.55102041 13.16326531 13.7755102 14.3877551 15. ] -#+end_example - -Et enfin, voici un exemple de sortie graphique: -#+begin_src python :results output file :session :var matplot_lib_filename="./cosxsx.png" :exports results -import matplotlib.pyplot as plt - -plt.figure(figsize=(10,5)) -plt.plot(x,numpy.cos(x)/x) -plt.tight_layout() - -plt.savefig(matplot_lib_filename) -print(matplot_lib_filename) -#+end_src +[[file:montecarlo.png]] + +Il est alors aisé d'obtenir une approximation (pas terrible) de $\pi$ en +comptant combien de fois, en moyenne, $X^2 + Y^2$ est inférieur à 1 : +#+BEGIN_SRC python :results value :exports both :session + 4*np.mean(accept) +#+END_SRC #+RESULTS: -[[file:./cosxsx.png]] - -Vous remarquerez le paramètre ~:exports results~ qui indique que le code -ne doit pas apparaître dans la version finale du document. Nous vous -recommandons dans le cadre de ce MOOC de ne pas changer ce paramètre -(indiquer ~both~) car l'objectif est que vos analyses de données soient -parfaitement transparentes pour être reproductibles. - -Attention, la figure ainsi générée n'est pas stockée dans le document -org. C'est un fichier ordinaire, ici nommé ~cosxsx.png~. N'oubliez pas -de le committer si vous voulez que votre analyse soit lisible et -compréhensible sur GitLab. - -Enfin, n'oubliez pas que nous vous fournissons dans les ressources de -ce MOOC une configuration avec un certain nombre de raccourcis -claviers permettant de créer rapidement les blocs de code python (en -faisant ~
#+HTML_HEAD: -* Quelques explications -Ceci est un document org-mode avec quelques exemples de code -python. Une fois ouvert dans emacs, ce document peut aisément être -exporté au format HTML, PDF, et Office. Pour plus de détails sur -org-mode vous pouvez consulter https://orgmode.org/guide/. +#+BEGIN_SRC python :exports none :session + import numpy +#+END_SRC -Lorsque vous utiliserez le raccourci =C-c C-e h o=, ce document sera -compilé en html. Tout le code contenu sera ré-exécuté, les résultats -récupérés et inclus dans un document final. Si vous ne souhaitez pas -ré-exécuter tout le code à chaque fois, il vous suffit de supprimer -le # et l'espace qui sont devant le ~#+PROPERTY:~ au début de ce -document. +* Les données -Comme nous vous l'avons montré dans la vidéo, on inclue du code -python de la façon suivante (et on l'exécute en faisant ~C-c C-c~): - -#+begin_src python :results output :exports both -print("Hello world!") +#+begin_src python :results output :exports code :session + x = [14.0, 7.6, 11.2, 12.8, 12.5, 9.9, 14.9, 9.4, 16.9, 10.2, 14.9, 18.1, 7.3, + 9.8, 10.9,12.2, 9.9, 2.9, 2.8, 15.4, 15.7, 9.7, 13.1, 13.2, 12.3, 11.7, 16.0, + 12.4, 17.9, 12.2, 16.2, 18.7, 8.9, 11.9, 12.1, 14.6, 12.1, 4.7, 3.9, 16.9, 16.8, + 11.3, 14.4, 15.7, 14.0, 13.6, 18.0, 13.6, 19.9, 13.7, 17.0, 20.5, 9.9, 12.5, + 13.2, 16.1, 13.5, 6.3, 6.4, 17.6, 19.1, 12.8, 15.5, 16.3, 15.2, 14.6, 19.1, + 14.4, 21.4, 15.1, 19.6, 21.7, 11.3, 15.0, 14.3, 16.8, 14.0, 6.8, 8.2, 19.9, + 20.4, 14.6, 16.4, 18.7, 16.8, 15.8, 20.4, 15.8, 22.4, 16.2, 20.3, 23.4, 12.1, + 15.5, 15.4, 18.4, 15.7, 10.2, 8.9, 21.0] #+end_src #+RESULTS: -: Hello world! - -Voici la même chose, mais avec une session python, donc une -persistance d'un bloc à l'autre (et on l'exécute toujours en faisant -~C-c C-c~). -#+begin_src python :results output :session :exports both -import numpy -x=numpy.linspace(-15,15) -print(x) -#+end_src + + +* Quelques statistiques + +** Moyenne + +#+BEGIN_SRC python :results value :exports both :session + "%.2f" % numpy.mean(x) +#+END_SRC #+RESULTS: -#+begin_example -[-15. -14.3877551 -13.7755102 -13.16326531 -12.55102041 - -11.93877551 -11.32653061 -10.71428571 -10.10204082 -9.48979592 - -8.87755102 -8.26530612 -7.65306122 -7.04081633 -6.42857143 - -5.81632653 -5.20408163 -4.59183673 -3.97959184 -3.36734694 - -2.75510204 -2.14285714 -1.53061224 -0.91836735 -0.30612245 - 0.30612245 0.91836735 1.53061224 2.14285714 2.75510204 - 3.36734694 3.97959184 4.59183673 5.20408163 5.81632653 - 6.42857143 7.04081633 7.65306122 8.26530612 8.87755102 - 9.48979592 10.10204082 10.71428571 11.32653061 11.93877551 - 12.55102041 13.16326531 13.7755102 14.3877551 15. ] -#+end_example - -Et enfin, voici un exemple de sortie graphique: -#+begin_src python :results output file :session :var matplot_lib_filename="./cosxsx.png" :exports results -import matplotlib.pyplot as plt - -plt.figure(figsize=(10,5)) -plt.plot(x,numpy.cos(x)/x) -plt.tight_layout() - -plt.savefig(matplot_lib_filename) -print(matplot_lib_filename) -#+end_src +: 14.11 + + +** Ecart-type + +#+BEGIN_SRC python :results value :exports both :session + "%.2f" % numpy.std(x, ddof=1) +#+END_SRC #+RESULTS: -[[file:./cosxsx.png]] - -Vous remarquerez le paramètre ~:exports results~ qui indique que le code -ne doit pas apparaître dans la version finale du document. Nous vous -recommandons dans le cadre de ce MOOC de ne pas changer ce paramètre -(indiquer ~both~) car l'objectif est que vos analyses de données soient -parfaitement transparentes pour être reproductibles. - -Attention, la figure ainsi générée n'est pas stockée dans le document -org. C'est un fichier ordinaire, ici nommé ~cosxsx.png~. N'oubliez pas -de le committer si vous voulez que votre analyse soit lisible et -compréhensible sur GitLab. - -Enfin, n'oubliez pas que nous vous fournissons dans les ressources de -ce MOOC une configuration avec un certain nombre de raccourcis -claviers permettant de créer rapidement les blocs de code python (en -faisant ~
#+HTML_HEAD: -* Quelques explications - -Ceci est un document org-mode avec quelques exemples de code -python. Une fois ouvert dans emacs, ce document peut aisément être -exporté au format HTML, PDF, et Office. Pour plus de détails sur -org-mode vous pouvez consulter https://orgmode.org/guide/. - -Lorsque vous utiliserez le raccourci =C-c C-e h o=, ce document sera -compilé en html. Tout le code contenu sera ré-exécuté, les résultats -récupérés et inclus dans un document final. Si vous ne souhaitez pas -ré-exécuter tout le code à chaque fois, il vous suffit de supprimer -le # et l'espace qui sont devant le ~#+PROPERTY:~ au début de ce -document. - -Comme nous vous l'avons montré dans la vidéo, on inclue du code -python de la façon suivante (et on l'exécute en faisant ~C-c C-c~): - -#+begin_src python :results output :exports both -print("Hello world!") -#+end_src +#+BEGIN_SRC python :exports none :session + import numpy + import matplotlib.pyplot as plt +#+END_SRC #+RESULTS: -: Hello world! -Voici la même chose, mais avec une session python, donc une -persistance d'un bloc à l'autre (et on l'exécute toujours en faisant -~C-c C-c~). -#+begin_src python :results output :session :exports both -import numpy -x=numpy.linspace(-15,15) -print(x) +* Les données + +#+begin_src python :results output :exports code :session + x = [14.0, 7.6, 11.2, 12.8, 12.5, 9.9, 14.9, 9.4, 16.9, 10.2, 14.9, 18.1, 7.3, + 9.8, 10.9,12.2, 9.9, 2.9, 2.8, 15.4, 15.7, 9.7, 13.1, 13.2, 12.3, 11.7, 16.0, + 12.4, 17.9, 12.2, 16.2, 18.7, 8.9, 11.9, 12.1, 14.6, 12.1, 4.7, 3.9, 16.9, 16.8, + 11.3, 14.4, 15.7, 14.0, 13.6, 18.0, 13.6, 19.9, 13.7, 17.0, 20.5, 9.9, 12.5, + 13.2, 16.1, 13.5, 6.3, 6.4, 17.6, 19.1, 12.8, 15.5, 16.3, 15.2, 14.6, 19.1, + 14.4, 21.4, 15.1, 19.6, 21.7, 11.3, 15.0, 14.3, 16.8, 14.0, 6.8, 8.2, 19.9, + 20.4, 14.6, 16.4, 18.7, 16.8, 15.8, 20.4, 15.8, 22.4, 16.2, 20.3, 23.4, 12.1, + 15.5, 15.4, 18.4, 15.7, 10.2, 8.9, 21.0] #+end_src #+RESULTS: -#+begin_example -[-15. -14.3877551 -13.7755102 -13.16326531 -12.55102041 - -11.93877551 -11.32653061 -10.71428571 -10.10204082 -9.48979592 - -8.87755102 -8.26530612 -7.65306122 -7.04081633 -6.42857143 - -5.81632653 -5.20408163 -4.59183673 -3.97959184 -3.36734694 - -2.75510204 -2.14285714 -1.53061224 -0.91836735 -0.30612245 - 0.30612245 0.91836735 1.53061224 2.14285714 2.75510204 - 3.36734694 3.97959184 4.59183673 5.20408163 5.81632653 - 6.42857143 7.04081633 7.65306122 8.26530612 8.87755102 - 9.48979592 10.10204082 10.71428571 11.32653061 11.93877551 - 12.55102041 13.16326531 13.7755102 14.3877551 15. ] -#+end_example -Et enfin, voici un exemple de sortie graphique: -#+begin_src python :results output file :session :var matplot_lib_filename="./cosxsx.png" :exports results -import matplotlib.pyplot as plt +* Plot +#+begin_src python :results output file :session :var matplot_lib_filename="./plot.png" :exports results plt.figure(figsize=(10,5)) -plt.plot(x,numpy.cos(x)/x) +plt.plot(x) plt.tight_layout() plt.savefig(matplot_lib_filename) @@ -71,23 +45,19 @@ print(matplot_lib_filename) #+end_src #+RESULTS: -[[file:./cosxsx.png]] +[[file:./plot.png]] + -Vous remarquerez le paramètre ~:exports results~ qui indique que le code -ne doit pas apparaître dans la version finale du document. Nous vous -recommandons dans le cadre de ce MOOC de ne pas changer ce paramètre -(indiquer ~both~) car l'objectif est que vos analyses de données soient -parfaitement transparentes pour être reproductibles. +* Histogramme -Attention, la figure ainsi générée n'est pas stockée dans le document -org. C'est un fichier ordinaire, ici nommé ~cosxsx.png~. N'oubliez pas -de le committer si vous voulez que votre analyse soit lisible et -compréhensible sur GitLab. +#+begin_src python :results output file :session :var matplot_lib_filename="./hist.png" :exports results +plt.figure(figsize=(10,5)) +plt.hist(x) +plt.tight_layout() -Enfin, n'oubliez pas que nous vous fournissons dans les ressources de -ce MOOC une configuration avec un certain nombre de raccourcis -claviers permettant de créer rapidement les blocs de code python (en -faisant ~
#+HTML_HEAD: @@ -11,83 +10,223 @@ #+HTML_HEAD: #+HTML_HEAD: -* Quelques explications -Ceci est un document org-mode avec quelques exemples de code -python. Une fois ouvert dans emacs, ce document peut aisément être -exporté au format HTML, PDF, et Office. Pour plus de détails sur -org-mode vous pouvez consulter https://orgmode.org/guide/. +* Identification du système -Lorsque vous utiliserez le raccourci =C-c C-e h o=, ce document sera -compilé en html. Tout le code contenu sera ré-exécuté, les résultats -récupérés et inclus dans un document final. Si vous ne souhaitez pas -ré-exécuter tout le code à chaque fois, il vous suffit de supprimer -le # et l'espace qui sont devant le ~#+PROPERTY:~ au début de ce -document. +#+BEGIN_SRC shell :exports both + uname -a +#+END_SRC -Comme nous vous l'avons montré dans la vidéo, on inclue du code -python de la façon suivante (et on l'exécute en faisant ~C-c C-c~): +#+RESULTS: +: Linux aluminium 4.14.0-3-amd64 #1 SMP Debian 4.14.17-1 (2018-02-14) x86_64 GNU/Linux -#+begin_src python :results output :exports both -print("Hello world!") -#+end_src +#+BEGIN_SRC python :results value :session :exports both + import sys + sys.version +#+END_SRC #+RESULTS: -: Hello world! +: 3.7.3rc1 (default, Mar 13 2019, 11:01:15) +: [GCC 8.3.0] + + +* Récupération des données -Voici la même chose, mais avec une session python, donc une -persistance d'un bloc à l'autre (et on l'exécute toujours en faisant -~C-c C-c~). -#+begin_src python :results output :session :exports both -import numpy -x=numpy.linspace(-15,15) -print(x) -#+end_src +On prend les données distribuées par [[https://ourworldindata.org][Our World In Data]], elles-mêmes issues de +l'[[https://www.ecdc.europa.eu/][ECDC]] : + +#+BEGIN_SRC python :results output :session :exports none + import numpy as np + import matplotlib.pyplot as plt + import pandas as pd + plt.style.use('seaborn') +#+END_SRC + +#+RESULTS: + +#+BEGIN_SRC python :results value :session :exports both + url = "https://covid.ourworldindata.org/data/ecdc/full_data.csv" + data_raw = pd.read_csv(url, index_col='date') + data_raw.tail(3) +#+END_SRC + +#+RESULTS: +: location new_cases new_deaths total_cases total_deaths +: date +: 2020-04-06 Zimbabwe 0 0 9 1 +: 2020-04-07 Zimbabwe 0 0 9 1 +: 2020-04-08 Zimbabwe 1 0 10 1 + + +* Evolution comparée de l'épidémie dans quelques pays + +Je préfère fonder l'analyse sur le nombre de décès, qui me semble plus fiable et +comparable d'un pays à l'autre que, par exemple, le nombre de cas dépistés. Afin +de recaler l'axe temporel, on positionne $t=0$ au moment où le nombre total de +morts dépasse 100. + +La logique est la suivante : on considère le nombre $n^c$ (de décès ou de cas) +dans un pays $c$. On suppose que l'évolution de $n^c$ en fonction du temps $t$ +suit une loi de croissance exponentielle (au moins dans une phase initiale). On +s'attend ainsi à avoir une expression du type : +\[ + n^c(t) = n^c_0 \, e^{\lambda^c \left(t-t^c_0\right)}. +\] +Dans cette expression, la valeur initiale $n^c_0 = n^c(t^c_0)$ est liée à la +taille de la population considérée. La constante de temps $\lambda^c$ en revanche, +n'a pas vraiment de raison de l'être : elle reflète simplement le rythme auquel +l'épidémie se propage. + +Si l'on se donne un seuil $\tau$ unique, on peut fixer pour chaque pays +l'instant initial $t^c_0$ de sorte que +\[ + n^c(t^c_0) = \tau. +\] +En effectuant un recalage de l'instant initial, c'est à dire en effectuant le +changement de variable temporelle $\tilde{t}^c=t-t^c_0$, on obtient +\[ + n^c(\tilde{t}^c) = \tau \, e^{\lambda^c\,\tilde{t}^c}, +\] +ce qui signifie que tous les pays deviennent comparables, indépendamment de leur +taille : on s'est ramené à comparer les constantes de temps $\lambda^c$ entre +elles. + +En pratique, on prend +\[ + t^c_0 = \min \left\{t ; n^c(t) > \tau \right\}. +\] + +#+BEGIN_SRC python :results value :exports both :session + # Parameters + column = "total_deaths" + threshold = 100 + countries = ["France", "Italy", "Spain", "Belgium", "Germany", "Switzerland"] + + total = {} + for country in countries: + total[country] = data_raw.query(f"location == '{country}' & {column} > {threshold}") + + total[countries[0]] +#+END_SRC #+RESULTS: #+begin_example -[-15. -14.3877551 -13.7755102 -13.16326531 -12.55102041 - -11.93877551 -11.32653061 -10.71428571 -10.10204082 -9.48979592 - -8.87755102 -8.26530612 -7.65306122 -7.04081633 -6.42857143 - -5.81632653 -5.20408163 -4.59183673 -3.97959184 -3.36734694 - -2.75510204 -2.14285714 -1.53061224 -0.91836735 -0.30612245 - 0.30612245 0.91836735 1.53061224 2.14285714 2.75510204 - 3.36734694 3.97959184 4.59183673 5.20408163 5.81632653 - 6.42857143 7.04081633 7.65306122 8.26530612 8.87755102 - 9.48979592 10.10204082 10.71428571 11.32653061 11.93877551 - 12.55102041 13.16326531 13.7755102 14.3877551 15. ] + location new_cases new_deaths total_cases total_deaths +date +2020-03-16 France 924 36 5423 127 +2020-03-17 France 1210 21 6633 148 +2020-03-18 France 1097 27 7730 175 +2020-03-19 France 1404 69 9134 244 +2020-03-20 France 1861 128 10995 372 +2020-03-21 France 1617 78 12612 450 +2020-03-22 France 1847 112 14459 562 +2020-03-23 France 1559 112 16018 674 +2020-03-24 France 3838 186 19856 860 +2020-03-25 France 2446 240 22302 1100 +2020-03-26 France 2931 231 25233 1331 +2020-03-27 France 3922 365 29155 1696 +2020-03-28 France 3809 299 32964 1995 +2020-03-29 France 4611 319 37575 2314 +2020-03-30 France 2599 292 40174 2606 +2020-03-31 France 4376 418 44550 3024 +2020-04-01 France 7578 499 52128 3523 +2020-04-02 France 4861 509 56989 4032 +2020-04-03 France 2116 471 59105 4503 +2020-04-04 France 5233 2004 64338 6507 +2020-04-05 France 4267 1053 68605 7560 +2020-04-06 France 1873 518 70478 8078 +2020-04-07 France 3912 833 74390 8911 +2020-04-08 France 3777 1417 78167 10328 #+end_example -Et enfin, voici un exemple de sortie graphique: -#+begin_src python :results output file :session :var matplot_lib_filename="./cosxsx.png" :exports results -import matplotlib.pyplot as plt +#+BEGIN_SRC python :session :results output file :exports results :var pltfile="total.png" + plt.figure(figsize=(10,5)) + for country in total.keys(): + x = total[country]['total_deaths'] + plt.plot(np.arange(len(x)), x, label=country, linewidth=3) + plt.xlabel(f"Number of days since {column} > {threshold}") + plt.yscale("log") + plt.legend(fontsize='x-large') + plt.tight_layout() + plt.savefig(pltfile) + plt.close() + print(pltfile) +#+END_SRC + +#+RESULTS: +[[file:total.png]] -plt.figure(figsize=(10,5)) -plt.plot(x,numpy.cos(x)/x) -plt.tight_layout() -plt.savefig(matplot_lib_filename) -print(matplot_lib_filename) -#+end_src +* Vers un passage du pic ? + +On s'intéresse ici à l'atteinte d'un pic du nombre de décès quotidiens liés à +Covid-19. Là encore, l'axe temporel est recalé afin de permettre une +comparaison entre pays : $t=0$ est positionné au moment où le nombre quotidien +de morts dépasse 15. + +#+BEGIN_SRC python :session :results value :exports both + threshold = 15 + + daily = {} + for country in countries: + daily[country] = data_raw.query(f"location == '{country}' & new_deaths > {threshold}") + + daily[countries[0]] +#+END_SRC + +#+RESULTS: +#+begin_example + location new_cases new_deaths total_cases total_deaths +date +2020-03-14 France 785 18 3661 79 +2020-03-16 France 924 36 5423 127 +2020-03-17 France 1210 21 6633 148 +2020-03-18 France 1097 27 7730 175 +2020-03-19 France 1404 69 9134 244 +2020-03-20 France 1861 128 10995 372 +2020-03-21 France 1617 78 12612 450 +2020-03-22 France 1847 112 14459 562 +2020-03-23 France 1559 112 16018 674 +2020-03-24 France 3838 186 19856 860 +2020-03-25 France 2446 240 22302 1100 +2020-03-26 France 2931 231 25233 1331 +2020-03-27 France 3922 365 29155 1696 +2020-03-28 France 3809 299 32964 1995 +2020-03-29 France 4611 319 37575 2314 +2020-03-30 France 2599 292 40174 2606 +2020-03-31 France 4376 418 44550 3024 +2020-04-01 France 7578 499 52128 3523 +2020-04-02 France 4861 509 56989 4032 +2020-04-03 France 2116 471 59105 4503 +2020-04-04 France 5233 2004 64338 6507 +2020-04-05 France 4267 1053 68605 7560 +2020-04-06 France 1873 518 70478 8078 +2020-04-07 France 3912 833 74390 8911 +2020-04-08 France 3777 1417 78167 10328 +#+end_example + +Le pic observé pour la France au 04/04 (jour 20 sur le graphique) correspond à la prise en compte +instantanée de tous les décès en EHPAD, dont le décompte est disponible à partir +du 01/02 (/cf./ [[https://www.gouvernement.fr/info-coronavirus/carte-et-donnees][données du gouvernement]]) + +#+BEGIN_SRC python :session :results output file :exports results :var pltfile="daily.png" + plt.figure(figsize=(10,5)) + for country in daily.keys(): + x = daily[country]['new_deaths'] + plt.plot(np.arange(len(x)), x, label=country, linewidth=3) + plt.xlabel(f"Number of days since daily deaths > {threshold}") + plt.yscale("log") + plt.legend(fontsize='x-large') + plt.tight_layout() + plt.savefig(pltfile) + plt.close() + print(pltfile) +#+END_SRC #+RESULTS: -[[file:./cosxsx.png]] - -Vous remarquerez le paramètre ~:exports results~ qui indique que le code -ne doit pas apparaître dans la version finale du document. Nous vous -recommandons dans le cadre de ce MOOC de ne pas changer ce paramètre -(indiquer ~both~) car l'objectif est que vos analyses de données soient -parfaitement transparentes pour être reproductibles. - -Attention, la figure ainsi générée n'est pas stockée dans le document -org. C'est un fichier ordinaire, ici nommé ~cosxsx.png~. N'oubliez pas -de le committer si vous voulez que votre analyse soit lisible et -compréhensible sur GitLab. - -Enfin, n'oubliez pas que nous vous fournissons dans les ressources de -ce MOOC une configuration avec un certain nombre de raccourcis -claviers permettant de créer rapidement les blocs de code python (en -faisant ~0] -data -#+end_src - -#+RESULTS: -: Date Count Temperature Pressure Malfunction -: 1 11/12/81 6 70 50 1 -: 8 2/03/84 6 57 200 1 -: 9 4/06/84 6 63 200 1 -: 10 8/30/84 6 70 200 1 -: 13 1/24/85 6 53 200 2 -: 20 10/30/85 6 75 200 2 -: 22 1/12/86 6 58 200 1 +# Les vols où aucun incident n'est relevé n'apportant aucune information +# sur l'influence de la température ou de la pression sur les +# dysfonctionnements, nous nous concentrons sur les expériences où au +# moins un joint a été défectueux. + +# #+begin_src python :results value :session *python* :exports both +# data = data[data.Malfunction>0] +# data +# #+end_src + +# #+RESULTS: +# #+begin_example +# Date Count Temperature Pressure Malfunction +# 0 4/12/81 6 66 50 0 +# 1 11/12/81 6 70 50 1 +# 2 3/22/82 6 69 50 0 +# 3 11/11/82 6 68 50 0 +# 4 4/04/83 6 67 50 0 +# 5 6/18/82 6 72 50 0 +# 6 8/30/83 6 73 100 0 +# 7 11/28/83 6 70 100 0 +# 8 2/03/84 6 57 200 1 +# 9 4/06/84 6 63 200 1 +# 10 8/30/84 6 70 200 1 +# 11 10/05/84 6 78 200 0 +# 12 11/08/84 6 67 200 0 +# 13 1/24/85 6 53 200 2 +# 14 4/12/85 6 67 200 0 +# 15 4/29/85 6 75 200 0 +# 16 6/17/85 6 70 200 0 +# 17 7/29/85 6 81 200 0 +# 18 8/27/85 6 76 200 0 +# 19 10/03/85 6 79 200 0 +# 20 10/30/85 6 75 200 2 +# 21 11/26/85 6 76 200 0 +# 22 1/12/86 6 58 200 1 +# #+end_example Très bien, nous avons une variabilité de température importante mais la pression est quasiment toujours égale à 200, ce qui devrait @@ -145,26 +164,27 @@ logmodel.summary() #+begin_example Generalized Linear Model Regression Results ============================================================================== -Dep. Variable: Frequency No. Observations: 7 -Model: GLM Df Residuals: 5 +Dep. Variable: Frequency No. Observations: 23 +Model: GLM Df Residuals: 21 Model Family: Binomial Df Model: 1 -Link Function: logit Scale: 1.0 -Method: IRLS Log-Likelihood: -3.6370 -Date: Fri, 20 Jul 2018 Deviance: 3.3763 -Time: 16:56:08 Pearson chi2: 0.236 -No. Iterations: 5 +Link Function: logit Scale: 1.0000 +Method: IRLS Log-Likelihood: -3.9210 +Date: Thu, 09 Apr 2020 Deviance: 3.0144 +Time: 12:50:13 Pearson chi2: 5.00 +No. Iterations: 6 +Covariance Type: nonrobust =============================================================================== coef std err z P>|z| [0.025 0.975] ------------------------------------------------------------------------------- -Intercept -1.3895 7.828 -0.178 0.859 -16.732 13.953 -Temperature 0.0014 0.122 0.012 0.991 -0.238 0.240 +Intercept 5.0850 7.477 0.680 0.496 -9.570 19.740 +Temperature -0.1156 0.115 -1.004 0.316 -0.341 0.110 =============================================================================== #+end_example -L'estimateur le plus probable du paramètre de température est 0.0014 -et l'erreur standard de cet estimateur est de 0.122, autrement dit on -ne peut pas distinguer d'impact particulier et il faut prendre nos -estimations avec des pincettes. +# L'estimateur le plus probable du paramètre de température est 0.0014 +# et l'erreur standard de cet estimateur est de 0.122, autrement dit on +# ne peut pas distinguer d'impact particulier et il faut prendre nos +# estimations avec des pincettes. * Estimation de la probabilité de dysfonctionnant des joints toriques La température prévue le jour du décollage est de 31°F. Essayons @@ -187,36 +207,43 @@ print(matplot_lib_filename) #+RESULTS: [[file:proba_estimate_python.png]] -Comme on pouvait s'attendre au vu des données initiales, la -température n'a pas d'impact notable sur la probabilité d'échec des -joints toriques. Elle sera d'environ 0.2, comme dans les essais -précédents où nous il y a eu défaillance d'au moins un joint. Revenons -à l'ensemble des données initiales pour estimer la probabilité de -défaillance d'un joint: - -#+begin_src python :results output :session *python* :exports both -data = pd.read_csv("shuttle.csv") -print(np.sum(data.Malfunction)/np.sum(data.Count)) -#+end_src - -#+RESULTS: -: 0.06521739130434782 - -Cette probabilité est donc d'environ $p=0.065$, sachant qu'il existe -un joint primaire un joint secondaire sur chacune des trois parties du -lançeur, la probabilité de défaillance des deux joints d'un lançeur -est de $p^2 \approx 0.00425$. La probabilité de défaillance d'un des -lançeur est donc de $1-(1-p^2)^3 \approx 1.2%$. Ça serait vraiment -pas de chance... Tout est sous contrôle, le décollage peut donc avoir -lieu demain comme prévu. - -Seulement, le lendemain, la navette Challenger explosera et emportera -avec elle ses sept membres d'équipages. L'opinion publique est -fortement touchée et lors de l'enquête qui suivra, la fiabilité des -joints toriques sera directement mise en cause. Au delà des problèmes -de communication interne à la NASA qui sont pour beaucoup dans ce -fiasco, l'analyse précédente comporte (au moins) un petit -problème... Saurez-vous le trouver ? Vous êtes libre de modifier cette -analyse et de regarder ce jeu de données sous tous les angles afin -d'expliquer ce qui ne va pas. - +D'après cette analyse, la probabilité d'échec des joints toriques serait +supérieure à 80% à pour une température de 31°F. + +# Comme on pouvait s'attendre au vu des données initiales, la +# température n'a pas d'impact notable sur la probabilité d'échec des +# joints toriques. Elle sera d'environ 0.2, comme dans les essais +# précédents où nous il y a eu défaillance d'au moins un joint. Revenons +# à l'ensemble des données initiales pour estimer la probabilité de +# défaillance d'un joint: + +# #+begin_src python :results output :session *python* :exports both +# data = pd.read_csv("shuttle.csv") +# print(np.sum(data.Malfunction)/np.sum(data.Count)) +# #+end_src + +# #+RESULTS: +# : 0.06521739130434782 + +# Cette probabilité est donc d'environ $p=0.065$, sachant qu'il existe +# un joint primaire un joint secondaire sur chacune des trois parties du +# lançeur, la probabilité de défaillance des deux joints d'un lançeur +# est de $p^2 \approx 0.00425$. La probabilité de défaillance d'un des +# lançeur est donc de $1-(1-p^2)^3 \approx 1.2%$. Ça serait vraiment +# pas de chance... Tout est sous contrôle, le décollage peut donc avoir +# lieu demain comme prévu. + +# Seulement, le lendemain, la navette Challenger explosera et emportera +# avec elle ses sept membres d'équipages. L'opinion publique est +# fortement touchée et lors de l'enquête qui suivra, la fiabilité des +# joints toriques sera directement mise en cause. Au delà des problèmes +# de communication interne à la NASA qui sont pour beaucoup dans ce +# fiasco, l'analyse précédente comporte (au moins) un petit +# problème... Saurez-vous le trouver ? Vous êtes libre de modifier cette +# analyse et de regarder ce jeu de données sous tous les angles afin +# d'expliquer ce qui ne va pas. + +# Local Variables: +# org-babel-python-command: "/home/francois/.julia/conda/3/bin/python" +# org-confirm-babel-evaluate: nil +# End: diff --git a/module2/exo5/freq_temp_python.png b/module2/exo5/freq_temp_python.png index 93cb9e626441d23f6dff59ed252d7b14eb37abdb..dc6464064ba17d067083743f62aacdcf9980fc54 100644 Binary files a/module2/exo5/freq_temp_python.png and b/module2/exo5/freq_temp_python.png differ diff --git a/module2/exo5/proba_estimate_python.png b/module2/exo5/proba_estimate_python.png index 77fc4b275dd8815b1ab91cd3b67b1beb93e00748..ead96326d1a0082ba7106fe789cc9916cfb7051c 100644 Binary files a/module2/exo5/proba_estimate_python.png and b/module2/exo5/proba_estimate_python.png differ