Kim Jong-Un inspects missile test

Kim Jong-Un inspects missile test

Preface

The War on the rocks has published very interesting analysis of the US BMD using GBS kill probability for KN-14 1. Meanwhile this paper has some drawbacks. First, KN-14 is postulated as 100% reliable missile. With 2 test events it is not the case. Second, North Korea has not got full scale test bed for long range ICBM test. Accuracy of RV for Korean missiles is unknown even for Kim. Third, atmospheric reentry test at the speed of 5-7 km per second has no records so far. At last, all estimates have been made by authors with point not interval method. From statistics point of view it is the most simple and primitive estimate of reliability for any stochastic system.

Problem

We want to estimate US BMD system reliability for the KN-11 as a real target 2. To this end we apply Bayesian binomial test and Markov chain model 3.

Data

We use open source data 4.

Bayesian binomial model

We try Bayesian binomial test for our data converting confidence into credibility.

library(expm)
## Loading required package: Matrix
## 
## Attaching package: 'expm'
## The following object is masked from 'package:Matrix':
## 
##     expm
library(diagram)
## Loading required package: shape
library(ggplot2)
library(rjags)
## Loading required package: coda
## Linked to JAGS 4.2.0
## Loaded modules: basemod,bugs
library(mcmc)
library(stringr)
library(bayesplot)
## This is bayesplot version 1.2.0
model_string <- "model{

# Likelihood
Y1 ~ dbinom(theta1,n1)
Y2 ~ dbinom(theta2,n2)
Y3 ~ dbinom(theta3,n3)
Y4 ~ dbinom(theta4,n4)
# Prior
theta1 ~ dbeta(a, b)
theta2 ~ dbeta(a, b)
theta3 ~ dbeta(a, b)
theta4 ~ dbeta(a, b)
}"

Y1 <- 4 # KN-11 failed test events
n1 <- 12# KN-11 all test events
Y2 <- 8# GMD failed test events
n2 <- 18# GMD all test events
Y3 <- 5# THAAD failed test events
n3 <- 11# THAAD all test events
Y4 <- 6# AEGIS failed test events
n4 <- 37# AEGIS all test events


a <- 1
b <- 1
model <- jags.model(textConnection(model_string),                    inits=list(.RNG.name="base::Wichmann-Hill",.RNG.seed= 12341),
                    data = list(Y1=Y1,n1=n1,Y2=Y2,
                                n2=n2,Y3=Y3,n3=n3,
                                Y4=Y4,n4=n4,a=a,b=b))
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 4
##    Unobserved stochastic nodes: 4
##    Total graph size: 14
## 
## Initializing model
update(model, 10000, progress.bar="none"); # Burnin for 10000 samples

samp <- coda.samples(model, 
                     variable.names=c("theta1","theta2","theta3","theta4"), 
                     n.iter=20000, progress.bar="none")

summary(samp)
## 
## Iterations = 10001:30000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 20000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##          Mean      SD  Naive SE Time-series SE
## theta1 0.3572 0.12274 0.0008679      0.0008679
## theta2 0.4506 0.10768 0.0007614      0.0007810
## theta3 0.4636 0.13350 0.0009440      0.0009440
## theta4 0.1794 0.06056 0.0004283      0.0004331
## 
## 2. Quantiles for each variable:
## 
##           2.5%    25%    50%    75%  97.5%
## theta1 0.13980 0.2676 0.3507 0.4396 0.6086
## theta2 0.24542 0.3760 0.4487 0.5236 0.6660
## theta3 0.21107 0.3688 0.4613 0.5569 0.7270
## theta4 0.07767 0.1352 0.1743 0.2177 0.3118
mcmc_hist(samp,pars = c("theta1","theta2","theta3","theta4"),binwidth = 0.05)

mcmc_intervals(samp,pars = c("theta1","theta2","theta3","theta4"))

mcmc_combo(samp,pars = c("theta1","theta2","theta3","theta4"))

