Problem
Main assumptions
Preliminary point estimate
Is it enough?
Interval estimates
Power of the test
Bayesian estimate
- Current results
- Future results
Conclusions

Ju Lang-2

library(binom)
library(pwr)
library(ggplot2)
library(rjags)

## Loading required package: coda

## Linked to JAGS 4.1.0

## Loaded modules: basemod,bugs

library(mcmc)
library(stringr)
source(file = "bayes_prp_tst.R")
source(file = "binom_test.R")
source(file = "r_jags.R")
source(file = "generic.R")

Problem

We want to estimate progress in SLBM program made by China in recent years compared to that made by France. So we choose chinese JL-2 and french M51 SLBMs.

M51

Main assumptions

We consider missile test flight as 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 data of a real experiment. Let’s see it here: https://missilethreat.csis.org/missile/jl-2/, https://missilethreat.csis.org/missile/m51/.

Preliminary point estimate

The data for JL-2 and M51 test flights is a real experiment with two possible outcomes: success and failure, relevant with Bernoulli trial. So far JL-2 has 6 trials with 5 successful outcomes and M51 - 7 trials with 6 successful outcomes. So we got our first point estimates of the probability of success for two SLBMs: \(P_{JL2}=5/6=0.83\) and \(P_{M51}=6/7=0.86\)

Is it enough?

Good question! As of 2016, China deploys four Jin-class nuclear ballistic missile submarines each armed with 12 JL-2 SLBMs.

SLBM underwater launch

We suppose that the submarine on patrol being on high alert must be able to launch no less than 90% of missiles on board. So we have got for chinese and french submarines low limits as \(Nch=12*0.9=11\) and \(Nfr=16*0.9=14\). What is the probability for successful launch of no less than 90% of SLBMs for chinese \(P(11<=X<=12)\) and french \(P(14<=X<=16)\) submarines?

sum(dbinom(11:12,12,0.83))

## [1] 0.3696076

sum(dbinom(14:16,16,0.86))

## [1] 0.6074478

The result is unsatisfactory for chinese \(P(11<=X<=12)=0.37\) and very modest though satisfactory for french one \(P(14<=X<=16)=0.61\).

Interval estimates

First, we postulate null H₀| P=0.95 and alternative H₁| P<0.95 hypothesis for JL-2 and M51. Then we go through the test to prove or reject null hypothesis, applying binom.test function.

binom.test(x = 5,n = 6,p = 0.95,alternative = "less")

## 
##  Exact binomial test
## 
## data:  5 and 6
## number of successes = 5, number of trials = 6, p-value = 0.2649
## alternative hypothesis: true probability of success is less than 0.95
## 95 percent confidence interval:
##  0.0000000 0.9914876
## sample estimates:
## probability of success 
##              0.8333333

Next, we repeat this test for M51.

binom.test(x = 6,n = 7,p = 0.95,alternative = "less")

## 
##  Exact binomial test
## 
## data:  6 and 7
## number of successes = 6, number of trials = 7, p-value = 0.3017
## alternative hypothesis: true probability of success is less than 0.95
## 95 percent confidence interval:
##  0.0000000 0.9926992
## sample estimates:
## probability of success 
##              0.8571429

As we can see both JL-2 and M51 have got statistically significant results i.e. H₀| P=0.95 is not rejected.

So we can apply null hypothesis for our formulas above in terms of H₀| P=0.95:

sum(dbinom(11:12,12,0.95))

## [1] 0.8816401

sum(dbinom(14:16,16,0.95))

## [1] 0.9570621

The results are more than satisfactory for both chinese and french submarines though french submarine has got \(96%\) chance compared to chinese \(88%\) for launching more than 90% of SLBMs. This is a real handicap.

Power of the test

To be honest 95% for successful launch of missile is not the goal of China and France at all. The goal is 99%. How much trials do both countries need to get this result? We use pwr.p.test function.

pwr.p.test(h = ES.h(5/6,0.99),power = 0.8,alternative = "less")

## 
##      proportion power calculation for binomial distribution (arcsine transformation) 
## 
##               h = -0.6407338
##               n = 15.0596
##       sig.level = 0.05
##           power = 0.8
##     alternative = less

