library(binom)
library(pwr)
library(ggplot2)
Problem
We want to compare SLBM reliability in terms of probability of launch, provided we have got knowledge of some historical data for D5, M51, JL-2, SSNX30, KN-11.
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).
- We need some example of a real experiment. Let’s see it here: https://en.wikipedia.org/wiki/Submarine-launched_ballistic_missile.
- We suppose that the submarine on patrol being on high alert must be able to launch no less than 75% of missiles on board.
- We use Bayesian inference made by MCMC framework: https://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo.
Data
library(rjags)
## Loading required package: coda
## Linked to JAGS 4.2.0
## Loaded modules: basemod,bugs
library(bayesplot)
## This is bayesplot version 1.2.0
N1 <- 24 # number of D5 SLBM onboard SSBN
Y1 <- c(18:24) # number of D5 SLBM to be launched with success
N2 <- 16 # number of M51 SLBM onboard SSBN
Y2 <- c(12:16) # number of M51 SLBM to be launched with success
N3 <- 12 # number of JL2 SLBM onboard SSBN
Y3 <- c(9:12) # number of JL2 SLBM to be launched with success
N4 <- 16 # number of SSNX30 SLBM onboard SSBN
Y4 <- c(12:16) # number of SSNX30 SLBM to be launched with success
N5 <- 2 # number of KN-11 SLBM onboard SSBN
Y5 <- c(1,2) # number of KN-11 SLBM to be launched with success
a1 <- 155 #prior success for Trident D5 SLBM: 2016
b1 <- 6 #prior falure for Trident D5 SLBM: 2016
a2 <- 6 #prior success for M51 SLBM: 2015
b2 <- 1 #prior falure for M51 SLBM: 2015
a3 <- 5 #prior success for JL2 SLBM: 2015
b3 <- 1 #prior falure for JL2 SLBM: 2015
a4 <- 18 #prior success for SSNX30 SLBM: 2016
b4 <- 8 #prior falure for SSNX30 SLBM: 2016
a5 <- 8 #prior success for KN-11 SLBM: 2016
b5 <- 4 #prior falure for KN-11 SLBM: 2016
data_list <- list(Y1=Y1,N1=N1,Y2=Y2,N2=N2,Y3=Y3,N3=N3,Y4=Y4,N4=N4,Y5=Y5,N5=N5,
a1=a1,b1=b1,a2=a2,b2=b2,a3=a3,b3=b3,a4=a4,b4=b4,a5=a5,b5=b5) # data for the model
Model and posterior sampling
inits1 <- list(.RNG.name="base::Super-Duper", .RNG.seed=12345)
model <- jags.model(file = "binom.bug", data = data_list,n.chains = 4,inits = inits1)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 23
## Unobserved stochastic nodes: 5
## Total graph size: 47
##
## Initializing model
update(model, 10000, progress.bar="none"); # Burnin for 10000 samples
samp <- coda.samples(model,
variable.names=c("PD5","PM51","PJL2","PSSNX30","PKN11","diff1","diff2","diff3","diff4"),
n.iter=20000, progress.bar="none")
model
## JAGS model:
##
## model{
##
## for(i in 1:length(Y1)){
## Y1[i] ~ dbinom(PD5,N1) # posterior distribution of success trials for D5
## }
## PD5 ~ dbeta(a1, b1) # prior distribution for D5 probability of success
##
## for(i in 1:length(Y2)){
## Y2[i] ~ dbinom(PM51,N2) # posterior distribution of success trials for M51
## }
## PM51 ~ dbeta(a2, b2) # prior distribution for M51 probability of success
##
## for(i in 1:length(Y3)){
## Y3[i] ~ dbinom(PJL2,N3) # posterior distribution of success trials for JL2
## }
## PJL2 ~ dbeta(a3, b3) # prior distribution for JL2 probability of success
##
## for(i in 1:length(Y4)){
## Y4[i] ~ dbinom(PSSNX30,N4) # posterior distribution of success trials for SSNX30
## }
## PSSNX30 ~ dbeta(a4, b4) # prior distribution for SSNX30 probability of success
##
## for(i in 1:length(Y5)){
## Y5[i] ~ dbinom(PKN11,N5) # posterior distribution of success trials for KN-11
## }
## PKN11 ~ dbeta(a5, b5) # prior distribution for KN-11 probability of success
##
## diff1 = PD5 - PM51 # difference between posterior distributions of D5 and M51
## diff2 = PD5 - PJL2 # difference between posterior distributions of D5 and JL-2
## diff3 = PD5 - PSSNX30 # difference between posterior distributions D5 and SSNX30
## diff4 = PD5 - PKN11 # difference between posterior distributions D5 and KN-11
## }
## Fully observed variables:
## N1 N2 N3 N4 N5 Y1 Y2 Y3 Y4 Y5 a1 a2 a3 a4 a5 b1 b2 b3 b4 b5
summary(samp)
##
## Iterations = 10001:30000
## Thinning interval = 1
## Number of chains = 4
## 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
## PD5 0.91783 0.01517 5.365e-05 5.365e-05
## PJL2 0.86983 0.04562 1.613e-04 1.597e-04
## PKN11 0.68787 0.11236 3.973e-04 4.013e-04
## PM51 0.87328 0.03564 1.260e-04 1.260e-04
## PSSNX30 0.83043 0.03606 1.275e-04 1.275e-04
## diff1 0.04455 0.03858 1.364e-04 1.364e-04
## diff2 0.04800 0.04811 1.701e-04 1.685e-04
## diff3 0.08740 0.03920 1.386e-04 1.386e-04
## diff4 0.22996 0.11338 4.008e-04 4.039e-04
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## PD5 0.88634 0.90804 0.91851 0.92832 0.9454
## PJL2 0.76848 0.84155 0.87442 0.90331 0.9449
## PKN11 0.44899 0.61394 0.69581 0.76935 0.8824
## PM51 0.79574 0.85109 0.87635 0.89875 0.9348
## PSSNX30 0.75356 0.80736 0.83269 0.85572 0.8951
## diff1 -0.02447 0.01751 0.04240 0.06882 0.1277
## diff2 -0.03406 0.01346 0.04414 0.07874 0.1525
## diff3 0.01475 0.06014 0.08560 0.11277 0.1701
## diff4 0.03322 0.14774 0.22176 0.30490 0.4694
Results
summary(samp)
##
## Iterations = 10001:30000
## Thinning interval = 1
## Number of chains = 4
## 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
## PD5 0.91783 0.01517 5.365e-05 5.365e-05
## PJL2 0.86983 0.04562 1.613e-04 1.597e-04
## PKN11 0.68787 0.11236 3.973e-04 4.013e-04
## PM51 0.87328 0.03564 1.260e-04 1.260e-04
## PSSNX30 0.83043 0.03606 1.275e-04 1.275e-04
## diff1 0.04455 0.03858 1.364e-04 1.364e-04
## diff2 0.04800 0.04811 1.701e-04 1.685e-04
## diff3 0.08740 0.03920 1.386e-04 1.386e-04
## diff4 0.22996 0.11338 4.008e-04 4.039e-04
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## PD5 0.88634 0.90804 0.91851 0.92832 0.9454
## PJL2 0.76848 0.84155 0.87442 0.90331 0.9449
## PKN11 0.44899 0.61394 0.69581 0.76935 0.8824
## PM51 0.79574 0.85109 0.87635 0.89875 0.9348
## PSSNX30 0.75356 0.80736 0.83269 0.85572 0.8951
## diff1 -0.02447 0.01751 0.04240 0.06882 0.1277
## diff2 -0.03406 0.01346 0.04414 0.07874 0.1525
## diff3 0.01475 0.06014 0.08560 0.11277 0.1701
## diff4 0.03322 0.14774 0.22176 0.30490 0.4694
mcmc_areas(samp,pars = c("PD5","PM51","PJL2","PSSNX30","PKN11"))
mcmc_areas(samp,pars = c("diff1","diff2","diff3","diff4"))
mcmc_hist(samp,pars = c("PD5","PM51","PJL2","PSSNX30","PKN11"))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
mcmc_hist(samp,pars = c("diff1","diff2","diff3","diff4"))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
mcmc_intervals(samp,pars = c("PD5","PM51","PJL2","PSSNX30","PKN11"))
mcmc_intervals(samp,pars = c("diff1","diff2","diff3","diff4"))
Power of test
library(binom)
binom.power(p.alt = 0.85,n = 27,p = 0.5,alpha = 0.05) #D5
## alpha
## 0.9657898
binom.power(p.alt = 0.85,n = 7,p = 0.5,alpha = 0.05) #M51
## alpha
## 0.6569511
binom.power(p.alt = 0.85,n = 6,p = 0.5,alpha = 0.05) #JL2
## alpha
## 0.6172653
binom.power(p.alt = 0.85,n = 26,p = 0.5,alpha = 0.05) #SSNX30
## alpha
## 0.9614512
binom.power(p.alt = 0.85,n = 12,p = 0.5,alpha = 0.05) #KN-11
## alpha
## 0.8035035
binom.plot(n=30,type="levelplot")
## Loading required package: lattice
Conclusion
- US Trident D5 leads SLBM race in terms of reliability as 95% credible interval \(0.88634{\le}P_{D5}(18\le{X}{\le24}){\le}0.9454\) while North Korean KN-11 is outsider \(0.44899{\le}P_{KN11}(1\le{X}{\le2}){\le}0.8824\).
- French M51 \(0.79574{\le}P_{M51}(12\le{X}{\le16}){\le}0.9348\) and Chinese JL-2 \(0.76848 {\le}P_{JL2}(9\le{X}{\le12}){\le}0.9449\) are very close to each other.
- SSNX30 95% credible interval for reliability \(0.75356{\le}P_{SSNX30}(12\le{X}{\le16}){\le}0.8951\) is near M51 and JL-2.
- Power of binomial test (flight test) demonstrates that D5 and SSNX30 have more than 90% probability to see difference between \(H_0:P=0.5\) and \(H_1:P=0.85\).
- M51 together with JL-2 need more additional trials to produce reliable power of test.
- North Korean KN-11 has 80% power of test.