Ground-Based Interceptor

Ground-Based Interceptor

Ground-Based Midcourse Defense

Ground-Based Midcourse Defense (GMD) is the United States’ anti-ballistic missile system for intercepting incoming warheads in space, during the midcourse phase of ballistic trajectory flight. It is a major component of the American missile defense strategy to counter ballistic missiles, including intercontinental ballistic missiles (ICBMs) carrying nuclear, chemical, biological or conventional warheads. The system is deployed in military bases in the states of Alaska and California.

For more see:

https://en.wikipedia.org/wiki/Ground-Based_Midcourse_Defense

https://missilethreat.csis.org/system/gmd/

GBMD intercept test events data

We can import GBMD intercept test events from the open source WIKI data like this, applying substitution FAILURE=FALSE and SUCCESS=TRUE.

library(rvest)
## Loading required package: xml2
library(htmlTable)
library(DT)
url <- "https://en.wikipedia.org/wiki/Ground-Based_Midcourse_Defense"
gbmd_test <- url %>%
  read_html()%>%
  html_nodes(xpath='//*[@id="mw-content-text"]/div/table[1]') %>%
  html_table()
gbmd_test_data <- data.frame(cbind(gbmd_test[[1]]$Date,gbmd_test[[1]]$Result))

gbmd_test_data[,2] <- ifelse(gbmd_test[[1]]$Result=="Failure",FALSE,TRUE)

names(gbmd_test_data)=c("Date","Result")
datatable(gbmd_test_data,caption = "GBMD intercept events")
#gbmd_test_data %>% htmlTable

Reliability point estimate

summary(gbmd_test_data$Result)
##    Mode   FALSE    TRUE 
## logical       8      12
NS <- sum(gbmd_test_data$Result)
NALL <-length(gbmd_test_data$Result) 
prop_success <- NS/NALL
prop_success
## [1] 0.6

As we see GBMD has got 60% reliability according to the test results to date. Can we make out something else applying current data? Let’s see!

Reliability confidence interval

binom.test(NS,NALL,0.85)
## 
##  Exact binomial test
## 
## data:  NS and NALL
## number of successes = 12, number of trials = 20, p-value =
## 0.005921
## alternative hypothesis: true probability of success is not equal to 0.85
## 95 percent confidence interval:
##  0.3605426 0.8088099
## sample estimates:
## probability of success 
##                    0.6

Proportion of success model

We apply GBMD intercept test data for the proportion of success as function of the current intercept test event using beta prior as Bayesian likelihood.

#The prop_model function - Rasmus Bååth R code

# This function takes a number of successes and failuers coded as a TRUE/FALSE
# or 0/1 vector. This should be given as the data argument.
# The result is a visualization of the how a Beta-Binomial
# model gradualy learns the underlying proportion of successes 
# using this data. The function also returns a sample from the
# posterior distribution that can be further manipulated and inspected.
# The default prior is a Beta(1,1) distribution, but this can be set using the
# prior_prop argument.

# Make sure the packages tidyverse and ggridges are installed, otherwise run:
# install.packages(c("tidyverse", "ggridges"))

