Preface

The Bulletin of atomic scientists published some interesting paper titled “How US nuclear force modernization is undermining strategic stability: The burst-height compensating super-fuze”. See http://thebulletin.org/how-us-nuclear-force-modernization-undermining-strategic-stability-burst-height-compensating-super10578

It seems very strange to read this paper after real news stories about failed UK and US Trident D5 flights. We’ll try to explore some grounds of this paper using Bayesian statistics together with MCMC framework.

Problem

We want to estimate progress in strategic nuclear programs made by US DOD in recent years comparing kill target potential of US submarine with 24 Trident D5 armed with MK4A and MK4 MIRV s.

Data

We use open source Trident D5 flight test data (https://rpubs.com/alex-lev/200551) together with BAS diagram for MK4/MK4A kill probability (above).

Main assumptions

We consider random experiment with exactly two possible outcomes, “success” and “failure”, in which the probability of success is the same every time the experiment is conducted. (http://en.wikipedia.org/wiki/Bernoulli_trial).

Point estimate

Given BAS diagram we can make logical inference about SSBN kill potential with 24 D5 each with 12 MIRV as \(N_{MK4}=24*12*0.5=144\) and \(N_{MK4A}=24*12*0.85=244.8\). That means absolute increase by 100 or 70% killed targets for one of 14 US nuclear SSBNs. Good news for Donald Trump as commander in chief who is facing North Korean and Iranian ICBM threat! The bad news is that this point estimate is made with assumption of 100% reliability of SLBM launch.

Interval estimate

Frequentist case

Here we postulate null hypothesis about SLBM (Trident D5) reliability \(H_{0}=0.95\) and alternative \(H_{1}\ne0.95\). Let’s prove it for two time periods: post flight test and current.

fit1 <- binom.test(23,27,0.95)
fit2 <- binom.test(155,161,0.95)
fit1
## 
##  Exact binomial test
## 
## data:  23 and 27
## number of successes = 23, number of trials = 27, p-value = 0.04374
## alternative hypothesis: true probability of success is not equal to 0.95
## 95 percent confidence interval:
##  0.6626891 0.9581126
## sample estimates:
## probability of success 
##              0.8518519
fit2
## 
##  Exact binomial test
## 
## data:  155 and 161
## number of successes = 155, number of trials = 161, p-value =
## 0.5874
## alternative hypothesis: true probability of success is not equal to 0.95
## 95 percent confidence interval:
##  0.9206528 0.9862029
## sample estimates:
## probability of success 
##              0.9627329

We got two 95% confident intervals for D5 reliability: \(0.663\le P_{0}\le0.958\) and \(0.921\le P_{1}\le0.986\), where \(P_0\) is the reliability for the post flight test interval and \(P_1\) - current. In fact we failed to prove null \(H_0:P_0=0.95\). Anyway we can proceed our inference with these interval estimates to produce mean expected value for killed targets. Let’s do it.

NMK4low <- fit1$conf.int[1]*0.5*12*24
NMK4upper <- fit1$conf.int[2]*0.5*12*24
NMK4Alow <- fit2$conf.int[1]*0.85*12*24
NMK4Aupper <- fit2$conf.int[2]*0.85*12*24
NMK4low
## [1] 95.42724
NMK4upper
## [1] 137.9682
NMK4Alow
## [1] 225.3758
NMK4Aupper
## [1] 241.4225

So we got two intervals for ME: \(95\le N_{MK4}\le138\) and \(225\le N_{MK4A}\le241\), where \(N_{MK4}\) is the ME for the post flight test interval and \(MK4A\) - current. In fact the main question is the value of probability that no less than 216 (75%) out of 288 MIRV s would hit and kill targets. Let’s do calculations.

sum(dbinom(216:288,288,fit1$conf.int[1]*0.5))
## [1] 8.74943e-48
sum(dbinom(216:288,288,fit1$conf.int[2]*0.5))
## [1] 6.284611e-21
sum(dbinom(216:288,288,fit2$conf.int[1]*0.85))
## [1] 0.9190525
sum(dbinom(216:288,288,fit2$conf.int[2]*0.85))
## [1] 0.9999548

Now we have produced very interesting estimates of the probability \(P({216 \le N\le 288})\) for two periods: \(P({216 \le N_{MK4}\le 288})=0\) and \(0.414\le P({216 \le N_{MK4A}\le 288})\le0.998\). So it’s time to ask Bayes!

Bayesian case

This time we’ll use JAGS framework for MCMC with Bayesian estimates of posterior. For more see http://mcmc-jags.sourceforge.net/.

Model

model {

prec <- 1/(0.001^2) # precision as 1/sigma**2

# priors for binomial and normal distribution

theta1 ~ dbeta(155,6) # current D5 reliability

theta2 ~ dbeta(23,4) # initial D5 reliability

y1 ~ dbin(theta1,N1) # success for current D5 reliability

y2 ~ dbin(theta2,N2) # success for initial D5 reliability

y1_ ~ dnorm(0.85,prec) # MK4A kill prior

y2_ ~ dnorm(0.5,prec) # MK4 kill prior

# posterior distributions

y11 <- theta1*y1_ # posterior for current D5 + MK4A reliability

y22 <- theta2*y2_ # posterior for initial D5 + MK4 reliability

NW1 ~ dbin(y11,288) # posterior for number MK4A would kill target

NW2 ~ dbin(y22,288) # posterior for number MK4 would kill target

# odds ratios for priors and posteriors

oddsratio = (theta1/(1-theta1))/(theta2/(1-theta2))

oddsratio1 = (y1_/(1-y1_))/(y2_/(1-y2_))

oddsratio2 = (y11/(1-y11))/(y22/(1-y22))

oddsratio3 = NW1/NW2
}

Posterior sample

library(rjags)
## Loading required package: coda
## Linked to JAGS 4.1.0
## Loaded modules: basemod,bugs
library(bayesplot)
## This is bayesplot version 1.1.0
N1 <- 24 #number of SLBM to be launched by SSBN (100%)
N2<- 24 #number of SLBM to be launched by SSBN (100%)

model <- jags.model(file = "model.bug", data = list(N1=N1,N2=N2),n.chains=3)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 0
##    Unobserved stochastic nodes: 8
##    Total graph size: 67
## 
## Initializing model
update(model, 10000, progress.bar="none"); # Burnin for 10000 samples

samp <- coda.samples(model, 
        variable.names=c("theta1","theta2",
                         "y1_","y2_","y11","y22",                                           "NW1","NW2",
                               "oddsratio","oddsratio1",
                               "oddsratio2","oddsratio3"), 
        n.iter=30000, progress.bar="none")

summary(samp)
## 
## Iterations = 10001:40000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 30000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                Mean        SD  Naive SE Time-series SE
## NW1        235.6878 7.486e+00 2.495e-02      2.505e-02
## NW2        122.6254 1.276e+01 4.253e-02      4.238e-02
## oddsratio    5.6280 4.497e+00 1.499e-02      1.502e-02
## oddsratio1   5.6672 4.955e-02 1.652e-04      1.657e-04
## oddsratio2   6.1765 1.050e+00 3.501e-03      3.519e-03
## oddsratio3   1.9443 2.256e-01 7.519e-04      7.526e-04
## theta1       0.9627 1.491e-02 4.971e-05      4.971e-05
## theta2       0.8519 6.724e-02 2.241e-04      2.243e-04
## y11          0.8183 1.271e-02 4.237e-05      4.237e-05
## y1_          0.8500 9.936e-04 3.312e-06      3.325e-06
## y22          0.4259 3.363e-02 1.121e-04      1.122e-04
## y2_          0.5000 9.962e-04 3.321e-06      3.333e-06
## 
## 2. Quantiles for each variable:
## 
##                2.5%      25%      50%      75%    97.5%
## NW1        220.0000 231.0000 236.0000 241.0000 250.0000
## NW2         96.0000 114.0000 123.0000 132.0000 146.0000
## oddsratio    1.0198   2.7464   4.4359   7.0834  17.4513
## oddsratio1   5.5708   5.6339   5.6667   5.7004   5.7653
## oddsratio2   4.5615   5.4512   6.0214   6.7321   8.6624
## oddsratio3   1.5944   1.7879   1.9120   2.0672   2.4742
## theta1       0.9284   0.9539   0.9646   0.9736   0.9861
## theta2       0.6984   0.8112   0.8609   0.9015   0.9568
## y11          0.7891   0.8108   0.8199   0.8276   0.8383
## y1_          0.8481   0.8493   0.8500   0.8507   0.8520
## y22          0.3491   0.4056   0.4305   0.4508   0.4784
## y2_          0.4980   0.4993   0.5000   0.5007   0.5020
color_scheme_set("brightblue")
#posterior distribution for number of targets killed given D5 #successful launch: NW1 - current D5 reliability with MK4A, 
#NW1 - initial D5 reliability with MK4 
mcmc_intervals(samp, par=c("NW1","NW2"))

#posterior distribution for probability target would be killed 
#given D5 successful launch: y11 - current D5 reliability with MK4A, y22 - initial D5 reliability with MK4 
mcmc_intervals(samp, par=c("y11","y22"))

#prior distribution for probability  D5 successful launch:
#theta1 - current D5 reliability, theta2 - initial D5 reliability 
mcmc_intervals(samp, par=c("theta1","theta2"))

#prior distribution for probability  MIRV would kill target:
#y1_ - MK4A, y2_ - MK4 
mcmc_intervals(samp, par=c("y1_","y2_"))

color_scheme_set("pink")
mcmc_hist(samp, binwidth=2, par=c("NW1","NW2"))

mcmc_hist(samp, par=c("y11","y22"))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

mcmc_hist(samp, par=c("theta1","theta2"))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

mcmc_hist(samp, par=c("y1_","y2_"))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

mcmc_hist(samp, par=c("oddsratio","oddsratio1","oddsratio2","oddsratio3"))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#color_scheme_set("blue")
#scatterplot for correlations
#mcmc_scatter(samp, pars = c("y11","NW1"))
#mcmc_scatter(samp, pars = c("y22","NW2"))

#color_scheme_set("green")
#diagnostics for sampling
#mcmc_trace(samp, pars = c("NW1","NW2"))
#mcmc_trace(samp, pars = c("theta1","theta2"))


#sample transformation of variables array
samp.dat <- as.data.frame(mcmc.list(samp)[[1]])
Y11 <- as.array(samp.dat$y11)
Y22 <- as.array(samp.dat$y22)

NW1_ <- as.array(rep(0,5000))
NW2_ <- as.array(rep(0,5000))

#
for (i in 1:5000) {
  NW1_[i] <- sum(dbinom(216:288,288,Y11[i]))
}

hist(NW1_,xlim = c(0.98,1),breaks = 300,main = "MK4A, P(216<=N<=288)", xlab = "Probability", col="red")

summary(NW1_)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.5765  0.9956  0.9988  0.9937  0.9997  1.0000
quantile(NW1_,c(0.025,0.975))
##      2.5%     97.5% 
## 0.9502535 0.9999542
for (i in 1:5000) {
  NW2_[i] <- sum(dbinom(144:288,288,Y22[i]))
}

hist(NW2_,xlim = c(0.01,0.5),main = "MK4, P(144<=N<=288)", xlab = "Probability",col="blue")

summary(NW2_)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0000000 0.0007351 0.0103300 0.0420700 0.0535400 0.4450000
quantile(NW2_,c(0.025,0.975))
##         2.5%        97.5% 
## 8.919839e-08 2.527411e-01

Conclusions

  1. Applying frequentist methods we got 95% confident intervals for number of targets killed given successful launch: \(95\le N_{MK4}\le138\) and \(225\le N_{MK4A}\le241\).
  2. The value of probability that no less than 216 (75%) out of 288 MIRV s would kill targets has 95% confident intervals: \(P({216 \le N_{MK4}\le 288})=0\) and \(0.414\le P({216 \le N_{MK4A}\le 288})\le0.998\).
  3. Applying Bayesian methods we have got 95% credible intervals for number of target killed given successful launch: \(96\le N_{MK4}\le146\) and \(220\le N_{MK4A}\le250\).
  4. The value of probability that no less than 216 (75%) out of 288 MIRV s would kill targets \(P({216 \le N\le 288})\) has 95% credible interval for MK4A as \(0.9588866 \le P({216 \le N_{MK4A}\le 288})\le 0.9999554\). For MK4 interval \(144 \le N_{MK4}\le 288\) probability looks like \(0\le P({144 \le N_{MK4}\le 288})\le 0.2577\).
  5. The mean values for killed targets \(N_{MK4A}=236\), \(N_{MK4}=123\) give mean effect ratio \(N_{MK4A}/N_{MK4}=1.94\) with 95% credible interval \(1.593\le P(N_{MK4A}/N_{MK4})\le2.468\).
  6. The 95% credible intervals for kill probability given success launch look like: \(0.789\le P_{kill,D5+MK4A}\le 0.838\) and \(0.349\le P_{kill,D5+MK4}\le 0.478\).
  7. The truth is that all these estimates, both frequentist and Bayesian, have been made with assumption of 100% probability of correctly received EAM by SSBN radioman and that is bad news: who knows what signal to noise ratio would be under stressed environment for submerged nuclear submarine on patrol.

P.S. The most thrilling and exciting news is that Lawrence Livermore National Laboratory began reestimating nuclear weapons effects using old films of 210 atmospheric tests in 50-s and 60-s. Unfortunetly, that’s not a fake story!