pwr.p.test(h = ES.h(6/7,0.99),power = 0.8,alternative = "less")

## 
##      proportion power calculation for binomial distribution (arcsine transformation) 
## 
##               h = -0.5748585
##               n = 18.70881
##       sig.level = 0.05
##           power = 0.8
##     alternative = less

So we got 15 additional trials for JL-2 and 19 for M51. What is the power of test for JL-2 and M51?

binom.plot(n = 10,np = 500,conf.level = 0.95,type = "levelplot")

## Loading required package: lattice

binom.plot(c(6, 7), type = "xyplot")

pwr.p.test(h = ES.h(p1 = 5/6,p2 = 0.99),sig.level =0.95, n=6,alternative = "less")

## 
##      proportion power calculation for binomial distribution (arcsine transformation) 
## 
##               h = -0.6407338
##               n = 6
##       sig.level = 0.95
##           power = 0.9993462
##     alternative = less

pwr.p.test(h = ES.h(p1 = 6/7,p2 = 0.99),sig.level =0.95, n=7, alternative = "less")

## 
##      proportion power calculation for binomial distribution (arcsine transformation) 
## 
##               h = -0.5748585
##               n = 7
##       sig.level = 0.95
##           power = 0.9992267
##     alternative = less

As we can see test flight programs has got \(99.9%\) power so far. Significant results for both SLBMs.

Bayesian estimate

