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).
C3I multiplier for Trident D5
Alexander Levakov, Senior Research Fellow, Ph.D
March, 2017
“A piece of paper makes you an officer, a radio makes you a commander.” Omar N. Bradley
Preface
In the previous topic we have made comparison for MK4 and MK4A kill probability given successful launch of SLBM Trident D5.
See https://rpubs.com/alex-lev/257362
Problem
Now we want to estimate C3I multiplier i.e. probability of target would be killed by SLBM given SSBN correctly received EAM from TACAMO and successfully launched SLBM.
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.
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
data {
prec <- 1/(0.01^2) # precision
prec2 <- 1/(0.0001^2) # precision
NSLBM <- 24 # number of SLBM onboard SSBN
NMIRV <- 12 # number of MIRV onboard SLBM
N0 <- 3 # number of EAM repeats in communications session
N1 <- 2 # number of EAM correct received in communications session
EAM <- 256 # number of EAM characters
}
model {
PS0 ~ dbeta(23,4) # prior for probability of success launch of SLBM in 1990
NSL0 ~ dbin(PS0 , NSLBM) # number of success launches of 24 SLBM in 1990
PS1 ~ dbeta(155,6) # prior for for probability of success launch of SLBM in 2016
NSL1 ~ dbin(PS1 , NSLBM) # number of success launches of 24 SLBM in 2016
PK0 ~ dnorm(0.5,prec) # prior for kill probability of MIRV in 1990
NWK0 ~ dbin(PK0 , NMIRV) # number of targets killed by SLBM in 1990
PK1 ~ dnorm(0.85,prec) # prior for kill probability of MIRV in 2016
NWK1 ~ dbin(PK1 , NMIRV) # number of targets killed by SLBM in 2016
P1 ~ dnorm(0.9995,prec2) # prior for probability of correctly received character of EAM
PEAM <- pow(P1,EAM) # probability for correctly received EAM
NEAM ~ dbin(PEAM, N0) # number of correctly received EAM
NEAMS <- step(NEAM - 2.5) # number of correctly received EAM more than twice
PEAMS ~ dbeta(N1,N0-N1) # prior for probability of N1 received EAM out of N0 attempts
PEAMS2 <- PEAM*PEAMS # integral probability of correctly and authenticated EAM
PSLK0 <- PS0*PK0 # kill probability given success launch of SLBM without EAM in 1990
PSLK1 <- PS1*PK1 # kill probability given success launch of SLBM without EAM in 2016
PSLK0_EAM <- PSLK0*PEAMS2 # kill probability given success launch of SLBM and correct EAM in 1990
PSLK1_EAM <- PSLK1*PEAMS2 # kill probability given success launch of SLBM and correct EAM in 2016
NWKP0~dbin(PSLK0,NSL0*NWK0) # number of targets killed by SSBN given success launch of SLBM without EAM in 1990
NWKP1~dbin(PSLK1,NSL1*NWK1) # number of targets killed by SSBN given success launch of SLBM without EAM in 2016
NWKP0_EAM~dbin(PSLK0_EAM,NSL0*NWK0) # number of targets killed by SSBN given success launch of SLBM with EAM in 1990
NWKP1_EAM~dbin(PSLK1_EAM,NSL1*NWK1) # number of targets killed by SSBN given success launch of SLBM with EAM in 2016
}
Posterior sample
library(rjags)## Loading required package: coda## Linked to JAGS 4.1.0## Loaded modules: basemod,bugslibrary(bayesplot)## This is bayesplot version 1.1.0model <- jags.model(file = "model_example_3.bug", n.chains=4)## Compiling data graph
## Resolving undeclared variables
## Allocating nodes
## Initializing
## Reading data back into data table
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 0
## Unobserved stochastic nodes: 15
## Total graph size: 44
##
## Initializing modelupdate(model, 1000, progress.bar="text"); # Burnin for 10000 samples
samp <- coda.samples(model,
variable.names=c("NSL0","NSL1","NWK0","NWK1",
"NWKP0", "NWKP1", "NWKP0_EAM",
"NWKP1_EAM","NEAM","NEAMS",
"P1","PEAM","PEAMS","PEAMS2",
"PS0","PS1","PK0","PK1", "PSLK0","PSLK1","PSLK0_EAM","PSLK1_EAM"),
n.iter=5000, progress.bar="text")
summary(samp)##
## Iterations = 1001:6000
## Thinning interval = 1
## Number of chains = 4
## Sample size per chain = 5000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## NEAM 2.6425 5.620e-01 3.974e-03 3.934e-03
## NEAMS 0.6836 4.651e-01 3.289e-03 3.263e-03
## NSL0 20.4455 2.334e+00 1.650e-02 1.641e-02
## NSL1 23.1193 9.791e-01 6.923e-03 6.932e-03
## NWK0 6.0039 1.730e+00 1.223e-02 1.211e-02
## NWK1 10.2052 1.256e+00 8.882e-03 8.453e-03
## NWKP0 52.5481 1.864e+01 1.318e-01 1.299e-01
## NWKP0_EAM 30.8392 1.631e+01 1.153e-01 1.144e-01
## NWKP1 193.1895 2.668e+01 1.886e-01 1.874e-01
## NWKP1_EAM 113.2778 4.412e+01 3.120e-01 3.140e-01
## P1 0.9995 1.001e-04 7.078e-07 7.116e-07
## PEAM 0.8802 2.257e-02 1.596e-04 1.605e-04
## PEAMS 0.6659 2.360e-01 1.669e-03 1.646e-03
## PEAMS2 0.5862 2.084e-01 1.474e-03 1.453e-03
## PK0 0.5000 1.001e-02 7.080e-05 7.080e-05
## PK1 0.8500 9.976e-03 7.054e-05 7.055e-05
## PS0 0.8520 6.617e-02 4.679e-04 4.744e-04
## PS1 0.9630 1.472e-02 1.041e-04 1.047e-04
## PSLK0 0.4260 3.414e-02 2.414e-04 2.413e-04
## PSLK0_EAM 0.2496 9.112e-02 6.443e-04 6.480e-04
## PSLK1 0.8185 1.575e-02 1.114e-04 1.104e-04
## PSLK1_EAM 0.4797 1.708e-01 1.208e-03 1.191e-03
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## NEAM 1.00000 2.0000 3.0000 3.0000 3.0000
## NEAMS 0.00000 0.0000 1.0000 1.0000 1.0000
## NSL0 15.00000 19.0000 21.0000 22.0000 24.0000
## NSL1 21.00000 23.0000 23.0000 24.0000 24.0000
## NWK0 3.00000 5.0000 6.0000 7.0000 9.0000
## NWK1 7.00000 9.0000 10.0000 11.0000 12.0000
## NWKP0 19.00000 39.0000 51.0000 65.0000 92.0000
## NWKP0_EAM 5.00000 19.0000 29.0000 41.0000 67.0000
## NWKP1 136.00000 176.0000 195.0000 213.0000 241.0000
## NWKP1_EAM 25.00000 82.0000 117.0000 147.0000 189.0000
## P1 0.99931 0.9994 0.9995 0.9996 0.9997
## PEAM 0.83717 0.8648 0.8800 0.8953 0.9253
## PEAMS 0.15323 0.5022 0.7059 0.8645 0.9865
## PEAMS2 0.13470 0.4410 0.6216 0.7604 0.8754
## PK0 0.48017 0.4933 0.5001 0.5067 0.5194
## PK1 0.83054 0.8433 0.8500 0.8567 0.8694
## PS0 0.70077 0.8122 0.8599 0.9013 0.9560
## PS1 0.92922 0.9542 0.9649 0.9737 0.9860
## PSLK0 0.34919 0.4051 0.4296 0.4507 0.4818
## PSLK0_EAM 0.05769 0.1860 0.2619 0.3218 0.3923
## PSLK1 0.78513 0.8085 0.8194 0.8296 0.8467
## PSLK1_EAM 0.11089 0.3608 0.5093 0.6217 0.7206Results
samp.dat <- as.data.frame(mcmc.list(samp)[[1]])
NSL0=samp.dat$NSL0
NSL1=samp.dat$NSL1
PS0=samp.dat$PS0
PS1=samp.dat$PS1
NWK0=samp.dat$NWK0
NWK1=samp.dat$NWK1
PK0=samp.dat$PK0
PK1=samp.dat$PK1
NWKP0=samp.dat$NWKP0
NWKP1=samp.dat$NWKP1
PSLK0=samp.dat$PSLK0
PSLK1=samp.dat$PSLK1
NWKP0_EAM=samp.dat$NWKP0_EAM
NWKP1_EAM=samp.dat$NWKP1_EAM
PSLK0_EAM=samp.dat$PSLK0_EAM
PSLK1_EAM=samp.dat$PSLK1_EAM
NEAMS=samp.dat$NEAMS
NEAM=samp.dat$NEAM
PEAM=samp.dat$PEAM
PEAMS=samp.dat$PEAMS
PEAMS2=samp.dat$PEAMS2
color_scheme_set("brightblue")
mcmc_intervals(samp, par=c("P1", "PEAM","PEAMS","PEAMS2"))
color_scheme_set("green")
mcmc_intervals(samp, par=c("PS0","PS1","PK0","PK1"))
color_scheme_set("blue")
mcmc_intervals(samp, par=c("NSL0","NSL1","NWK0","NWK1"))
color_scheme_set("red")
mcmc_intervals(samp, par=c("NWKP0", "NWKP1", "NWKP0_EAM", "NWKP1_EAM"))
par(mfrow=c(2,2))
hist(NWKP0_EAM,freq=F,col="red")
lines(density(NWKP0_EAM))
plot(PS0,PSLK0_EAM)
plot(PS1,PSLK1_EAM)
hist(NWKP1_EAM,freq=F,col="red")
lines(density(NWKP1_EAM))
par(mfrow=c(2,2))
hist(NWKP0,freq=F,col="blue")
lines(density(NWKP0))
plot(PS0,PSLK0)
plot(PS1,PSLK1)
hist(NWKP1,freq=F,col="blue")
lines(density(NWKP1))
par(mfrow=c(2,2))
hist(NSL0,xlim = c(10,24),breaks = 10,main = "1990", xlab = "Number of SLBM success by SSBN",col="blue",freq=F)
hist(NSL1,xlim = c(16,24),breaks = 5,main = "2016", xlab = "Number of of SLBM success by SSBN",col="red",freq=F)
hist(NWK0,breaks = 12,main = "1990", xlab = "Number of targets killed by SLBM",col="blue",freq=F)
hist(NWK1,breaks = 6,main = "2016", xlab = "Number of targets killed by SLBM",col="red",freq=F)
par(mfrow=c(2,2))
hist(NWKP0,breaks = 12,main = "1990", xlab = "Number of targets killed by SSBN",col="blue",freq=F)
hist(NWKP1,breaks = 12,main = "2016", xlab = "Number of targets killed by SSBN",col="red",freq=F)
hist(NWKP0_EAM,breaks = 12,main = "1990", xlab = "Number of targets killed by SSBN with EAM",col="blue",freq=F)
hist(NWKP1_EAM,breaks = 12,main = "2016", xlab = "Number of targets killed by SSBN with EAM",col="red",freq=F)
The effect of EAM
Applying results above we can get 95% credible interval for EAM effect as C3I multiplier:
REAM0 <- NWKP0_EAM/NWKP0
REAM1 <- NWKP1_EAM/NWKP1
hist(REAM0, col="blue",freq = F,main = "EAM effect in 1990")
quantile(REAM0,c(0.025,0.975),na.rm = T)## 2.5% 97.5%
## 0.125000 1.047619hist(REAM1, col="red",freq = F,main = "EAM effect in 2016")
quantile(REAM1,c(0.025,0.975),na.rm = T)## 2.5% 97.5%
## 0.1443713 0.9027046Conclusions
- The 95% credible interval for probability of correctly receiving EAM with 256 characters no less than 2 times out of 3 repeats in the ambient environment looks like: \(0.14115\le P_{EAM,2/3}\le 0.8743\) with the mean value \(0.6226\) given probability of correct character \(0.99930\le P_{1} \le0.9997\) with the mean \(0.9995\).
- The 95% credible interval for probability of target killed by MIRV given SLBM successful launch and correctly received EAM in 1990 and 2016 looks like: \(0.06007 \le P_{KMIRV|S,EAM,1990}\le 0.3918\) with the mean value \(0.2629\) and \(0.11546 \le P_{KMIRV|S,EAM,2016}\le 0.7187\) with the mean value \(0.5092\).
- The 95% credible interval for number of targets killed by SSBN given SLBM successful launch and correctly received EAM in 1990 and 2016 looks like: \(5 \le N_{KSSBN|S,EAM,1990}\le 67\) with the mean value \(29\) and \(26 \le N_{KSSBN|S,EAM,2016}\le 189\) with the mean value \(117\).
- The 95% credible interval for the effect of EAM in 1990 and 2016 looks like: \(0.1249298 \le R_{1990}\le 1.0492213\) and \(0.1563933 \le R_{2016}\le 0.9045693\).
- C3I does matter in evaluating SSBN target kill potential!