%% Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Faillure
%
% In this document we reperform some of the analysis provided in
% Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure by Siddhartha R. Dalal, Edward B. Fowlkes, Bruce Hoadley published in Journal of the American Statistical Association, Vol. 84, No. 408 (Dec., 1989), pp. 945-957 and available at http://www.jstor.org/stable/2290069.
% On the fourth page of this article, they indicate that the maximum likelihood estimates of the logistic regression using only temperature are: $\hat{\alpha}=5.085$ and $\hat{\beta}=-0.1156$ and $\hat{\beta}=-0.1156$ their asymptotic standard errors are $s_{\hat{\alpha}}=3.052$ and $s_{\hat{\beta}}=0.047$. The Goodness of fit indicated for this model was $G^2=18.086$ with 21 degrees of freedom. Our goal is to reproduce the computation behind these values and the Figure 4 of this article, possibly in a nicer looking way.
%
%% Technical information on the computer on which the analysis is run
% We will use Matlab.

% Peut afficher des informations personnelles
%system('systeminfo')
version

%
%% Loading and inspecting data
% Let's start by reading data.

data = readtable("./shuttle.csv")

% We know from our previous experience on this data set that filtering  data is a really bad idea. We will therefore process it as such.

data.Frequency = data.Malfunction./data.Count;
figure()
plot(data.Temperature,data.Frequency,"o")
ylim([0 1])

%
%% Logistic regression
% Let's assume O-rings independently fail with the same probability which  solely depends on temperature. A logistic regression should allow us to  estimate the influence of temperature.

data.Success = data.Count - data.Malfunction;
data.Intercept = ones(size(data.Success))
logmodel = fitglm([data.Intercept,data.Temperature],data.Frequency,'Distribution','binomial','VarNames',{'Intercept','Temperature','Frequency'})

% The maximum likelyhood estimator of the intercept and of Temperature are thus 𝛼̂ =5.0849 and 𝛽̂ =−0.1156. This corresponds to the values from the article of Dalal et al. The standard errors are 𝑠𝛼̂ =7.477 and 𝑠𝛽̂ =0.115, which is different from the 3.052 and 0.04702 reported by Dallal et al. The deviance is 3.01444 with 21 degrees of freedom. I cannot find any value similar to the Goodness of fit (𝐺2=18.086) reported by Dalal et al. There seems to be something wrong. Oh I know, I haven't indicated that  my observations are actually the result of 6 observations for each  rocket launch. Let's indicate these weights (since the weights are  always the same throughout all experiments, it does not change the  estimates of the fit but it does influence the variance estimates).

logmodel = fitglm([data.Intercept,data.Temperature],data.Frequency,'Distribution','binomial','Link','logit','Weights',data.Count,'VarNames',{'Intercept','Temperature','Frequency'})

% Good, now I have recovered the asymptotic standard errors 𝑠𝛼̂ =3.052 and 𝑠𝛽̂ =0.047. The Goodness of fit (Deviance) indicated for this model is 𝐺2=18.086
%  with 21 degrees of freedom (Df Residuals).
% I have therefore managed to fully replicate the results of the Dalal et al. article.
%
%% Predicting faillure probability
% The temperature when launching the shuttle was 31°F. Let's try to  estimate the failure probability for such temperature using our model.:

data_pred = table();
data_pred.Temperature = linspace(30,90,121)';
data_pred.Intercept = ones(size(data_pred.Temperature));
data_pred.Frequency = logmodel.predict(data_pred);
figure()
plot(data_pred.Temperature,data_pred.Frequency)
ylim([0 1])
hold on
scatter(data.Temperature,data.Frequency,"o")
hold off
grid
legend("Frequency")

% This figure is very similar to the Figure 4 of Dalal et al. I have managed to replicate the Figure 4 of the Dalal et al. article.
%
%% Computing and ploting uncertainty

figure()
xlim([30 90])
ylim([0 1])
plotregression(data.Temperature,data.Frequency)

%