DF41 range and reliability
Alexander Levakov, Senior Research Fellow, Ph.D
January, 2016
Research goal
We want to estimate potential DF-41 missile threat in terms of Range (km) and Reliability (%). For this purpose we use linear regression model (https://rpubs.com/alex-lev/209272) and Bayesian binomial test (https://rpubs.com/alex-lev/243536).
DF-41
Dongfeng-41
For more about DF-41 see: https://en.wikipedia.org/wiki/DF-41, https://missilethreat.csis.org/missile/df-41/, http://www.globalsecurity.org/jhtml/jframe.html#http://www.globalsecurity.org/wmd/world/china/images/df-41-image04.jpg.
Data
So we have got some open source information about DF-41: S - number of stages (3), D - diameter (2.25 m) and M - weight (80 ton).
Can we estimate potential range of DF-41 using our previous linear regression model based on SLBM data (https://rpubs.com/alex-lev/231516)?
load(file = "slbm.dat")
attach(slbm)
slbm## S M L D R P W
## 1 1 5400 10.40 0.58 150 975 1
## 2 1 13700 11.80 1.30 650 1597 1
## 3 1 19650 14.20 1.30 1420 1179 1
## 4 1 14200 8.89 1.50 2400 650 1
## 5 1 14200 8.89 1.50 3000 650 1
## 6 1 14200 9.65 1.50 3000 650 2
## 7 2 33300 13.00 1.80 7800 1100 1
## 8 2 33300 13.00 1.80 9100 1100 1
## 9 2 26900 10.60 1.54 3900 450 1
## 10 2 35300 14.10 1.80 6500 1600 2
## 11 2 35300 14.10 1.80 8000 1600 1
## 12 2 35300 14.10 1.80 6500 1600 2
## 13 3 90000 16.00 2.40 8300 2550 2
## 14 3 40300 14.80 1.90 8300 2300 2
## 15 3 40800 14.80 1.90 8300 2800 2
## 16 3 36800 11.50 2.00 9300 1150 2
## 17 2 13600 9.45 1.37 2800 500 1
## 18 2 16200 9.86 1.37 4630 760 2
## 19 2 29485 10.36 1.88 5600 2000 2
## 20 3 32000 10.36 1.88 7400 1360 2
## 21 3 57500 13.42 2.11 11000 2880 2
## 22 2 20000 10.70 1.50 3000 1360 1
## 23 2 19500 10.70 1.50 3200 1360 1
## 24 2 19950 10.40 1.50 3200 1000 1
## 25 2 14700 10.70 1.40 2500 600 1
## 26 2 23000 13.00 2.00 8000 700 2Linear regression model
Here is the model for \(R\) - range of missile: \(R=f(S,D,M) + E\), where \(E=N(0,{\sigma})\) - error of estimate. \({R} = e^{{3.442} + {0.322}{L} + {2.984}{D} - {0.00002517}{M}}\)
Frequentis case
fit.slbm<-lm(log(R)~S + D + M, data = slbm) #linear regression model
summary(fit.slbm) # summary of the model##
## Call:
## lm(formula = log(R) ~ S + D + M, data = slbm)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.82063 -0.23408 0.00621 0.18996 0.67479
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.442e+00 4.336e-01 7.938 6.72e-08 ***
## S 3.220e-01 1.547e-01 2.082 0.04922 *
## D 2.984e+00 3.766e-01 7.922 6.94e-08 ***
## M -2.517e-05 7.535e-06 -3.340 0.00296 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3461 on 22 degrees of freedom
## Multiple R-squared: 0.8829, Adjusted R-squared: 0.867
## F-statistic: 55.32 on 3 and 22 DF, p-value: 2.066e-10shapiro.test(fit.slbm$residuals) #test for normality of the residuals##
## Shapiro-Wilk normality test
##
## data: fit.slbm$residuals
## W = 0.98117, p-value = 0.8978#95% confident interval prediction for DF-41 Range (km)
exp(predict(fit.slbm,data.frame(S=3,D=2.25,M=80000),interval = "conf",level = 0.95))## fit lwr upr
## 1 9021.785 5450.503 14933.05Not bad! The true value for \(R\) is lying in the 95% confident interval \({5450.503}\leqslant{R}\leqslant{14933.05}\) with the mean value \(R=9021.785\) that is very close to the upper bound of DF-41 - \(R=15000\).
Bayesian case
library(bayesboot)
set.seed(12345)
lm.coefs <- function(d, w) {
coef( lm(log(R)~S + D + M, data = d, weights = w) )
}
blm <- bayesboot(slbm, lm.coefs, R = 10000, use.weights = TRUE)
plot(blm)
summary(blm)## Bayesian bootstrap
##
## Number of posterior draws: 10000
##
## Summary of the posterior (with 95% Highest Density Intervals):
## statistic mean sd hdi.low hdi.high
## (Intercept) 3.567428e+00 5.150899e-01 2.703757e+00 4.665412e+00
## S 3.123405e-01 1.323728e-01 5.670599e-02 5.774338e-01
## D 2.880673e+00 4.071315e-01 2.012368e+00 3.542198e+00
## M -2.224895e-05 7.933479e-06 -3.617302e-05 -5.798076e-06
##
## Quantiles:
## statistic q2.5% q25% median q75%
## (Intercept) 2.775375e+00 3.196509e+00 3.479382e+00 3.873703e+00
## S 7.635106e-02 2.237353e-01 3.029292e-01 3.905000e-01
## D 2.002252e+00 2.617766e+00 2.931485e+00 3.177034e+00
## M -3.523914e-05 -2.793057e-05 -2.321088e-05 -1.758096e-05
## q97.5%
## 4.754934e+00
## 6.042734e-01
## 3.539850e+00
## -4.053415e-06
##
## Call:
## bayesboot(data = slbm, statistic = lm.coefs, R = 10000, use.weights = TRUE)B0<-c()
B1<-c()
B2<-c()
B3<-c()
R<-c()
for (i in 1:10000) {
B0[i]<-blm$`(Intercept)`[i]
B1[i]<-blm$S[i]
B2[i]<-blm$D[i]
B3[i]<-blm$M[i]
R[i]<-exp(B0[i] + B1[i]*3 + B2[i]*2.25 + B3[i]*80000)
}
quantile(R,c(0.025,0.975))## 2.5% 97.5%
## 7476.808 18413.735mean(R)## [1] 10280.93median(R)## [1] 9358.418mean(R>15000)## [1] 0.0702hist(R,breaks = 50,freq = F,col="lightblue",xlim = c(5000,20000))
curve(dnorm(x,mean = mean(R),sd = sd(R)),add = T)
abline(v = mean(R),col="blue")
abline(v = median(R),col="red")
Binomial test
DF-41 has got 7 successful test flights so far. Does it mean 100% reliability? Let’s see!
Frequentist case
What is the 95% confidence interval for DF-41 reliability if our null hypothesis \(H_0|P=0.95\) is true?
binom.test(7,7,0.95,alternative = "greater")##
## Exact binomial test
##
## data: 7 and 7
## number of successes = 7, number of trials = 7, p-value = 0.6983
## alternative hypothesis: true probability of success is greater than 0.95
## 95 percent confidence interval:
## 0.6518363 1.0000000
## sample estimates:
## probability of success
## 1Well, not so bad! As frequentists we can be 95% confident that \(0.65<P_{DF41}{<}1\). What about Bayes?
Bayesian case
We try Bayesian binomial test for our data converting confidence into credibility.
library(rjags)## Loading required package: coda## Linked to JAGS 4.1.0## Loaded modules: basemod,bugslibrary(stringr)
source("lib/generic_functions.R")
source("lib/utility_functions.R")
source("lib/BayesianFirstAid-package.R")
source("lib/bayes_binom_test.R")
source("lib/bayes_cor_test.R")
source("lib/bayes_poisson_test.R")
source("lib/bayes_prop_test.R")
source("lib/bayes_t_test.R")
fit.bs <- bayes.binom.test(x = 7, n = 7, cred.mass=0.95, p = 0.95, n.iter = 5000)
summary(fit.bs)## Data
## number of successes = 7, number of trials = 7
##
## Model parameters and generated quantities
## theta: the relative frequency of success
## x_pred: predicted number of successes in a replication
##
## Measures
## mean sd HDIlo HDIup %<comp %>comp
## theta 0.889 0.098 0.69 1 0.666 0.334
## x_pred 6.230 1.025 4.00 7 0.000 1.000
##
## 'HDIlo' and 'HDIup' are the limits of a 95% HDI credible interval.
## '%<comp' and '%>comp' are the probabilities of the respective parameter being
## smaller or larger than 0.95.
##
## Quantiles
## q2.5% q25% median q75% q97.5%
## theta 0.637 0.842 0.917 0.964 0.997
## x_pred 4.000 6.000 7.000 7.000 7.000plot(fit.bs)
Bayes can’t give us 100% reliability for DF-41 launch so far. In fact we have got 95% credible interval \(0.684 {\leq} P_{DF41} {\leq}1)\) where \(P(P_{DF41}{>}0.95)=0.338\) and \(P(P_{DF41}{<}0.95)=0.662\) with the mean value \(P_{DF41}=0.887\).
Conclusions
- DF-41 potential range can be estimated as \({5450.503}\leqslant{R}\leqslant{14933.05}\) with the mean \(9021.785\) when applying frequentist point of view.
- Bayes gives 95% credible interval for DF-41 range as \(7470 {\leq} R_{DF41} {\leq} 18174\) with the mean \(10260\) and median \(9349\) while \(P(R>15000)=0.067\) that is 15000 km is too much for DF-41.
- DF-41 reliability for the successful launch is less than 95% with the probability 66% and greater 95% with the probability 34% with the mean value \(P_{DF41}=0.887\) so far.
- I believe Bayes!