There is another way to consider this problem - Bayesian statistics (http://en.wikipedia.org/wiki/Bayesian_statistics). We can use bayes.binom.test and bayes.prop.test functions for comparison.

Current results

bsJL2<-bayes.binom.test(5,6,0.95)
summary(bsJL2)

##   Data
## number of successes = 5, number of trials = 6
## 
##   Model parameters and generated quantities
## theta: the relative frequency of success
## x_pred: predicted number of successes in a replication
## 
##   Measures
##         mean    sd HDIlo HDIup %<comp %>comp
## theta  0.751 0.143 0.481 0.991  0.954  0.046
## x_pred 4.508 1.303 2.000 6.000  0.004  0.996
## 
## 'HDIlo' and 'HDIup' are the limits of a 95% HDI credible interval.
## '%<comp' and '%>comp' are the probabilities of the respective parameter being
## smaller or larger than 0.95.
## 
##   Quantiles
##        q2.5%  q25% median  q75% q97.5%
## theta  0.427 0.662  0.771 0.862  0.964
## x_pred 2.000 4.000  5.000 6.000  6.000

plot(bsJL2)

bsM51<-bayes.binom.test(6,7,0.95)
summary(bsM51)

##   Data
## number of successes = 6, number of trials = 7
## 
##   Model parameters and generated quantities
## theta: the relative frequency of success
## x_pred: predicted number of successes in a replication
## 
##   Measures
##         mean    sd HDIlo HDIup %<comp %>comp
## theta  0.779 0.130 0.523 0.985  0.943  0.057
## x_pred 5.454 1.366 3.000 7.000  0.001  0.999
## 
## 'HDIlo' and 'HDIup' are the limits of a 95% HDI credible interval.
## '%<comp' and '%>comp' are the probabilities of the respective parameter being
## smaller or larger than 0.95.
## 
##   Quantiles
##        q2.5%  q25% median  q75% q97.5%
## theta  0.479 0.698  0.799 0.878  0.969
## x_pred 2.000 5.000  6.000 7.000  7.000

plot(bsM51)

suc <- c( 5, 6)
all <- c( 6, 7)

fit<-bayes.prop.test(suc,all)
plot(fit)

Future results

Suppose both JL-2 and M51 would pass 15 and 19 successful trials in the future. What would be the power and credible intervals for JL-2 and M51?

binom.plot(n = 30,np = 500,conf.level = 0.95,type = "levelplot")

binom.plot(c(21, 26), type = "xyplot")

bsJL2<-bayes.binom.test(5+15,6+15,0.95)
summary(bsJL2)

##   Data
## number of successes = 20, number of trials = 21
## 
##   Model parameters and generated quantities
## theta: the relative frequency of success
## x_pred: predicted number of successes in a replication
## 
##   Measures
##          mean    sd  HDIlo  HDIup %<comp %>comp
## theta   0.913 0.057  0.803  0.997  0.702  0.298
## x_pred 19.177 1.712 16.000 21.000  0.000  1.000
## 
## 'HDIlo' and 'HDIup' are the limits of a 95% HDI credible interval.
## '%<comp' and '%>comp' are the probabilities of the respective parameter being
## smaller or larger than 0.95.
## 
##   Quantiles
##         q2.5%   q25% median   q75% q97.5%
## theta   0.775  0.882  0.924  0.955  0.989
## x_pred 15.000 18.000 20.000 20.000 21.000

plot(bsJL2)

bsM51<-bayes.binom.test(6+19,7+19,0.95)
summary(bsM51)

##   Data
## number of successes = 25, number of trials = 26
## 
##   Model parameters and generated quantities
## theta: the relative frequency of success
## x_pred: predicted number of successes in a replication
## 
##   Measures
##          mean    sd  HDIlo  HDIup %<comp %>comp
## theta   0.929 0.047  0.837  0.998  0.608  0.392
## x_pred 24.148 1.769 21.000 26.000  0.000  1.000
## 
## 'HDIlo' and 'HDIup' are the limits of a 95% HDI credible interval.
## '%<comp' and '%>comp' are the probabilities of the respective parameter being
## smaller or larger than 0.95.
## 
##   Quantiles
##         q2.5%   q25% median   q75% q97.5%
## theta   0.812  0.903  0.938  0.964  0.991
## x_pred 20.000 23.000 25.000 25.000 26.000

plot(bsM51)

suc <- c( 5+15, 6+19)
all <- c( 6+15, 7+19)

fit<-bayes.prop.test(suc,all)
plot(fit)

JL2.all<-as.data.frame(bsJL2$mcmc_samples[1:5000,])
pp<-JL2.all$theta
dd5<-c()
for (i in 1:5000) {dd5[i]=sum(dbinom(11:12,12,pp[i]))}
hist(dd5,freq = T,xlim = c(0.1,1.0),breaks = 50,main="Histogram for distribution P(11<X<12) JL-2 SLBM launch", xlab = "Probability",col = "lightblue")

summary(dd5)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.007562 0.572300 0.771000 0.720100 0.905500 0.999900

M51.all<-as.data.frame(bsM51$mcmc_samples[1:5000,])
pp<-M51.all$theta
dd5<-c()
for (i in 1:5000) {dd5[i]=sum(dbinom(14:16,16,pp[i]))}
hist(dd5,freq = T,xlim = c(0.1,1.0),breaks = 50,main="Histogram for distribution P(14<X<16) M51 SLBM launch", xlab = "Probability",col = "lightblue")

summary(dd5)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.01794 0.80040 0.92750 0.86300 0.98210 1.00000

Conclusions

The current results of flight tests are satisfactory both for JL-2 and M51 in terms of power 99%.
M51 has got handicap in this race with JL-2 so far in terms of overall submarine readiness on patrol - 96% chance compared to chinese 88% for launching more than 90% of SLBMs.
In the future after 15 and 19 consecutive successful trials both SLBMs would have good chances to finish this race neck and neck.
Bayesian estimates give 95% credible intervals for JL-2 as \(0.481 < P_{JL2} < 0.991\) and M51 as \(0.523 < P_ {M51}< 0.985\) with the mean values \(P_{JL2} = 0.751\) and \(P_ {M51}=0.779\) so far for the current situation with trials.
Bayesian estimates give 95% credible intervals for JL-2 as \(0.803 < P_{JL2} < 0.997\) and M51 as \(0.837 < P_ {M51}< 0.998\) with the mean values \(P_{JL2} = 0.913\) and \(P_ {M51}=0.0.929\) after 15 and 19 consecutive successful trials.
The winner takes it all. Vive la France!

JL-2 and M51. The winner takes it all

Alexander Levakov, Senior Research Fellow, Ph.D

January, 2017