# Example usage:
# data <- c(TRUE, FALSE, TRUE, TRUE, FALSE, TRUE, TRUE)
# prop_model(data)
prop_model <- function(data = c(), prior_prop = c(1, 1), n_draws = 10000,
                       gr_name="Proportion graph") {
  library(tidyverse)
 
  data <- as.logical(data)
  # data_indices decides what densities to plot between the prior and the posterior
  # For 20 datapoints and less we're plotting all of them.
  data_indices <- round(seq(0, length(data), length.out = min(length(data) + 1, 40)))
  
  # dens_curves will be a data frame with the x & y coordinates for the 
  # denities to plot where x = proportion_success and y = probability
  proportion_success <- c(0, seq(0, 1, length.out = 100), 1)
  dens_curves <- map_dfr(data_indices, function(i) {
    value <- ifelse(i == 0, "Prior", ifelse(data[i], "Success", "Failure"))
    label <- paste0("n=", i)
    probability <- dbeta(proportion_success,
                         prior_prop[1] + sum(data[seq_len(i)]),
                         prior_prop[2] + sum(!data[seq_len(i)]))
    probability <- probability / max(probability)
    data_frame(value, label, proportion_success, probability)
  })
  # Turning label and value into factors with the right ordering for the plot
  dens_curves$label <- fct_rev(factor(dens_curves$label, levels =  paste0("n=", data_indices )))
  dens_curves$value <- factor(dens_curves$value, levels = c("Prior", "Success", "Failure"))
  
  graph_label <- paste("Prior likelihood distribution Beta(a =", 
                       as.character(prior_prop[1]),", b =",
                                    as.character(prior_prop[2]),")") 
  
  p <- ggplot(dens_curves, aes(x = proportion_success, y = label,
                               height = probability, fill = value)) +
    ggridges::geom_density_ridges(stat="identity", color = "white", alpha = 0.8,
                                  panel_scaling = TRUE, size = 1) +
    scale_y_discrete("", expand = c(0.01, 0)) +
    scale_x_continuous("Proportion of success") +
    scale_fill_manual(values = hcl(120 * 2:0 + 15, 100, 65), name = "", drop = FALSE,
                      labels =  c("Prior   ", "Success   ", "Failure   ")) +
    ggtitle(paste0(gr_name, ": ", sum(data),  " successes, ", sum(!data), " failures"),
            subtitle = graph_label) +
    labs(caption = "based on Rasmus Bååth R code") +
    theme_light() +
    theme(legend.position = "top")
  print(p)
  
  # Returning a sample from the posterior distribution that can be further 
  # manipulated and inspected
  posterior_sample <- rbeta(n_draws, prior_prop[1] + sum(data), prior_prop[2] + sum(!data))
  invisible(posterior_sample)
}

Uniform prior Beta(1,1)

prop_model(data = gbmd_test_data$Result,prior_prop = c(1,1),gr_name = "GBMD intercept test events")

Gaussian prior Beta(2,2)

prop_model(data = gbmd_test_data$Result,prior_prop = c(2,2),gr_name = "GBMD intercept test events")

Bayesian STAN beta-bernouilli model

See relevant math stuff here:

  1. https://en.wikipedia.org/wiki/Bayesian_inference
  2. https://en.wikipedia.org/wiki/Beta_distribution
  3. https://en.wikipedia.org/wiki/Bernoulli_distribution

Posterior sample with Beta(2,2)

y <- as.integer(gbmd_test_data$Result)
N <- length(y)
data_list <- list(N=N,y=y,a=2,b=2)
fit.ber <- stan(model_code = stan_code,data = data_list,chains = 2,verbose = FALSE,
            iter = 5000,warmup = 1000,seed=12345)

rstan::traceplot(fit.ber)

s1 <- as.array(fit.ber)
pp_par <- paste("pp[",1:NALL,"]",sep = "")
pp_par2 <- paste("theta[",1:NALL,"]",sep = "")
bayesplot::mcmc_intervals(x = s1,pars = pp_par)

bayesplot::mcmc_intervals(x = s1,pars = pp_par2)

Reliability crystall ball magic

last_par_pp<- paste("pp[",NALL,"]",sep = "")
last_par_theta<- paste("theta[",NALL,"]",sep = "")
bayesplot::mcmc_hist(x = s1,pars = c(last_par_pp,last_par_theta))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

plot(fit.ber,pars = c(last_par_pp,last_par_theta))
## ci_level: 0.8 (80% intervals)
## outer_level: 0.95 (95% intervals)

print(fit.ber,pars = c(last_par_pp,last_par_theta))
## Inference for Stan model: 8110aaf154812cced3e43a9d64740f79.
## 2 chains, each with iter=5000; warmup=1000; thin=1; 
## post-warmup draws per chain=4000, total post-warmup draws=8000.
## 
##           mean se_mean   sd 2.5%  25%  50%  75% 97.5% n_eff Rhat
## pp[20]    0.59       0 0.10 0.38 0.52 0.59 0.66  0.77  8000    1
## theta[20] 0.77       0 0.06 0.64 0.73 0.78 0.82  0.88  8000    1
## 
## Samples were drawn using NUTS(diag_e) at Wed Mar 27 08:42:51 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

