Commit 440bda27 authored by Dorian Goepp's avatar Dorian Goepp

Exercises 2, 3 and 5

parent 7808267b
#+TITLE: Votre titre #+TITLE: Analyse de données simples
#+AUTHOR: Votre nom #+AUTHOR: Dorian Goepp
#+DATE: La date du jour #+DATE: 2022-10-25
#+LANGUAGE: fr #+LANGUAGE: fr
# #+PROPERTY: header-args :eval never-export # #+PROPERTY: header-args :eval never-export
...@@ -82,3 +82,18 @@ faisant ~<r~ ou ~<R~ suivi de ~Tab~). ...@@ -82,3 +82,18 @@ faisant ~<r~ ou ~<R~ suivi de ~Tab~).
Maintenant, à vous de jouer! Vous pouvez effacer toutes ces Maintenant, à vous de jouer! Vous pouvez effacer toutes ces
informations et les remplacer par votre document computationnel. informations et les remplacer par votre document computationnel.
* Mon travail
#+begin_src R :results output :session *R* :exports both
donnees = c(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)
summary(donnees)
sd(donnees)
#+end_src
#+RESULTS:
:
: Min. 1st Qu. Median Mean 3rd Qu. Max.
: 2.80 11.85 14.50 14.11 16.80 23.40
:
: [1] 4.334094
...@@ -82,3 +82,22 @@ faisant ~<r~ ou ~<R~ suivi de ~Tab~). ...@@ -82,3 +82,22 @@ faisant ~<r~ ou ~<R~ suivi de ~Tab~).
Maintenant, à vous de jouer! Vous pouvez effacer toutes ces Maintenant, à vous de jouer! Vous pouvez effacer toutes ces
informations et les remplacer par votre document computationnel. informations et les remplacer par votre document computationnel.
* Mon analyse
#+begin_src R :results output :session *R* :exports both
donnees = c(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_src R :results output graphics :file (org-babel-temp-file "figure" ".png") :exports both :width 600 :height 400 :session *R*
plot(donnees, type="l")
#+end_src
#+RESULTS:
#+begin_src R :results output graphics :file (org-babel-temp-file "figure" ".png") :exports both :width 600 :height 400 :session *R*
hist(donnees)
#+end_src
#+RESULTS:
...@@ -38,34 +38,36 @@ Nous commençons donc par charger ces données: ...@@ -38,34 +38,36 @@ Nous commençons donc par charger ces données:
#+begin_src R :results output :session *R* :exports both #+begin_src R :results output :session *R* :exports both
data = read.csv("shuttle.csv",header=T) data = read.csv("shuttle.csv",header=T)
data data
original_data = data
#+end_src #+end_src
#+RESULTS: #+RESULTS:
#+begin_example #+begin_example
Date Count Temperature Pressure Malfunction
1 4/12/81 6 66 50 0 Date Count Temperature Pressure Malfunction
2 11/12/81 6 70 50 1 1 4/12/81 6 66 50 0
3 3/22/82 6 69 50 0 2 11/12/81 6 70 50 1
4 11/11/82 6 68 50 0 3 3/22/82 6 69 50 0
5 4/04/83 6 67 50 0 4 11/11/82 6 68 50 0
6 6/18/82 6 72 50 0 5 4/04/83 6 67 50 0
7 8/30/83 6 73 100 0 6 6/18/82 6 72 50 0
8 11/28/83 6 70 100 0 7 8/30/83 6 73 100 0
9 2/03/84 6 57 200 1 8 11/28/83 6 70 100 0
10 4/06/84 6 63 200 1 9 2/03/84 6 57 200 1
11 8/30/84 6 70 200 1 10 4/06/84 6 63 200 1
12 10/05/84 6 78 200 0 11 8/30/84 6 70 200 1
13 11/08/84 6 67 200 0 12 10/05/84 6 78 200 0
14 1/24/85 6 53 200 2 13 11/08/84 6 67 200 0
15 4/12/85 6 67 200 0 14 1/24/85 6 53 200 2
16 4/29/85 6 75 200 0 15 4/12/85 6 67 200 0
17 6/17/85 6 70 200 0 16 4/29/85 6 75 200 0
18 7/29/85 6 81 200 0 17 6/17/85 6 70 200 0
19 8/27/85 6 76 200 0 18 7/29/85 6 81 200 0
20 10/03/85 6 79 200 0 19 8/27/85 6 76 200 0
21 10/30/85 6 75 200 2 20 10/03/85 6 79 200 0
22 11/26/85 6 76 200 0 21 10/30/85 6 75 200 2
23 1/12/86 6 58 200 1 22 11/26/85 6 76 200 0
23 1/12/86 6 58 200 1
#+end_example #+end_example
Le jeu de données nous indique la date de l'essai, le nombre de joints Le jeu de données nous indique la date de l'essai, le nombre de joints
...@@ -85,6 +87,7 @@ data ...@@ -85,6 +87,7 @@ data
#+end_src #+end_src
#+RESULTS: #+RESULTS:
:
: Date Count Temperature Pressure Malfunction : Date Count Temperature Pressure Malfunction
: 2 11/12/81 6 70 50 1 : 2 11/12/81 6 70 50 1
: 9 2/03/84 6 57 200 1 : 9 2/03/84 6 57 200 1
...@@ -104,12 +107,37 @@ plot(data=data, Malfunction/Count ~ Temperature, ylim=c(0,1)) ...@@ -104,12 +107,37 @@ plot(data=data, Malfunction/Count ~ Temperature, ylim=c(0,1))
#+end_src #+end_src
#+RESULTS: #+RESULTS:
[[file:freq_temp.png]]
À première vue, ce n'est pas flagrant mais bon, essayons quand même À 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 d'estimer l'impact de la température $t$ sur la probabilité de
dysfonctionnements d'un joint. dysfonctionnements d'un joint.
** Admettons qu'en fait les tests sans défaut de joint soient importants
Comment les représenterait-on alors ?
#+begin_src R :results output graphics :file (org-babel-temp-file "figure" ".png") :exports both :width 600 :height 400 :session *R*
plot(data=original_data, Malfunction/Count ~ Temperature, ylim=c(0, 1))
#+end_src
#+RESULTS:
Toutes les observations en dessous de 66°F ont rencontré au moins une
défaillance de joint torique. Ce n'est probablement pas une
information à négliger !
** Et si la pression avait un rôle ?
#+begin_src R :results output graphics :file (org-babel-temp-file "figure" ".png") :exports both :width 600 :height 400 :session *R*
plot(data=original_data, Malfunction/Count ~ Pressure, ylim=c(0, 1))
#+end_src
#+RESULTS:
Une fois encore, nous avons très peu de données, mais il ne semble
rien y avoir à en tirer.
* Estimation de l'influence de la température * Estimation de l'influence de la température
Supposons que chacun des 6 joints toriques est endommagé avec la même Supposons que chacun des 6 joints toriques est endommagé avec la même
...@@ -120,6 +148,10 @@ température $t$ suit une loi binomiale de paramètre $n=6$ et ...@@ -120,6 +148,10 @@ 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 $p=p(t)$. Pour relier $p(t)$ à $t$, on va donc effectuer une
régression logistique. régression logistique.
> Les joints peuvent défaillir en même temps s'ils sont le produit d'une
même série de production, ce qui limite le réalisme de l'hypothèse
d'indépendance des probabilités de défaut d'un joint.
#+begin_src R :results output :session *R* :exports both #+begin_src R :results output :session *R* :exports both
logistic_reg = glm(data=data, Malfunction/Count ~ Temperature, weights=Count, logistic_reg = glm(data=data, Malfunction/Count ~ Temperature, weights=Count,
family=binomial(link='logit')) family=binomial(link='logit'))
...@@ -156,6 +188,50 @@ et l'erreur standard de cet estimateur est de 0.049, autrement dit on ...@@ -156,6 +188,50 @@ et l'erreur standard de cet estimateur est de 0.049, autrement dit on
ne peut pas distinguer d'impact particulier et il faut prendre nos ne peut pas distinguer d'impact particulier et il faut prendre nos
estimations avec des pincettes. estimations avec des pincettes.
** Refaisons l'analyse avec toutes les données
De même que dans la section précédente, prenons en compte les
événements sans défaut de joint.
#+begin_src R :results output :session *R* :exports both
logistic_reg_full = glm(data=original_data, Malfunction/Count ~ Temperature, weights=Count,
family=binomial(link='logit'))
summary(logistic_reg)
#+end_src
#+RESULTS:
#+begin_example
Call:
glm(formula = Malfunction/Count ~ Temperature, family = binomial(link = "logit"),
data = data, weights = Count)
Deviance Residuals:
2 9 10 11 14 21 23
-0.3015 -0.2836 -0.2919 -0.3015 0.6891 0.6560 -0.2850
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.389528 3.195752 -0.435 0.664
Temperature 0.001416 0.049773 0.028 0.977
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1.3347 on 6 degrees of freedom
Residual deviance: 1.3339 on 5 degrees of freedom
AIC: 18.894
Number of Fisher Scoring iterations: 4
#+end_example
Cette fois-ci, avec la même erreur standard, on obtient un coefficient
bien supérieur, donc une corrélation bien plus forte entre la
température et les défaillances.
Je ne suis pas cappable d'évaluer la signification de ces résultats,
n'ayant pas eu le temps d'étudier la régression logistique. J'admets
donc, pour pouvoir continuer, que l'incertitude n'est pas trop élevée.
* Estimation de la probabilité de dysfonctionnant des joints toriques * Estimation de la probabilité de dysfonctionnant des joints toriques
La température prévue le jour du décollage est de 31°F. Essayons La température prévue le jour du décollage est de 31°F. Essayons
d'estimer la probabilité de dysfonctionnement des joints toriques à d'estimer la probabilité de dysfonctionnement des joints toriques à
...@@ -185,6 +261,7 @@ sum(data_full$Malfunction)/sum(data_full$Count) ...@@ -185,6 +261,7 @@ sum(data_full$Malfunction)/sum(data_full$Count)
#+end_src #+end_src
#+RESULTS: #+RESULTS:
:
: [1] 0.06521739 : [1] 0.06521739
Cette probabilité est donc d'environ $p=0.065$, sachant qu'il existe Cette probabilité est donc d'environ $p=0.065$, sachant qu'il existe
...@@ -205,3 +282,27 @@ problème... Saurez-vous le trouver ? Vous êtes libre de modifier cette ...@@ -205,3 +282,27 @@ 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 analyse et de regarder ce jeu de données sous tous les angles afin
d'expliquer ce qui ne va pas. d'expliquer ce qui ne va pas.
** En prenant en compte les situations sans défaillance
#+begin_src R :results output graphics :file "proba_estimate_all.png" :exports both :width 600 :height 400 :session *R*
# shuttle=shuttle[shuttle$r!=0,]
tempv = seq(from=30, to=90, by = .5)
rmv <- predict(logistic_reg_full, list(Temperature=tempv), type="response")
plot(tempv, rmv, type="l", ylim=c(0, 1))
points(data=data, Malfunction/Count ~ Temperature)
#+end_src
#+RESULTS:
Et maintenant la probabilité de défaillance d'un joint:
#+begin_src R :results output :session *R* :exports both
rmv[1]
#+end_src
#+RESULTS:
: 1
: 0.834373
The probability is here far too high to consider taking off, admitting
that the pressure for the takeoff day was 200 PSI.
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