samp.dat <- as.data.frame(mcmc.list(samp)[[1]])
summary(samp.dat)
##      theta1           theta2           theta3            theta4       
##  Min.   :0.0451   Min.   :0.1057   Min.   :0.06498   Min.   :0.03151  
##  1st Qu.:0.2676   1st Qu.:0.3760   1st Qu.:0.36876   1st Qu.:0.13516  
##  Median :0.3507   Median :0.4487   Median :0.46131   Median :0.17431  
##  Mean   :0.3572   Mean   :0.4506   Mean   :0.46364   Mean   :0.17943  
##  3rd Qu.:0.4396   3rd Qu.:0.5236   3rd Qu.:0.55689   3rd Qu.:0.21773  
##  Max.   :0.8299   Max.   :0.8386   Max.   :0.89051   Max.   :0.46727
icbm_F_post<-samp.dat$theta1
icbm_S_post<-1-samp.dat$theta1
thaad_F_post<-samp.dat$theta2
thaad_H_post<-1-samp.dat$theta2
gmd_F_post<-samp.dat$theta3
gmd_H_post<-1-samp.dat$theta3
aegis_F_post<-samp.dat$theta4
aegis_H_post<-1-samp.dat$theta4

# Markov chain states for 3 layers BMD system and ICBM as target
pS1S5<-icbm_F_post # KN-11 failed to start successfully
pS1S2<-icbm_S_post # KN-11 started successfully and went boosting
pS2S3<-aegis_F_post # AEGIS failed to hit KN-11 on boost and ICBM went to midterm phase
pS2S5<-aegis_H_post # AEGIS  hit KN-11 on boost 
pS3S4<-gmd_F_post # GMD failed to hit KN-11 on midterm and ICBM went reentry phase
pS3S5<-gmd_H_post # GMD  hit KN-11 on midterm
pS4S6<-thaad_F_post # THAAD failed to hit KN-11 RV on impact - ICBM passed BMD
pS4S5<-thaad_H_post # THAAD hit KN-11 RV on impact - ICBM failed passing BMD

Markov chain model

KN-11 killed by 1 AEGIS or 1 GBD or 1 THAAD

res1<-c()
res2<-c()

for (i in 1:5000) {
  
  MC<-matrix(c(0.00,pS1S2[i],0,0,pS1S5[i],0,0,0,pS2S3[i],0,pS2S5[i],0,0,0,0,pS3S4[i],pS3S5[i],0,0,0,0,0,pS4S5[i],pS4S6[i],0,0,0,0,1,0,0,0,0,0,0,1),nrow=6,ncol=6,byrow=T)
  # MC
  MCT<-t(MC)
  # MCT
  
  res<-solve(diag(nrow=4,ncol=4)-MC[1:4,1:4])%*%MC[1:4,5:6]
  # print(res[1,2])
  res2[i]<-res[1,2]
  res1[i]<-res[1,1]
  
}
i
## [1] 5000
summary(res2)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.001502 0.014274 0.020909 0.023903 0.030329 0.106657
summary(res1)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.8933  0.9697  0.9791  0.9761  0.9857  0.9985
hist(res2, col="lightblue", main="Successfull impact of KN-11 RV", xlab = "Probability", freq = F)

hist(res1, col="red", main="KN-11 failed to pass through BMD with 1 AEGIS + 1 GBD + 1 THAAD", xlab = "Probability", freq = F)

stateNames <- c("S1","S2","S3", "S4", "S5", "S6")
row.names(MCT) <- stateNames; colnames(MC) <- stateNames
plotmat(MCT,pos = NULL, 
        lwd = 1, box.lwd = 2, 
        cex.txt = 0.8, 
        box.size = 0.1, 
        box.type = "circle", 
        box.prop = 0.5,
        box.col = "light yellow",
        arr.length=.1,
        arr.width=.1,
        self.cex = .4,
        self.shifty = -.01,
        self.shiftx = .13,
        main = "Markov chain for BMD: KN-11 killed by 1 AEGIS or 1 GBD or 1 THAAD")

