#+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]]