Commit 19885907 authored by Samuel MEYNARD's avatar Samuel MEYNARD

Début d'étude ici jusqu'au point 2

parent ef18fe6a
......@@ -12,10 +12,65 @@
* En demandant à la lib maths
Mon ordinateur m'indique que π vaut approximativement:
#+begin_src python :results output :exports both
from math import *
Mon ordinateur m'indique que $\pi$ vaut /approximativement/:
#+begin_src python :results value :session *python* :exports both
from math import pi
pi
#+end_src
#+RESULTS:
: 3.141592653589793
* En utilisant la méthode des aiguilles de Buffon
Mais calculé avec la *méthode* des _aiguilles de Buffon_, on obtiendrait comme *approximation* :
#+begin_src python :results value :session *python* :exports both
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:
: 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∼U(0,1) et Y∼U(0,1) alors P[X2+Y2≤1]=\pi/4$ (voir méthode de Monte Carlo sur Wikipedia). Le code suivant illustre ce fait :
#+begin_src python :results output file :var matplot_lib_filename="figure_pi_mc2.png" :exports both :session *python*y
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
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 output :session *python* :exports both
4*np.mean(accept)
#+end_src
#+RESULTS:
: 3.112
:
:
......@@ -60,7 +60,7 @@ print(x)
#+end_example
Finally, an example for graphical output:
#+begin_src python :results output file :session :var matplot_lib_filename="./cosxsx.png" :exports results
#+begin_src python :results value :session :var matplot_lib_filename="./cosxsx.png" :exports results
import matplotlib.pyplot as plt
plt.figure(figsize=(10,5))
......@@ -72,7 +72,8 @@ print(matplot_lib_filename)
#+end_src
#+RESULTS:
[[file:./cosxsx.png]]
[[./cosxsx.png]]
Note the parameter ~:exports results~, which indicates that the code
will not appear in the exported document. We recommend that in the
......
#+TITLE: Votre titre
#+TITLE: Exo2
#+AUTHOR: Votre nom
#+DATE: La date du jour
#+LANGUAGE: fr
......@@ -11,83 +11,108 @@
#+HTML_HEAD: <script type="text/javascript" src="http://www.pirilampo.org/styles/lib/js/jquery.stickytableheaders.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="http://www.pirilampo.org/styles/readtheorg/js/readtheorg.js"></script>
* 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!")
* Donnée
#+begin_src python :results value :exports both :session
data = [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]
data
#+end_src
#+RESULTS:
: Hello world!
| 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 |
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)
#+begin_src python :results value :session :var matplot_lib_filename="fig1.png" :exports both
from numpy import mean, min, max, median, std
import matplotlib
import matplotlib.pyplot as plt
d_mean = mean(data)
d_min = min(data)
d_max = max(data)
d_median = median(data)
d_std = std(data, ddof=1)
f"Mean : {d_mean}\nMin: {d_min}\nMax: {d_max}\nMediane: {d_median}\nEcart Type: {d_std}"
#+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
: Mean : 14.113000000000001
: Min: 2.8
: Max: 23.4
: Mediane: 14.5
: Ecart Type: 4.334094455301447
#+BEGIN_SRC python :session :results file
import matplotlib
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
fig=plt.figure(figsize=(3,2))
plt.hist(data)
plt.savefig('test.png')
'./test.png'
#+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 ~<p~, ~<P~ ou ~<PP~ suivi de ~Tab~).
Maintenant, à vous de jouer! Vous pouvez effacer toutes ces
informations et les remplacer par votre document computationnel.
[[file:]]
[[file:Traceback (most recent call last):
File "<stdin>", line 7, in <module>
File "<string>", line 4, in <module>
File "/Library/Python/3.7/site-packages/matplotlib/pyplot.py", line 2610, in hist
if data is not None else {}), **kwargs)
File "/Library/Python/3.7/site-packages/matplotlib/__init__.py", line 1565, in inner
return func(ax, *map(sanitize_sequence, args), **kwargs)
File "/Library/Python/3.7/site-packages/matplotlib/axes/_axes.py", line 6559, in hist
cbook._check_in_list(['left', 'mid', 'right'], align=align)
File "/Library/Python/3.7/site-packages/matplotlib/cbook/__init__.py", line 2145, in _check_in_list
.format(v, k, ', '.join(map(repr, values))))
ValueError: 'edge' is not a valid value for align; supported values are 'left', 'mid', 'right'
]]
[[file:Traceback (most recent call last):
File "<stdin>", line 7, in <module>
File "<string>", line 4, in <module>
File "/Library/Python/3.7/site-packages/matplotlib/pyplot.py", line 2610, in hist
if data is not None else {}), **kwargs)
File "/Library/Python/3.7/site-packages/matplotlib/__init__.py", line 1565, in inner
return func(ax, *map(sanitize_sequence, args), **kwargs)
File "/Library/Python/3.7/site-packages/matplotlib/axes/_axes.py", line 6559, in hist
cbook._check_in_list(['left', 'mid', 'right'], align=align)
File "/Library/Python/3.7/site-packages/matplotlib/cbook/__init__.py", line 2145, in _check_in_list
.format(v, k, ', '.join(map(repr, values))))
ValueError: 'edge' is not a valid value for align; supported values are 'left', 'mid', 'right'
]]
File "<stdin>", line 7, in <module>
File "<string>", line 2, in <module>
ModuleNotFoundError: No module named 'hist'
]]
File "<stdin>", line 7, in <module>
File "<string>", line 4, in <module>
TypeError: bar() missing 1 required positional argument: 'height'
]]
File "<stdin>", line 7, in <module>
File "<string>", line 4, in <module>
TypeError: bar() missing 1 required positional argument: 'height'
]]
File "<stdin>", line 7, in <module>
File "<string>", line 4, in <module>
TypeError: bar() missing 1 required positional argument: 'height'
]]
File "<stdin>", line 7, in <module>
File "<string>", line 4, in <module>
TypeError: bar() missing 1 required positional argument: 'height'
]]
File "<stdin>", line 7, in <module>
File "<string>", line 3, in <module>
TypeError: bar() missing 1 required positional argument: 'height'
]]
File "<stdin>", line 7, in <module>
File "<string>", line 3, in <module>
TypeError: ba .r() missing 1 required positional argument: 'height'
]]
File "<stdin>", line 7, in <module>
File "<string>", line 3, in <module>
TypeError: bar() missing 1 required positional argument: 'height'
]]
File "<stdin>", line 7, in <module>
File "<string>", line 3, in <module>
NameError: name 'data' is not defined
]]
File "<stdin>", line 7, in <module>
File "<string>", line 1, in <module>
NameError: name 'plt' is not defined
]]
[[file:./test.png]]
......@@ -171,7 +171,7 @@ La température prévue le jour du décollage est de 31°F. Essayons
d'estimer la probabilité de dysfonctionnement des joints toriques à
cette température à partir du modèle que nous venons de construire:
#+begin_src python :results output file :var matplot_lib_filename="proba_estimate_python.png" :exports both :session *python*
#+begin_src python :results output file :var matplot_lib_filename="proba_estimate_python.png" :exports both :session *python*
import matplotlib.pyplot as plt
data_pred = pd.DataFrame({'Temperature': np.linspace(start=30, stop=90, num=121), 'Intercept': 1})
......
module2/exo5/freq_temp_python.png