For this case (1 AEGIS + 1 GBD + 1 THAAD) we have got 95% credible interval for the event KN-11 failed to pass through BMD as \(0.8933{\le}P_{KN11*}{\le}0.9985\) where \((P_{KN11*}>0.95)=0.947\) and \((P_{KN11*}<0.95)=0.053\)

KN-11 killed by 2 AEGIS or 2 GBD or 1 THAAD

res1<-c()
res2<-c()

for (i in 1:5000) {
  
  MC<-matrix(c(0.00,pS1S2[i],0,0,pS1S5[i],0,0,0,(pS2S3[i])^2,0,1-(pS2S3[i])^2,
               0,0,0,0,(pS3S4[i])^2,1-(pS3S4[i])^2,0,0,0,0,0,pS4S5[i],pS4S6[i],0,0,0,0,1,0,0,0,0,0,0,1),nrow=6,ncol=6,byrow=T)
  # MC
  MCT<-t(MC)
  # MCT
  
  res<-solve(diag(nrow=4,ncol=4)-MC[1:4,1:4])%*%MC[1:4,5:6]
  # print(res[1,2])
  res2[i]<-res[1,2]
  res1[i]<-res[1,1]
  
}
i
## [1] 5000
summary(res2)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 2.687e-05 8.123e-04 1.580e-03 2.376e-03 3.021e-03 2.510e-02
summary(res1)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.9749  0.9970  0.9984  0.9976  0.9992  1.0000
hist(res2, col="lightblue", main="Successfull impact of KN-11 RV", xlab = "Probability", freq = F)

hist(res1, col="red", main="KN-11 failed to pass through BMD with 2 AEGIS + 2 BMD + 1 THAAD", xlab = "Probability", freq = F)

stateNames <- c("S1","S2","S3", "S4", "S5", "S6")
row.names(MCT) <- stateNames; colnames(MC) <- stateNames
plotmat(MCT,pos = NULL, 
        lwd = 1, box.lwd = 2, 
        cex.txt = 0.8, 
        box.size = 0.1, 
        box.type = "circle", 
        box.prop = 0.5,
        box.col = "light yellow",
        arr.length=.1,
        arr.width=.1,
        self.cex = .4,
        self.shifty = -.01,
        self.shiftx = .13,
        main = "Markov chain for BMD: KN-11 killed by 2 AEGIS or 2 GBD or 1 THAAD")

For this case (2 AEGIS + 2 GBD + 1 THAAD) we have got 95% credible interval for the event KN-11 failed to pass through BMD as \(0.9749{\le}P_{KN11**}{\le}0.9999\) where \((P_{KN11**}>0.99)=0.9796\) and \((P_{KN11**}<0.99)=0.0204\)

So two AEGIS and two GBD missiles launched for one KN-11 target do matter even for three layers of BMD!

Conclusions

  1. Bayesian 95% credible interval for the event KN-11 failed to pass through BMD with 1 AEGIS or 1 GBD or 1 THAAD looks like \(0.8933{\le}P_{KN11*}{\le}0.9985\) where \((P_{KN11*}>0.95)=0.947\) and \((P_{KN11*}<0.95)=0.053\).
  2. Bayesian 95% credible interval for the event KN-11 failed to pass through BMD with 2 AEGIS or 2 GBD or 1 THAAD looks like \(0.9749{\le}P_{KN11**}{\le}0.9999\) where \((P_{KN11**}>0.99)=0.9796\) and \((P_{KN11**}<0.99)=0.0204\).
  3. North Korean missile threat is not a fake story though KN-11 and KN-14 reliability is not 100% case.
  4. US BMD statistics for test events (AEGIS, BMD, THAAD) should be updated if one wants to have real 95% credible interval for reliability of this system in particular case.

References

  1. https://warontherocks.com/2017/10/deadly-overconfidence-trump-thinks-missile-defenses-work-against-north-korea-and-that-should-scare-you/↩

  2. https://rpubs.com/alex-lev/269151↩

  3. https://rpubs.com/alex-lev/228238↩

  4. https://www.defense.gov/News/Special-Reports/BMDR/↩