Preface

North Korea has released pictures and videos of its controversial launch of a new missile which it claims can hit anywhere in the mainland United States¹.

Problem

We want to estimate Hwasong-15 Range (km), applying open source data on SLBM ² and Hwasong-15 ³.

Data

We use aggregated and transformed SLBM data to explore the relationship between M - Mass (kg), R - Range (km), S - Stage (1,2,3), D - Diameter (m), L - Length (m) of ballistic missile. See the previous topic on KN-11 ⁴.

load(file = "slbm.dat")
attach(slbm)
slbm[1:3,]

##   S     M    L    D    R    P W
## 1 1  5400 10.4 0.58  150  975 1
## 2 1 13700 11.8 1.30  650 1597 1
## 3 1 19650 14.2 1.30 1420 1179 1

Bayesian linear regression

This time we’ll use JAGS framework for MCMC with Bayesian estimates of posterior ⁵,⁶.

Model

model_string <- "model{

# Likelihood

for(i in 1:n){
R[i]   ~ dnorm(mu[i],inv.var)
mu[i] <- beta[1] + beta[2]*S[i] + beta[3]*D[i] + beta[4]*M[i]    
}

# Prior for beta
for(j in 1:4){
beta[j] ~ dnorm(0,0.0001)
}

# Prior for the inverse variance
inv.var   ~ dgamma(0.001, 0.001)
sigma     <- 1/sqrt(inv.var)

}"

Posterior sample

library(rjags)

## Loading required package: coda

## Linked to JAGS 4.3.0

## Loaded modules: basemod,bugs

library(bayesplot)

## This is bayesplot version 1.4.0

## - Plotting theme set to bayesplot::theme_default()

## - Online documentation at mc-stan.org/bayesplot

n <- length(slbm[,1])
model <- jags.model(textConnection(model_string), 
                    data = list(R=R,S=S,D=D,M=M,n=n),
                    inits=list(.RNG.name= "base::Wichmann-Hill",
                    .RNG.seed= 12341))

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 26
##    Unobserved stochastic nodes: 5
##    Total graph size: 173
## 
## Initializing model

update(model, 10000, progress.bar="none"); # Burnin for 10000 samples

num.iter <- 20000
samp <- coda.samples(model, 
                     variable.names=c("beta","sigma"), 
                     n.iter=num.iter, 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
## beta[1]   20.9494  99.75493 0.7053739      0.7378177
## beta[2]   53.2238  99.64124 0.7045700      0.8491243
## beta[3]   38.6979  98.79655 0.6985971      0.7893841
## beta[4]    0.1678   0.01544 0.0001092      0.0001439
## sigma   2193.0419 332.89736 2.3539398      2.5253523
## 
## 2. Quantiles for each variable:
## 
##             2.5%       25%       50%       75%     97.5%
## beta[1] -175.809  -45.2761   20.3851   88.8852  217.2929
## beta[2] -142.302  -13.5470   52.5946  121.2410  246.0507
## beta[3] -154.500  -27.9163   38.5901  104.8859  232.7657
## beta[4]    0.138    0.1575    0.1677    0.1779    0.1984
## sigma   1659.235 1956.5066 2152.5594 2386.1029 2946.0285

plot(samp)

mcmc_hist(samp, par=c("beta[1]", "beta[2]","beta[3]","sigma"))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

samp.dat <- as.data.frame(mcmc.list(samp)[[1]])
str(samp.dat)

## 'data.frame':    20000 obs. of  5 variables:
##  $ beta[1]: num  48.8 -131.6 -12.5 117 -10.4 ...
##  $ beta[2]: num  -178.3 17.7 41 -136 36.5 ...
##  $ beta[3]: num  95.3 -20.6 148.9 164.1 105.6 ...
##  $ beta[4]: num  0.188 0.185 0.167 0.169 0.181 ...
##  $ sigma  : num  2442 2365 2120 1880 2126 ...

B1 <- samp.dat$`beta[1]`
B2 <- samp.dat$`beta[2]`
B3 <- samp.dat$`beta[3]`
B4 <- samp.dat$`beta[4]`
SS <- samp.dat$sigma

RR <- c()
for (i in 1:10000) {
RR[i] <- B1[i] + B2[i]*3 + B3[i]*2.4 + B4[i]*70000
}

hist(RR,col = "red",main = "95% credible Hwasong-15 range distribution",xlab = "Range, km")

SE <- SS/sqrt(num.iter)

hist(SE,col = "blue",main = "95% credible interval for Standard Error",xlab = "SE, km")

quantile(RR,c(0.025,0.975))

##     2.5%    97.5% 
## 10171.41 13900.51

quantile(RR+2*SE,c(0.025,0.975))

##     2.5%    97.5% 
## 10202.67 13934.92

quantile(RR-2*SE,c(0.025,0.975))

##     2.5%    97.5% 
## 10140.14 13872.78

mean(RR)

## [1] 12023.88

Applying Bayesian regression we have got 95% credible interval for Hwasong-15 Range as \(P( 10171{\le}R_{HW15}{\le}13900)=0.95\) or including \(2*SE\) interval i.e. \(P(10140{\le}R_{HW15}{\le}13934)=0.95\) ⁷.

Conclusions

Applying Bayesian regression we have got 95% credible interval for Hwasong-15 Range as \(P(10140{\le}R_{HW15}{\le}13934)=0.95\) with the mean value \(M(R_{HW15})=12023\) km, provided \(D=2.4\) m, \(W=70\) t, \(S=3\).
Bayesian regression provides more robust and credible results for Hwasong-15 Range.
North Korean missile threat to the USA is more than real ⁸.

Hwasong-15

Alexander Levakov, Senior Research Fellow, Ph.D

December, 2017