French M51 SLBM and Triomphant SSBN

French M51 SLBM and Triomphant SSBN

Note:It is the updated version of the previous report on M51 (see https://rpubs.com/alex-lev/439745) as far as reliability issue concerned.

Case story

The following abstract is a quote: The French Triomphant class nuclear ballistic missile submarine (SSBN) Le Téméraire test-fired an M51 submarine-launched ballistic missile in the Atlantic off of Finistère, France in the early hours of June 12th, 2020. Some sort of a test appeared to be in the works just three days ago when Le Téméraire was spotted sailing out of port with huge test instrumentation masts attached that are commonly fitted to submarines prior to developmental ballistic missile launches. Then, last evening, our good friend @aircraftspots began tracking a U.S. Air Force RC-135S Cobra Ball ballistic missile and rocket tracking aircraft flying out over the Caribbean. Not long after, a French Falcon 50 maritime patrol aircraft showed up in the area, indicating a launch was likely imminent.

For more see: https://www.thedrive.com/the-war-zone/34043/france-just-test-fired-a-submarine-launched-ballistic-missile-in-the-atlantic

About M51

The M51 SLBM is a submarine-launched ballistic missile, built by EADS Astrium Space Transportation, and deployed with the French Navy. Designed to replace the M45 SLBM (In French terminology the MSBS - Mer-Sol-Balistique-Stratégique “Sea-ground-Strategic ballistic”), it was first deployed in 2010.

For more see:

https://military.wikia.org/wiki/M51_(missile)

https://en.wikipedia.org/wiki/M51_(missile)

https://www.naval-group.com/en/news/dcns-starts-the-last-adaptation-programme-on-the-french-ssbns-for-the-m51-missile/

Flight test events data

M51_test_data <- as.logical(c("TRUE","TRUE","TRUE","TRUE","FALSE","TRUE","TRUE"))

Reliability point estimate

summary(M51_test_data)
##    Mode   FALSE    TRUE 
## logical       1       6
NS <- sum(M51_test_data)
NALL <-length(M51_test_data) 
prop_success <- NS/NALL
prop_success
## [1] 0.8571429

As we see M51 has got 86% reliability (\(P=0.86\)) 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 M51 is 85% reliable in spite of test data. So we got alternative hypothesis \(H_1:P\neq 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 = 6, number of trials = 7, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.85
## 95 percent confidence interval:
##  0.4212768 0.9963897
## sample estimates:
## probability of success 
##              0.8571429

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 M51 reliability \(0.489\le P \le 0.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 M51 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 = M51_test_data,prior_prop = c(1,1),gr_name = "M51 test events")
## ── Attaching packages ───────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.2.1     ✔ purrr   0.3.3
## ✔ tibble  2.1.3     ✔ dplyr   0.8.5
## ✔ tidyr   1.0.0     ✔ stringr 1.4.0
## ✔ readr   1.3.1     ✔ forcats 0.4.0
## ── Conflicts ──────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## Warning: `data_frame()` is deprecated, use `tibble()`.
## This warning is displayed once per session.

Gaussian prior Beta(2,2)

prop_model(data = M51_test_data,prior_prop = c(2,2),gr_name = "M51 test events")

Bayesian STAN beta-Bernoulli 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.19.2, 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.7.0
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
##    * Does _not_ affect other ggplot2 plots
##    * See ?bayesplot_theme_set for details on theme setting
library(mcmc)


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)

knitr::opts_chunk$set(echo = TRUE)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
y <- as.integer(M51_test_data)
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)
## 'pars' not specified. Showing first 10 parameters by default.

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[7]    0.73       0 0.13 0.44 0.64 0.74 0.82  0.93  7830    1
## theta[7] 0.83       0 0.08 0.64 0.78 0.85 0.90  0.96 13156    1
## 
## Samples were drawn using NUTS(diag_e) at Sun Jun 14 16:19:57 2020.
## 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 M51 reliability.

Conclusions

1.The point estimate of M51 reliability looks like \(P=0.857\) in accordance with the data we have got so far.

2.The 95% confidence interval for M51 reliability looks like \(0.421 \le P \le0.996\) with the mean value \(P=0.857\).

3.Proportion model makes our belief in M51 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 M51 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 M51 reliability looks like 95% credible interval \(0.64 \le P \le 0.96\) with the mean value \(P=0.85\).

6.All the estimates made above are based on the assumption that M51 SLBM flight test data are correct in terms of Success and Failure events.

7.The most interesting and curious fact about M51 is that France deployed SLBM without declaring preliminary test flight program as it was with Trident II D5 (see https://rpubs.com/alex-lev/377492). The ongoing flight test events (seven during the past ten years with one failed) with some sort of modernization from one event to another make M51’s reliability some kind of mystery, when nothing is clear about final and stable version of SLBM as well as of nuclear posture that France persues boldly along with USA and UK.