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

* Chargement des données
  Nous commençons donc par charger ces données:
  #+begin_src python :results value :session *python* :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 lanceur principal), la
  température (en Fahrenheit) et la pression (en psi), et enfin le
  nombre de dysfonctionnements relevés. 

* Inspection graphique des données

  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 *python* 
    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]]

  Comment la fréquence d'échecs varie-t-elle avec la pression ?
  #+begin_src python :results output file :var matplot_lib_filename="freq_pressure_python.png" :exports both :session *python* 
    import matplotlib.pyplot as plt

    plt.clf()
    data["Frequency"]=data.Malfunction/data.Count
    data.plot(x="Pressure",y="Frequency",kind="scatter",ylim=[0,1])
    plt.grid(True)

    plt.savefig(matplot_lib_filename)
    print(matplot_lib_filename)
  #+end_src

  #+RESULTS:
  [[file:freq_pressure_python.png]]

  Avec si peu de points de pression, on ne peut pas étudier grand chose. De
  plus, visuellement, aucune relation n'a l'air de sauter aux yeux.

  À 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 *python* :exports both
    import statsmodels.api as sm

    data["Intercept"]=1

    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:                   23
  Model:                            GLM   Df Residuals:                       21
  Model Family:                Binomial   Df Model:                            1
  Link Function:                  logit   Scale:                          1.0000
  Method:                          IRLS   Log-Likelihood:                -3.9210
  Date:                Mon, 20 Jul 2020   Deviance:                       3.0144
  Time:                        11:22:26   Pearson chi2:                     5.00
  No. Iterations:                     6                                         
  Covariance Type:            nonrobust                                         
  ===============================================================================
                    coef    std err          z      P>|z|      [0.025      0.975]
  -------------------------------------------------------------------------------
  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.

* 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 *python* 
    import matplotlib.pyplot as plt

    data_pred = pd.DataFrame({'Temperature': np.linspace(start=30, stop=90, num=121), 'Intercept': 1})
    prediction = logmodel.get_prediction(data_pred[ ['Intercept','Temperature'] ]).summary_frame()
    data_pred['Frequency'] = prediction["mean"]
    data_pred.plot(x="Temperature",y="Frequency",kind="line",ylim=[0,1])
    plt.scatter(x=data["Temperature"], y=data["Frequency"])
    plt.fill_between(
        data_pred["Temperature"],
        prediction["mean_ci_lower"],
        prediction["mean_ci_upper"],
        alpha=0.2,
    )
    plt.grid(True)

    plt.savefig(matplot_lib_filename)
    print(matplot_lib_filename)
  #+end_src

  #+RESULTS:
  [[file:proba_estimate_python.png]]

  Hmm, l'intervalle de confiance semble très grand. Peut-on déduire quoi que ce
  soit de ces données ?