Missile defense for HUAM
Alexander Levakov, Senior Research Fellow, Ph.D
March, 2019
HUAM missile defense
1 About HUAM missile defence
We use CNN article about HUAM missile defense posture.
For more see:
https://edition.cnn.com/2017/03/07/asia/north-korea-japan-us-ballistic-missile-defense/index.html
2 Missiles Test Events Data
Here we use open source data discussed in the previous topics:
1.Hwasong-12 threat to HUAM - https://rpubs.com/alex-lev/478747
2.Demystifying THAAD Reliability - https://rpubs.com/alex-lev/477708
3.Demystifying AEGIS Reliability - https://rpubs.com/alex-lev/461808
No official information is available about PAC-3 IRBM intercept reliability so far. Here we use some questionable estimate of Patriot PAC-3 as 10 hits out of 11 trials.
3 Bayesian posterior samples for missiles reliability
library(rstan)
library(BH)
library(bayesplot)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
#STAN binomial model
print(stan_model(file = "binom.stan"))## S4 class stanmodel 'binom' coded as follows:
## data{
## int n; // Number of trials
## int x; // Number of success
## int a; // Alpha
## int b; // Beta
## }
##
## parameters{
## real<lower=0,upper=1> p; // unknown parameter
## }
##
## model{
## //p~uniform(0,1);
## p~beta(a+1,b+1);
## x~binomial(n,p);
## //x ~ beta_binomial(n, a+1, b+1);
##
##
## }
##
## generated quantities {
## real pp;
## pp = beta_binomial_rng(n, a+1, b+1);
## }
##
## #KN17 N=6,X=3
data_list_kn17 <- list(n=6,x=3,a=2,b=2)
fit.kn17 <- stan(file = "binom.stan",data = data_list_kn17,chains = 2,
iter = 5000,warmup = 1000,seed = 12345)
print(fit.kn17)## Inference for Stan model: binom.
## 2 chains, each with iter=5000; warmup=1000; thin=1;
## post-warmup draws per chain=4000, total post-warmup draws=8000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## p 0.50 0.00 0.14 0.24 0.41 0.50 0.60 0.76 3005 1
## pp 2.96 0.02 1.60 0.00 2.00 3.00 4.00 6.00 7997 1
## lp__ -8.81 0.01 0.71 -10.83 -8.98 -8.54 -8.37 -8.32 3556 1
##
## Samples were drawn using NUTS(diag_e) at Thu Mar 28 13:12:07 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).stan_hist(fit.kn17)
s1.kn17 <- extract(fit.kn17)
#AEGIS N=52,X=42
data_list_aegis <- list(n=52,x=42,a=2,b=2)
fit.aegis <- stan(file = "binom.stan",data = data_list_aegis,chains = 2,
iter = 5000,warmup = 1000,seed = 12345)
print(fit.aegis)## Inference for Stan model: binom.
## 2 chains, each with iter=5000; warmup=1000; thin=1;
## post-warmup draws per chain=4000, total post-warmup draws=8000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## p 0.78 0.00 0.05 0.66 0.74 0.78 0.81 0.87 3159 1
## pp 25.99 0.12 10.32 7.00 18.00 26.00 34.00 45.00 8000 1
## lp__ -31.37 0.01 0.70 -33.40 -31.54 -31.09 -30.92 -30.86 2885 1
##
## Samples were drawn using NUTS(diag_e) at Thu Mar 28 13:12:11 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).stan_hist(fit.aegis)
s1.aegis <- extract(fit.aegis)
#THAAD N=11,X=5
data_list_thaad <- list(n=11,x=5,a=2,b=2)
fit.thaad <- stan(file = "binom.stan",data = data_list_thaad,chains = 2,
iter = 5000,warmup = 1000,seed = 12345)
print(fit.thaad)## Inference for Stan model: binom.
## 2 chains, each with iter=5000; warmup=1000; thin=1;
## post-warmup draws per chain=4000, total post-warmup draws=8000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## p 0.47 0.00 0.12 0.25 0.39 0.47 0.55 0.70 2762 1
## pp 5.50 0.03 2.57 1.00 4.00 5.00 7.00 10.00 7978 1
## lp__ -12.26 0.01 0.69 -14.23 -12.43 -11.99 -11.81 -11.75 2953 1
##
## Samples were drawn using NUTS(diag_e) at Thu Mar 28 13:12:15 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).stan_hist(fit.thaad)
s1.thaad <- extract(fit.thaad)
#PATRIOT PAC-3 N=11,X=10
data_list_pac3 <- list(n=11,x=10,a=2,b=2)
fit.pac3 <- stan(file = "binom.stan",data = data_list_pac3,chains = 2,
iter = 5000,warmup = 1000,seed = 12345)
print(fit.pac3)## Inference for Stan model: binom.
## 2 chains, each with iter=5000; warmup=1000; thin=1;
## post-warmup draws per chain=4000, total post-warmup draws=8000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## p 0.77 0.00 0.10 0.55 0.70 0.78 0.84 0.93 2823 1
## pp 5.50 0.03 2.58 1.00 4.00 6.00 7.00 10.00 8000 1
## lp__ -9.79 0.01 0.72 -11.85 -9.97 -9.52 -9.33 -9.28 3068 1
##
## Samples were drawn using NUTS(diag_e) at Thu Mar 28 13:12:18 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).stan_hist(fit.pac3)
s1.pac3 <- extract(fit.pac3)
kn17 <- s1.kn17$p
aegis <- s1.aegis$p
thaad <- s1.thaad$p
pac3 <- s1.pac3$p4 Reliability of missile defense
4.1 One AEGIS, one THAAD and one PAC-3 missiles
library(ggplot2)
prob_hist <- function(prob,num1,num2,num3){
hist_label <- "Probability density of KN17 intercept by"
ggplot(data=as.data.frame(prob), aes(prob))+
geom_histogram(bins = 30,col="black",fill="green") +
geom_vline(xintercept = mean(prob), color = "red") +
geom_errorbarh(aes(y=-5, xmin=quantile(prob,0.1),
xmax=quantile(prob,0.9)),
data=as.data.frame(prob), col="#0094EA", size=3) +
ggtitle(paste(hist_label,num1," AEGIS,",num2," THAAD," ,num3," PAC-3",sep=" ")) + labs(x="Probability")
}
bmd_1_1 <- (1-kn17*(1-aegis)*(1-thaad)*(1-pac3))
quantile(bmd_1_1,c(0.1,0.5,0.9))## 10% 50% 90%
## 0.9742350 0.9881510 0.9952235prob_hist(bmd_1_1,"1","1","1")
4.2 Two AEGIS, one THAAD and one PAC-3 missiles
bmd_2_1 <- (1-kn17*(1-aegis)^2*(1-thaad)*(1-pac3))
quantile(bmd_2_1,c(0.1,0.5,0.9))## 10% 50% 90%
## 0.9934273 0.9974410 0.9991137prob_hist(bmd_2_1,"2","1","1")
4.3 Two AEGIS, two THAAD and two PAC-3 missiles
bmd_2_2 <- (1-kn17*(1-aegis)^2*(1-thaad)^2*(1-pac3)^2)
quantile(bmd_2_2,c(0.1,0.5,0.9))## 10% 50% 90%
## 0.9988094 0.9997078 0.9999455prob_hist(bmd_2_2,"2","2","2")
4.4 Four AEGIS, two THAAD and two PAC-3 missiles
bmd_4_2 <- (1-kn17*(1-aegis)^4*(1-thaad)^2*(1-pac3)^2)
quantile(bmd_4_2,c(0.1,0.5,0.9))## 10% 50% 90%
## 0.9999209 0.9999865 0.9999982prob_hist(bmd_4_2,"4","2","2")
5 Conclusion
- Probability for the KN17 intercept by single AEGIS, single THAAD and single PAC-3 missiles provided KN17 would have success launch event can be estimated as 95% credible interval \(0.974\le P_{1,1,1}\le0.995\) with the mean \(M[P_{1,1,1}]=0.988\).
- Probability for the KN17 intercept by two AEGIS, single THAAD and single PAC-3 missiles provided KN17 would have success launch event can be estimated as 95% credible interval \(0.993\le P_{2,1,1}\le0.999\) with the mean \(M[P_{2,1,1}]=0.997\).
- Probability for the KN17 intercept by two AEGIS, two THAAD and two PAC-3 missiles provided KN17 would have success launch event can be estimated as 95% credible interval \(0.9988\le P_{2,2,2}\le0.9999\) with the mean \(M[P_{2,2,2}]=0.9997\).
- Probability for the KN17 intercept by four AEGIS, two THAAD and two PAC-3 missiles provided KN17 would have success launch event can be estimated as 95% credible interval \(0.99992\le P_{4,2,2}\le0.99999\) with the mean \(M[P_{4,2,2}]=0.99998\).
- All the estimates made above are based on the assumption that KN17, AEGIS, THAAD and PAC-3 flight test data are correct in terms of \(Success\) and \(Failure\) events.