Standard Missile - 3 (SM-3)

Standard Missile - 3 (SM-3)

Aegis Ballistic Missile Defense System

The Aegis Ballistic Missile Defense System (Aegis BMD or ABMD) is a United States Department of Defense Missile Defense Agency program developed to provide missile defense against short to intermediate-range ballistic missiles. It is part of the United States national missile defense strategy. Aegis BMD (also known as Sea-Based Midcourse) is designed to intercept ballistic missiles post-boost phase and prior to reentry.

For more see:

https://en.wikipedia.org/wiki/Aegis_Ballistic_Missile_Defense_System

https://youtu.be/X1ACKtE62n4

AEGIS test events data

We can import AEGIS 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/Aegis_Ballistic_Missile_Defense_System"
aegis_test <- url %>%
  read_html()%>%
  html_nodes(xpath='//*[@id="mw-content-text"]/div/table[1]') %>%
  html_table()
aegis_test_data <- data.frame(cbind(aegis_test[[1]]$Date,aegis_test[[1]]$Result))

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

names(aegis_test_data)=c("Date","Result")
datatable(aegis_test_data)
#aegis_test_data %>% htmlTable

Reliability point estimate

summary(aegis_test_data$Result)
##    Mode   FALSE    TRUE 
## logical      10      42
NS <- sum(aegis_test_data$Result)
NALL <-length(aegis_test_data$Result) 
prop_success <- NS/NALL
prop_success
## [1] 0.8076923

As we see AEGIS has got 81% reliability according to the test results to date. Can we make out somethig 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 = 42, number of trials = 52, p-value = 0.4345
## alternative hypothesis: true probability of success is not equal to 0.85
## 95 percent confidence interval:
##  0.6746587 0.9037324
## sample estimates:
## probability of success 
##              0.8076923

Proportion of success model

We apply AEGIS test data for the proportion of success as function of the current test event using beta prior as Bayesian likelyhood.

#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 = aegis_test_data$Result,prior_prop = c(1,1),gr_name = "AEGIS test events")

Gaussian prior Beta(2,2)

prop_model(data = aegis_test_data$Result,prior_prop = c(2,2),gr_name = "AEGIS 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(aegis_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[52]    0.79       0 0.05 0.67 0.75 0.79 0.82  0.88  7671    1
## theta[52] 0.89       0 0.03 0.82 0.87 0.89 0.91  0.94  8000    1
## 
## Samples were drawn using NUTS(diag_e) at Sun Jan 27 14:48:57 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 conclusions about AEGIS reliability.

Conclusions

  1. The point estimate of AEGIS reliability looks like \(P=0.807\) in accordance with the data we have so far.
  2. The 95% confidence interval for AEGIS reliability looks like \(0.675\le P\le0.904\) with the mean value \(0.807\).
  3. Proportion model makes our belief in AEGIS 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 AEGIS 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 AEGIS reliability looks like 95% credible interval \(0.82\le P\le0.94\) with the mean value \(0.89\).
  6. All the estimates made above are based on the assumption that AEGIS flight test data are correct in terms of \(Success\) and \(Failure\) events.