Add new file

module2/exo5/exo5_python_fr.org

+#+TITLE: Analyse du risque de défaillance des joints toriques de la navette Challenger
+#+AUTHOR: Arnaud Legrand
+#+LANGUAGE: fr
+
+#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="http://www.pirilampo.org/styles/readtheorg/css/htmlize.css"/>
+#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="http://www.pirilampo.org/styles/readtheorg/css/readtheorg.css"/>
+#+HTML_HEAD: <script src="https://ajax.googleapis.com/ajax/libs/jquery/2.1.3/jquery.min.js"></script>
+#+HTML_HEAD: <script src="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.4/js/bootstrap.min.js"></script>
+#+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>
+
+#+LATEX_HEADER: \usepackage{a4}
+#+LATEX_HEADER: \usepackage[french]{babel}
+
+# #+PROPERTY: header-args  :session  :exports both
+
+Le 27 Janvier 1986, veille du décollage de la navette /Challenger/, eu
+lieu une télé-conférence de trois heures entre les ingénieurs de la
+Morton Thiokol (constructeur d'un des moteurs) et de la NASA. La
+discussion portait principalement sur les conséquences de la
+température prévue au moment du décollage de 31°F (juste en dessous de
+0°C) sur le succès du vol et en particulier sur la performance des
+joints toriques utilisés dans les moteurs. En effet, aucun test
+n'avait été effectué à cette température.
+
+L'étude qui suit reprend donc une partie des analyses effectuées cette
+nuit là et dont l'objectif était d'évaluer l'influence potentielle de
+la température et de la pression à laquelle sont soumis les joints
+toriques sur leur probabilité de dysfonctionnement. Pour cela, nous
+disposons des résultats des expériences réalisées par les ingénieurs
+de la NASA durant les 6 années précédant le lancement de la navette
+Challenger.
+
+
+* Première analyse
+
+** Chargement des données
+Nous commençons donc par charger ces données:
+#+begin_src python :results value :session first :exports both
+import numpy as np
+import pandas as pd
+data = pd.read_csv("shuttle.csv")
+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
+
+Le jeu de données nous indique la date de l'essai, le nombre de joints
+toriques mesurés (il y en a 6 sur le lançeur principal), la
+température (en Fahrenheit) et la pression (en psi), et enfin le
+nombre de dysfonctionnements relevés. 
+
+** Inspection graphique des données
+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 first :exports both
+data = data[data.Malfunction>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
+
+Très bien, nous avons une variabilité de température importante mais
+la pression est quasiment toujours égale à 200, ce qui devrait
+simplifier l'analyse.
+
+Comment la fréquence d'échecs varie-t-elle avec la température ?
+#+begin_src python :results output file :var matplot_lib_filename="freq_temp_python.png" :exports both :session first 
+import matplotlib.pyplot as plt
+
+plt.clf()
+data["Frequency"]=data.Malfunction/data.Count
+data.plot(x="Temperature",y="Frequency",kind="scatter",ylim=[0,1])
+plt.grid(True)
+
+plt.savefig(matplot_lib_filename)
+print(matplot_lib_filename)
+#+end_src
+
+#+RESULTS:
+[[file:freq_temp_python.png]]
+
+À première vue, ce n'est pas flagrant mais bon, essayons quand même
+d'estimer l'impact de la température $t$ sur la probabilité de
+dysfonctionnements d'un joint. 
+
+** Estimation de l'influence de la température
+
+Supposons que chacun des 6 joints toriques est endommagé avec la même
+probabilité et indépendamment des autres et que cette probabilité ne
+dépend que de la température. Si on note $p(t)$ cette probabilité, le
+nombre de joints $D$ dysfonctionnant lorsque l'on effectue le vol à
+température $t$ suit une loi binomiale de paramètre $n=6$ et
+$p=p(t)$. Pour relier $p(t)$ à $t$, on va donc effectuer une
+régression logistique.
+
+#+begin_src python :results value :session first :exports both
+import statsmodels.api as sm
+
+data["Success"]=data.Count-data.Malfunction
+data["Intercept"]=1
+
+
+# logit_model=sm.Logit(data["Frequency"],data[["Intercept","Temperature"]]).fit() 
+logmodel=sm.GLM(data['Frequency'], data[['Intercept','Temperature']], family=sm.families.Binomial(sm.families.links.logit)).fit()
+
+logmodel.summary()
+#+end_src
+
+#+RESULTS:
+#+begin_example
+                 Generalized Linear Model Regression Results                  
+==============================================================================
+Dep. Variable:              Frequency   No. Observations:                    7
+Model:                            GLM   Df Residuals:                        5
+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                                         
+===============================================================================
+                  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
+===============================================================================
+#+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.
+
+** Estimation de la probabilité de dysfonctionnant des joints toriques
+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 first 
+import matplotlib.pyplot as plt
+
+data_pred = pd.DataFrame({'Temperature': np.linspace(start=30, stop=90, num=121), 'Intercept': 1})
+data_pred['Frequency'] = logmodel.predict(data_pred[['Intercept','Temperature']])
+data_pred.plot(x="Temperature",y="Frequency",kind="line",ylim=[0,1])
+plt.scatter(x=data["Temperature"],y=data["Frequency"])
+plt.grid(True)
+
+plt.savefig(matplot_lib_filename)
+print(matplot_lib_filename)
+#+end_src
+
+#+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 first :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.
+
+
+
+* Nouvelle analyse
+
+Nous commençons donc par charger ces données:
+#+begin_src python :results value :session second :exports both
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import statsmodels.api as sm
+
+data = pd.read_csv("shuttle.csv")
+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
+
+
+Tracons, pour ces 3 jeux de données, le nombre de Malfunction en fonction de la température sur le même graphique, une couleur par pression. On le nomme "P_temp_python.png". 
+
+#+begin_src python :results output file :var matplot_lib_filename="P_temp_python.png" :exports both :session second
+
+dataP50 = data[data.Pressure==50]
+dataP100 = data[data.Pressure==100]
+dataP200 = data[data.Pressure==200]
+
+plt.clf()
+dataP50.plot(x="Temperature",y="Malfunction",kind="scatter",ylim=[-1,6],color='blue')
+dataP100.plot(x="Temperature",y="Malfunction",kind="scatter",ylim=[-1,6],color='red',ax=plt.gca())
+dataP200.plot(x="Temperature",y="Malfunction",kind="scatter",ylim=[-1,6],color='green',ax=plt.gca())
+plt.grid(True)
+plt.savefig(matplot_lib_filename)
+print(matplot_lib_filename)
+#+end_src
+
+#+RESULTS:
+[[file:P_temp_python.png]]
+
+
+** Regression logistique
+
+Nous allons maintenant effectuer 4 régressions logistiques
+
+*** Régression logistique sur les données P50
+
+#+begin_src python :results value :session second :exports both
+
+# Filtrer les données et créer une copie
+dataP50 = data[data.Pressure == 50].copy()
+
+# Calcul des colonnes nécessaires
+dataP50["Success"] = dataP50["Count"] - dataP50["Malfunction"]
+dataP50["Intercept"] = 1
+
+# Création du modèle avec une instance de sm.families.links.Logit()
+logmodelP50 = sm.GLM(
+    dataP50['Malfunction'],
+    dataP50[['Intercept', 'Temperature']],
+    family=sm.families.Binomial(link=sm.families.links.Logit())
+).fit()
+
+# Résumé du modèle
+logmodelP50.summary()
+
+# Générer des données pour la prédiction (de 25 à 90)
+data_pred50 = pd.DataFrame({
+    'Temperature': np.linspace(25, 90, 100),
+    'Intercept': 1
+})
+
+# Ajouter les prédictions au DataFrame
+data_pred50['Malfunction_Prob'] = logmodelP50.predict(data_pred50[['Intercept', 'Temperature']])
+
+# Tracer les points observés
+plt.figure()  # Nouvelle figure
+plt.scatter(dataP50['Temperature'], dataP50['Malfunction'], color='blue', label='Observations')
+
+# Tracer la courbe de régression logistique
+plt.plot(data_pred50['Temperature'], data_pred50['Malfunction_Prob'], color='red', label='Régression logistique')
+
+# Ajouter des labels et une légende
+plt.xlabel('Température')
+plt.ylabel('Probabilité de dysfonctionnement')
+plt.title('Régression logistique sur les données P50')
+plt.legend()
+plt.grid(True)
+
+# Créer une image regression_logistique_P50.png
+plt.savefig("regression_logistique_P50.png")
+print("regression_logistique_P50.png")
+#+end_src
+
+#+RESULTS:
+[[file:regression_logistique_P50.png]]
+
+
+*** Régression logistique sur les données P200
+
+#+begin_src python :results value :session second :exports both
+#+begin_src python :results value :session second :exports both
+# Filtrer les données et créer une copie
+dataP200 = data[data.Pressure == 200].copy()
+
+# Calcul des colonnes nécessaires
+dataP200["Success"] = dataP200["Count"] - dataP200["Malfunction"]
+dataP200["Intercept"] = 1
+dataP200["Proportion"] = dataP200["Malfunction"] / dataP200["Count"]
+
+# Création du modèle avec une instance de sm.families.links.Logit()
+logmodelP200 = sm.GLM(
+    dataP200['Proportion'],
+    dataP200[['Intercept', 'Temperature']],
+    family=sm.families.Binomial(link=sm.families.links.Logit())
+).fit()
+
+# Résumé du modèle
+logmodelP200.summary()
+
+# Générer des données pour la prédiction (de 25 à 90)
+data_pred200 = pd.DataFrame({
+    'Temperature': np.linspace(25, 90, 100),
+    'Intercept': 1
+})
+data_pred200['Proportion_Pred'] = logmodelP200.predict(data_pred200[['Intercept', 'Temperature']])
+
+# Tracer les points observés
+plt.figure()  # Nouvelle figure
+plt.scatter(dataP200['Temperature'], dataP200['Proportion'], color='blue', label='Observations')
+
+# Tracer la courbe de régression logistique
+plt.plot(data_pred200['Temperature'], data_pred200['Proportion_Pred'], color='red', label='Régression logistique')
+
+# Ajouter des labels et une légende
+plt.xlabel('Température')
+plt.ylabel('Proportion de dysfonctionnements')
+plt.title('Régression logistique sur les données P200')
+plt.legend()
+plt.grid(True)
+
+# Créer une image regression_logistique_P200.png
+plt.savefig("regression_logistique_P200.png")
+print("regression_logistique_P200.png")
+#+end_src
+
+
+
+*** Régression logistique sur toutes données
+
+#+begin_src python :results value :session second :exports both
+# Filtrer les données et créer une copie
+dataAll = data.copy()
+# Calcul des colonnes nécessaires
+dataAll["Success"] = dataAll["Count"] - dataAll["Malfunction"]
+dataAll["Intercept"] = 1
+# Création du modèle avec une instance de sm.families.links.Logit()
+logmodelAll = sm.GLM(
+    dataAll['Malfunction'],
+    dataAll[['Intercept', 'Temperature']],
+    family=sm.families.Binomial(link=sm.families.links.Logit())
+).fit()
+# Résumé du modèle
+logmodelAll.summary()
+# Générer des données pour la prédiction (de 25 à 90)
+data_predAll = pd.DataFrame({
+    'Temperature': np.linspace(25, 90, 100),
+    'Intercept': 1
+})
+# Ajouter les prédictions au DataFrame
+data_predAll['Malfunction_Prob'] = logmodelAll.predict(data_predAll[['Intercept', 'Temperature']])
+# Tracer les points observés
+plt.figure()  # Nouvelle figure
+plt.scatter(dataAll['Temperature'], dataAll['Malfunction'], color='blue', label='Observations')
+# Tracer la courbe de régression logistique
+plt.plot(data_predAll['Temperature'], data_predAll['Malfunction_Prob'], color='red', label='Régression logistique')
+# Ajouter des labels et une légende
+plt.xlabel('Température')
+plt.ylabel('Probabilité de dysfonctionnement')
+plt.title('Régression logistique sur toutes les données')
+# Ajouter une légende
+plt.legend()
+# Créer une image regression_logistique_All.png
+plt.savefig("regression_logistique_All.png")
+print("regression_logistique_All.png")
+#+end_src
+
+#+RESULTS:
+[[file:regression_logistique_All.png]]
+
+