RSM-56 Bulava

RSM-56 Bulava

RSM-56 Bulava

The RSM-56 Bulava (NATO reporting name SS-NX-30 or SS-N-32) is a submarine-launched ballistic missile (SLBM) developed for the Russian Navy and deployed in 2013 on the new Borei class of ballistic missile nuclear submarines. It is intended as the future cornerstone of Russia’s nuclear triad, and is the most expensive weapons project in the country.

A source in the Russian defense industry told TASS on June 29, 2018 that the D-30 missile system with the R-30 Bulava intercontinental ballistic missile had been accepted for service in the Russian Navy after its successful four-missile salvo launch tests in 2018.[9]

For more see: https://en.wikipedia.org/wiki/RSM-56_Bulava

BULAVA test events data

We can import BULAVA (RSM-56) test events from the open source WIKI data like this, applying substitution FAILURE=FALSE and SUCCESS=TRUE.

library(rvest)
## Loading required package: xml2
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(htmlTable)
library(DT)
url <- "https://en.wikipedia.org/wiki/RSM-56_Bulava"
bulava_test <- url %>%
  read_html()%>%
  html_nodes(xpath='//*[@id="mw-content-text"]/div/table[2]') %>%
  html_table()
bulava_test_data <- data.frame(cbind(bulava_test[[1]]$Date,bulava_test[[1]]$Result,bulava_test[[1]]$Position))

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

names(bulava_test_data)=c("Date","Result","Position")
bulava_test_data <- bulava_test_data %>% filter(Position=="Submerged") 
bulava_test_data %>% htmlTable
Date Result Position
1 23 September 2004 TRUE Submerged
2 21 December 2005 TRUE Submerged
3 7 September 2006 FALSE Submerged
4 25 October 2006 FALSE Submerged
5 29 June 2007 TRUE Submerged
6 18 September 2008 TRUE Submerged
7 28 November 2008 TRUE Submerged
8 23 December 2008 FALSE Submerged
9 15 July 2009 FALSE Submerged
10 9 December 2009 FALSE Submerged
11 7 October 2010 TRUE Submerged
12 29 October 2010 TRUE Submerged
13 27 June 2011 TRUE Submerged
14 27 August 2011 FALSE Submerged
15 28 October 2011 TRUE Submerged
16 23 December 2011 TRUE Submerged
17 6 September 2013 FALSE Submerged
18 9 September 2014 TRUE Submerged
19 29 October 2014 TRUE Submerged
20 28 November 2014 TRUE Submerged
21 15 November 2015 TRUE Submerged
22 27 September 2016 TRUE Submerged
23 26 June 2017 TRUE Submerged
24 22 May 2018 TRUE Submerged 
#datatable(bulava_test_data)

Reliability point estimate

We filtered BULAVA test data for Position=“Submerged” only.

summary(bulava_test_data$Result)
##    Mode   FALSE    TRUE 
## logical       7      17
NS <- sum(bulava_test_data$Result)
NALL <-length(bulava_test_data$Result) 
prop_success <- NS/NALL
prop_success
## [1] 0.7083333

As we see BULAVA has got 67% reliability (\(P=0.71\)) according to the test results to date. Can we make out something else applying current data? Let’s see!

Reliability confidence interval

Here we postulate our null hypothesis \(H_0:P=0.85\) i.e. we believe that BULAVA is 85% reliable in spite of test data. So we got alternative hypothesis \(H_1:P\ne 0.85\). Let’s prove it by binomial test.

binom.test(NS,NALL,0.85,alternative = "two.sided")
## 
##  Exact binomial test
## 
## data:  NS and NALL
## number of successes = 17, number of trials = 24, p-value = 0.07742
## alternative hypothesis: true probability of success is not equal to 0.85
## 95 percent confidence interval:
##  0.4890522 0.8738479
## sample estimates:
## probability of success 
##              0.7083333

As we see we can accept \(H_0:P=0.85\) i.e. we can’t reject \(H_0\). That’s good news. Actually we have got 95% confidence interval for BULAVA reliability \(0.489\le P\le0.874\). What is the true reliability? To this end we must estimate SLBM reliability throughout the time i.e. from test event to test event. Let’s see!

Proportion of success model

We apply BULAVA test data for the proportion of success as function of the current 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 = bulava_test_data$Result,prior_prop = c(1,1),gr_name = "BULAVA test events")

Gaussian prior Beta(2,2)

prop_model(data = bulava_test_data$Result,prior_prop = c(2,2),gr_name = "BULAVA test events")

As wee see our belief in BULAVA reliability is changing in time from event to event. That’s great! But is there any more possibility to get more robust estimate of SLBM reliability applying current data? Let’ see!

Bayesian STAN beta-bernouilli model

Here we apply Bayesian inference. For more see relevant math stuff:
1.https://en.wikipedia.org/wiki/Bayesian_inference
2.https://en.wikipedia.org/wiki/Beta_distribution
3.https://en.wikipedia.org/wiki/Bernoulli_distribution

library(rstan)
## Loading required package: StanHeaders
## rstan (Version 2.17.3, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## 
## Attaching package: 'rstan'
## The following object is masked from 'package:tidyr':
## 
##     extract
library(BH)
library(bayesplot)
## This is bayesplot version 1.4.0
## - Plotting theme set to bayesplot::theme_default()
## - Online documentation at mc-stan.org/bayesplot
library(mcmc)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

stan_code  <-  '
data { 
  int<lower=0> N; 
  int<lower=0,upper=1> y[N];
  
  real a;
  real b;
} 

transformed data{
  int<lower=0,upper=1> yy[N];
  for (k in 1:N)
  yy[k]=abs(1-y[k]);
}


parameters {
  real<lower=0,upper=1> theta[N];
} 
model {
  for (j in 1:N){
  theta[j] ~ beta(a+sum(y[1:j]),b+sum(yy[1:j]));
      for (n in 1:j) //
      y[j] ~ bernoulli(theta[j]);
  }
}

generated quantities {
  real<lower=0,upper=1> pp[N]; 
  for (j in 1:N)
  pp[j]=beta_rng(a+sum(y[1:j]),b+sum(yy[1:j]));
} 
'

Posterior sample with Beta(2,2)

y <- as.integer(bulava_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,
            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[24]    0.61       0 0.09 0.42 0.54 0.61 0.67  0.78  8000    1
## theta[24] 0.79       0 0.06 0.67 0.75 0.79 0.83  0.89  8000    1
## 
## Samples were drawn using NUTS(diag_e) at Sun Jan 27 14:52:49 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 BULAVA reliability.

Conclusions

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