Here is a moment of truth - we see significant difference between two distributions of the posterior samples: \(theta\), produced by the data through Bernoulli experiment with the prior Beta(a,b) distribution as Bayesian likelihood and \(pp\), produced by Beta(a,b) random distribution only. The last one is similar to the proportion model above. Both samples use Gaussian like \(Beta(2,2)\) distribution. Now we are ready to make the final conclusion about GBMD reliability.

Some GBMD missiles per target

Here we use previous posterior sample for estimating probability of ICBM intercept by two, three and four GBMD missiles launched simultaneously.

prob_hist <- function(prob,num){
  hist_label <- "Probability density of intercept by"
  ggplot(data=as.data.frame(prob), aes(prob))+
  geom_histogram(bins = 30,col="black",fill="green") + 
  geom_vline(xintercept = mean(prob), color = "red") + 
  geom_errorbarh(aes(y=-5, xmin=quantile(prob,0.1),
                     xmax=quantile(prob,0.9)), 
                 data=as.data.frame(prob), col="#0094EA", size=3) +
  ggtitle(paste(hist_label,num," GBMD missiles",sep=" ")) + labs(x="Probability")
}

s2 <- as.data.frame(fit.ber)
prob_gbmd <- s2$`theta[20]`
#Two missiles per target
prob_gbmd_2 <- 1-(1-prob_gbmd)^2
quantile(prob_gbmd_2,c(0.025,0.5,0.975))
##      2.5%       50%     97.5% 
## 0.8685942 0.9498334 0.9856257
prob_hist(prob_gbmd_2,"2")

#Three missiles per target
prob_gbmd_3 <- 1-(1-prob_gbmd)^3
quantile(prob_gbmd_3,c(0.025,0.5,0.975))
##      2.5%       50%     97.5% 
## 0.9523655 0.9887637 0.9982766
prob_hist(prob_gbmd_3,"3")

#Four missiles per target
prob_gbmd_4 <- 1-(1-prob_gbmd)^4
quantile(prob_gbmd_4,c(0.025,0.5,0.975))
##      2.5%       50%     97.5% 
## 0.9827325 0.9974833 0.9997934
prob_hist(prob_gbmd_4,"4")

Conclusion

  1. The point estimate of current GBMD intercept reliability looks like \(P=0.6\) in accordance with the data we have so far.
  2. The 95% confidence interval for GBMD intercept reliability looks like \(0.36\le P\le0.81\).
  3. Proportion model makes our belief in GBMD intercept reliability stronger and stronger from one event to another no matter what prior we use: \(Beta(1,1)\) or \(Beta(2,2)\).
  4. Bayesian estimate of the posterior sample for the GBMD intercept reliability \(theta\) produced by Bernoulli experiment with data and \(Beta(2,2)\) prior as likelihood is significantly greater than \(Beta(2,2)\) random distribution for the proportion of success.
  5. We believe that true GBMD intercept reliability for single missile looks like 95% credible interval \(0.64\le P_{1}\le0.88\) with the mean value \(M[P_{1}]=0.78\).
  6. GBMD intercept reliability by two missiles looks like 95% credible interval \(0.868\le P_{2}\le0.986\) with the mean value \(M[P_{2}]=0.95\).
  7. GBMD intercept reliability by three missiles looks like 95% credible interval \(0.952\le P_{3}\le0.998\) with the mean value \(M[P_{3}]=0.989\).
  8. GBMD intercept reliability by four missiles looks like 95% credible interval \(0.983\le P_{4}\le0.999\) with the mean value \(M[P_{4}]=0.997\).
  9. All the estimates made above are based on the assumption that GBMD intercept test data are correct in terms of \(Success\) and \(Failure\) events.