Problem

In the military science of ballistics, circular error probable (CEP) (also circular error probability or circle of equal probability) is a measure of a weapon system’s precision. It is defined as the radius of a circle, centered about the mean, whose boundary is expected to include the landing points of 50% of the rounds. See https://en.wikipedia.org/wiki/Circular_error_probable

We want to estimate M777 howitzer precision, using some open data: https://en.wikipedia.org/wiki/M777_howitzer. Testing at the Yuma Proving Ground by the US Army placed 13 of 14 Excalibur rounds, fired from up to 24 kilometers (15 mi), within 10 meters of their target,suggesting a circular error probable of about five meters. As of September 2015, nearly 770 Excalibur shells had been fired in combat.

library(binom)
library(pwr)
library(ggplot2)
library(BayesianFirstAid)
## Loading required package: rjags
## Loading required package: coda
## Linked to JAGS 4.1.0
## Loaded modules: basemod,bugs

Frequintist approach

We believe the true value for M777 Excalibur hit probability must be 99%. So our null hypothesis \(H0|P=0.99\), while alternative is \(H1|P<0.99\). We apply binom.test to verify it using frequentist approach with fixed parameters while data can vary.

binom.test(13,14,0.99)
## 
##  Exact binomial test
## 
## data:  13 and 14
## number of successes = 13, number of trials = 14, p-value = 0.1313
## alternative hypothesis: true probability of success is not equal to 0.99
## 95 percent confidence interval:
##  0.6613155 0.9981932
## sample estimates:
## probability of success 
##              0.9285714

So we can’t reject null hypothesis with p-value > 0.05: we have 13% probability of false rejecting null hypothesis when it is true. The 95% confidence interval for our estimate is 0.66 < P < 0.99. But is it really the case?

Bayesian approach

Now we apply Bayesian approach with fixed data while parameters can vary. Let’s do the estimate for credible interval using binom.bayes.

bs=binom.bayes(x = 13,n = 14,maxit = 10000)
print(bs)
##   method  x  n shape1 shape2 mean     lower     upper  sig
## 1  bayes 13 14   13.5    1.5  0.9 0.7523512 0.9997974 0.05
binom.bayes.densityplot(bs,500)

fit<-bayes.binom.test(x = 13,n = 14,p = 0.99,n.iter = 10000)
plot(fit)

summary(fit)
##   Data
## number of successes = 13, number of trials = 14
## 
##   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.875 0.080 0.718  0.995   0.99   0.01
## x_pred 12.255 1.629 9.000 14.000   0.00   1.00
## 
## '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.99.
## 
##   Quantiles
##        q2.5%   q25% median   q75% q97.5%
## theta  0.683  0.831   0.89  0.936  0.983
## x_pred 8.000 11.000  13.00 14.000 14.000
diagnostics(fit)
## 
## Iterations = 1:3334
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 3334 
## 
##   Diagnostic measures
##          mean    sd mcmc_se n_eff  Rhat
## theta   0.875 0.080   0.001  9944 1.001
## x_pred 12.255 1.629   0.017  9553 1.000
## 
## mcmc_se: the estimated standard error of the MCMC approximation of the mean.
## n_eff: a crude measure of effective MCMC sample size.
## Rhat: the potential scale reduction factor (at convergence, Rhat=1).
## 
##   Model parameters and generated quantities
## theta: The relative frequency of success
## x_pred: Predicted number of successes in a replication

We got almost the same result but more robust in sense of credibility - we got HDI (high density interval) for our estimate of model parameters given data: 0.718 < P < 0.996. More over we got the whole bunch of quantiles for model parameters: q2.5%=0.687, q25%=0.832, median=0.892, q75%=0.937, q97.5%=0.983. That’s it!

Power of test

Now we want to estimate the power of test.

pwr.p.test(ES.h(p1 = 0.8,p2 = 0.99),power=0.8,alternative = "less")
## 
##      proportion power calculation for binomial distribution (arcsine transformation) 
## 
##               h = -0.7269604
##               n = 11.69894
##       sig.level = 0.05
##           power = 0.8
##     alternative = less

As we can see the number of trials is significantly adequate.

Conclusions

  1. M777 Excalibur hit probability has the Bayesian estimate mean of 87.6% with the lower bound of 72% and upper bound of 99.4%.
  2. 14 trials proved to be adequate with the power of test of 80%.
  3. Bayesian estimate of binomial test gives us 95% credibility that out of 770 Excalibur shells been fired no less than 72% (554) had hit 10 meters circular area while 554/2=277 of fired shells had hit 5 meters circular area.
  4. Both approaches or paradigms (frequentist and Bayesian) could be applied for estimates of probability but one should take care making inferences on the results.