12.3 KB | W: | H:

module2/exo5/freq_temp_python.png

12.2 KB | W: | H:

module2/exo5/freq_temp_python.png
module2/exo5/freq_temp_python.png
module2/exo5/freq_temp_python.png
module2/exo5/freq_temp_python.png
  • 2-up
  • Swipe
  • Onion skin
module2/exo5/proba_estimate_python.png

14.3 KB | W: | H:

module2/exo5/proba_estimate_python.png

14.2 KB | W: | H:

module2/exo5/proba_estimate_python.png
module2/exo5/proba_estimate_python.png
module2/exo5/proba_estimate_python.png
module2/exo5/proba_estimate_python.png
  • 2-up
  • Swipe
  • Onion skin
......@@ -27,11 +27,18 @@ if sys.version_info.major < 3 or sys.version_info.minor < 6:
print("Veuillez utiliser Python 3.6 (ou plus) !")
#+END_SRC
#+RESULTS:
: Python 3.8.2 (default, Mar 11 2020, 00:29:50)
: [Clang 11.0.0 (clang-1100.0.33.17)] on darwin
: Type "help", "copyright", "credits" or "license" for more information.
#+BEGIN_SRC emacs-lisp :results output
(unless (featurep 'ob-python)
(print "Veuillez activer python dans org-babel (org-babel-do-languages) !"))
#+END_SRC
#+RESULTS:
** R 3.4
Nous n'utilisons que des fonctionnalités de base du langage R, une version antérieure devrait suffire.
......@@ -40,6 +47,14 @@ Nous n'utilisons que des fonctionnalités de base du langage R, une version ant
(print "Veuillez activer R dans org-babel (org-babel-do-languages) !"))
#+END_SRC
#+RESULTS:
#+begin_src R :results output graphics :file (org-babel-temp-file "figure" ".png") :exports both :width 600 :height 400 :session *R*
test
#+end_src
#+RESULTS:
* Préparation des données
Les données de l'incidence du syndrome grippal sont disponibles du site Web du [[http://www.sentiweb.fr/][Réseau Sentinelles]]. Nous les récupérons sous forme d'un fichier en format CSV dont chaque ligne correspond à une semaine de la période d'observation. Nous téléchargeons toujours le jeu de données complet (rien d'autre n'est proposé), qui commence en 1984 et se termine avec une semaine récente. L'URL est:
......@@ -80,6 +95,13 @@ Regardons ce que nous avons obtenu:
table[:5]
#+END_SRC
#+RESULTS:
| week | indicator | inc | inc_low | inc_up | inc100 | inc100_low | inc100_up | geo_insee | geo_name |
| 202011 | 3 | 101704 | 93652 | 109756 | 154 | 142 | 166 | FR | France |
| 202010 | 3 | 104977 | 96650 | 113304 | 159 | 146 | 172 | FR | France |
| 202009 | 3 | 110696 | 102066 | 119326 | 168 | 155 | 181 | FR | France |
| 202008 | 3 | 143753 | 133984 | 153522 | 218 | 203 | 233 | FR | France |
** Recherche de données manquantes
Il y a malheureusement beaucoup de façon d'indiquer l'absence d'un point de données. Nous testons ici seulement pour la présence de champs vides. Il faudrait aussi rechercher des valeurs non-numériques dans les colonnes à priori numériques. Nous ne le faisons pas ici, mais une vérification ultérieure capterait des telles anomalies.
......@@ -94,6 +116,9 @@ for row in table:
valid_table.append(row)
#+END_SRC
#+RESULTS:
: ['198919', '3', '0', '', '', '0', '', '', 'FR', 'France']
** Extraction des colonnes utilisées
Il y a deux colonnes qui nous intéressent: la première (~"week"~) et la troisième (~"inc"~). Nous vérifions leurs noms dans l'en-tête, que nous effaçons par la suite. Enfin, nous créons un tableau avec les deux colonnes pour le traitement suivant.
#+BEGIN_SRC python :results silent
......@@ -101,7 +126,7 @@ week = [row[0] for row in valid_table]
assert week[0] == 'week'
del week[0]
inc = [row[2] for row in valid_table]
assert inc[0] == 'inc
assert inc[0] == 'inc'
del inc[0]
data = list(zip(week, inc))
#+END_SRC
......@@ -111,6 +136,21 @@ Regardons les premières et les dernières lignes. Nous insérons ~None~ pour in
[('week', 'inc'), None] + data[:5] + [None] + data[-5:]
#+END_SRC
#+RESULTS:
| week | inc |
|--------+--------|
| 202011 | 101704 |
| 202010 | 104977 |
| 202009 | 110696 |
| 202008 | 143753 |
| 202007 | 183610 |
|--------+--------|
| 198448 | 78620 |
| 198447 | 72029 |
| 198446 | 87330 |
| 198445 | 135223 |
| 198444 | 68422 |
** Vérification
Il est toujours prudent de vérifier si les données semblent crédibles. Nous savons que les semaines sont données par six chiffres (quatre pour l'année et deux pour la semaine), et que les incidences sont des nombres entiers positifs.
#+BEGIN_SRC python :results output
......@@ -121,6 +161,8 @@ for week, inc in data:
print("Valeur suspecte dans la colonne 'inc': ", (week, inc))
#+END_SRC
#+RESULTS:
Pas de problème !
** Conversions
......@@ -140,6 +182,21 @@ str_data = [(str(date), str(inc)) for date, inc in converted_data]
[('date', 'inc'), None] + str_data[:5] + [None] + str_data[-5:]
#+END_SRC
#+RESULTS:
| date | inc |
|------------+--------|
| 1984-10-29 | 68422 |
| 1984-11-05 | 135223 |
| 1984-11-12 | 87330 |
| 1984-11-19 | 72029 |
| 1984-11-26 | 78620 |
|------------+--------|
| 2020-02-10 | 183610 |
| 2020-02-17 | 143753 |
| 2020-02-24 | 110696 |
| 2020-03-02 | 104977 |
| 2020-03-09 | 101704 |
** Vérification des dates
Nous faisons encore une vérification: nos dates doivent être séparées d'exactement une semaine, sauf autour du point manquant.
#+BEGIN_SRC python :results output
......@@ -149,6 +206,9 @@ for date1, date2 in zip(dates[:-1], dates[1:]):
print(f"Il y a {date2-date1} entre {date1} et {date2}")
#+END_SRC
#+RESULTS:
: Il y a 14 days, 0:00:00 entre 1989-05-01 et 1989-05-15
** Passage Python -> R
Nous passons au langage R pour inspecter nos données, parce que l'analyse et la préparation de graphiques sont plus concises en R, sans nécessiter aucune bibliothèque supplémentaire.
......@@ -164,17 +224,30 @@ data$date <- as.Date(data$date)
summary(data)
#+END_SRC
#+RESULTS:
:
: date inc
: Min. :1984-10-29 Min. : 0
: 1st Qu.:1993-09-06 1st Qu.: 5018
: Median :2002-07-08 Median : 16002
: Mean :2002-07-07 Mean : 62071
: 3rd Qu.:2011-05-09 3rd Qu.: 50759
: Max. :2020-03-09 Max. :1001824
** Inspection
Regardons enfin à quoi ressemblent nos données !
#+BEGIN_SRC R :results output graphics :file inc-plot.png
<
plot(data, type="l", xlab="Date", ylab="Incidence hebdomadaire")
#+END_SRC
Un zoom sur les dernières années montre mieux la situation des pics en hiver. Le creux des incidences se trouve en été.
#+BEGIN_SRC R :results output graphics :file inc-plot-zoom.png
plot(tail(data, 200), type="l", xlab="Date", ylab="Incidence hebdomadaire")
#+END_SRC
#+RESULTS:
* Étude de l'incidence annuelle
** Calcul de l'incidence annuelle
......@@ -201,6 +274,9 @@ inc_annuelle = data.frame(annee = annees,
head(inc_annuelle)
#+END_SRC
#+RESULTS:
: org_babel_R_eoe
** Inspection
Voici les incidences annuelles en graphique.
#+BEGIN_SRC R :results output graphics :file annual-inc-plot.png
......@@ -213,7 +289,19 @@ Une liste triée par ordre décroissant d'incidence annuelle permet de plus faci
head(inc_annuelle[order(-inc_annuelle$incidence),])
#+END_SRC
#+RESULTS:
: Error in head(inc_annuelle[order(-inc_annuelle$incidence), ]) :
: objet 'inc_annuelle' introuvable
Enfin, un histogramme montre bien que les épidémies fortes, qui touchent environ 10% de la population française, sont assez rares: il y en eu trois au cours des 35 dernières années.
#+BEGIN_SRC R :results output graphics :file annual-inc-hist.png
hist(inc_annuelle$incidence, breaks=10, xlab="Incidence annuelle", ylab="Nb d'observations", main="")
#+END_SRC
#+NAME: test
#+BEGIN_SRC R :file "test1.png" :results output graphics silent
x <- seq(-pi, pi, 0.1)
plot(x, sin(x))
#+END_SRC
[[file:test.svg]]
This diff is collapsed.
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment