Bayes and NASA
Alexander Levakov
March, 2015
library(bayesboot)
library(MCMCpack)## Loading required package: coda## Loading required package: MASS## ##
## ## Markov Chain Monte Carlo Package (MCMCpack)## ## Copyright (C) 2003-2016 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park## ##
## ## Support provided by the U.S. National Science Foundation## ## (Grants SES-0350646 and SES-0350613)
## ##library(ggplot2)Space Shuttle Challenger disaster
On Jan. 28, 1986, the Space Shuttle Challenger broke apart 73 seconds into its 10th mission. Forecasts for January 28 predicted an unusually cold morning, with temperatures close to −1 °C (30 °F), the minimum temperature permitted for launch. The Shuttle was never certified to operate in temperatures that low.
See: https://en.wikipedia.org/wiki/Space_Shuttle_Challenger_disaster
Inquiry of the accident
The Presidential Commission on the Space Shuttle Challenger Accident, also known as the Rogers Commission (after its chairman), was formed to investigate the disaster. In view of the findings, the Commission concluded that the cause of the Challenger accident was the failure of the pressure seal in the aft field joint of the right Solid Rocket Motor. The failure was due to a faulty design unacceptably sensitive to a number of factors. These factors were the effects of temperature, physical dimensions, the character of materials, the effects of reusability, processing, and the reaction of the joint to dynamic loading.
Data
Number of previous flights - 23.
Distress of booster joint
Yes=1 , No=0
y = c(1,1,1,1,0,0,0,0,0,0,0,0,1,1,0,0,0,1,0,0,0,0,0)
table(y)## y
## 0 1
## 16 7Temperature at launch time in F
x = c(53,57,58, 63,66, 67,67, 67,68, 69,70,70,70, 70,72, 73,75, 75,76, 76,78, 79,81)
summary(x)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 53.00 67.00 70.00 69.57 75.00 81.00hist(x,breaks = 3)
qqnorm(x)
shapiro.test(x)##
## Shapiro-Wilk normality test
##
## data: x
## W = 0.95368, p-value = 0.3482As we can see temperature at launch time for 23 flights is normally distributed.
Preliminary estimate
The data for Shuttle flights is a real experiment with two possible outcomes for booster joint status (1 - distress, 0 - normal). So we got our first estimate of the probability of normal booster joint status during lift off: \(P=16/23=0.695\).
plot(x,y,col=factor(y))
abline(v = 65)
grid()
As we can see temperature at 65 F point and above is more convenient and secure for shuttle launch provided the data we use for our estimates.
Logistic regression model
https://en.wikipedia.org/wiki/Logistic_regression
fit.sh <- glm(y ~ x, family=binomial)
summary(fit.sh)##
## Call:
## glm(formula = y ~ x, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0611 -0.7613 -0.3783 0.4524 2.2175
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 15.0429 7.3786 2.039 0.0415 *
## x -0.2322 0.1082 -2.145 0.0320 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 28.267 on 22 degrees of freedom
## Residual deviance: 20.315 on 21 degrees of freedom
## AIC: 24.315
##
## Number of Fisher Scoring iterations: 5confint(fit.sh)## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) 3.3305848 34.34215133
## x -0.5154718 -0.06082076Logistic model fits the data with significant p-values for B0 and B1 coefficients.
shut<-data.frame(x,y)
binomial_smooth <- function(...) {
geom_smooth(method = "glm", method.args = list(family = "binomial"), ...)
}
ggplot(shut,aes(x,y)) + geom_point(aes(colour=y)) +binomial_smooth()
The plot above demonstrates some possible odds for normal booster joint status during lift off provided temperature is under 65 F.
Bayesian logistic model
There is another way to consider this problem - Bayesian statistics (http://en.wikipedia.org/wiki/Bayesian_statistics). We can use bayesboot functions for this purpose:
glm.coefs <- function(d, w) {
coef( glm(y ~ x, data = shut, family=binomial, weights = w) )
}
set.seed(12345)
fit.bs <- bayesboot(shut, glm.coefs, R = 10000, use.weights = TRUE)
summary(fit.bs)## Bayesian bootstrap
##
## Number of posterior draws: 10000
##
## Summary of the posterior (with 95% Highest Density Intervals):
## statistic mean sd hdi.low hdi.high
## (Intercept) 17.4067694 7.697796 3.4024050 31.37601611
## x -0.2687944 0.118311 -0.4865798 -0.05909408
##
## Quantiles:
## statistic q2.5% q25% median q75% q97.5%
## (Intercept) 5.5968009 12.5346746 16.4640786 21.0286278 35.03862164
## x -0.5422121 -0.3228707 -0.2533843 -0.1934664 -0.09153434
##
## Call:
## bayesboot(data = shut, statistic = glm.coefs, R = 10000, use.weights = TRUE)plot(fit.bs)
posterior <- MCMClogit(y~x, data=shut,mcmc = 10000,seed = 12345,marginal.likelihood = "Laplace")
plot(posterior)
summary(posterior)##
## Iterations = 1001:11000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## (Intercept) 19.1011 8.994 0.08994 0.315128
## x -0.2927 0.132 0.00132 0.004629
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## (Intercept) 4.5819 12.669 18.2416 24.2756 40.07971
## x -0.6001 -0.369 -0.2805 -0.1974 -0.07953Comparing Classical and byasian models results
prob.cl<-1/(1+exp(-(3.33-0.51*30))) # q2.5% for clfssical model
prob.bs<-1/(1+exp(-(4.58-0.6*30))) # q2.5% for byasian model
prob.cl## [1] 6.331291e-06prob.bs## [1] 1.48514e-06Both models give 2.5% chance for normal booster joint status during lift off provided temperature is around 30 F.
Conclusion
The models used for analysis are limited with one factor (temperature F) and can be applied for demonstration purposes only.
Unfortunately NASA had neglected logistic regression (classical and byasian), making decision for Challenger launch with only 2.5% probability for success.