Demystifying M51
Alexander Levakov, Senior Research Fellow, Ph.D
November 14, 2018
M51 test fly
M51
The M51 SLBM is a submarine-launched ballistic missile, built by Airbus Defense & Space, 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://en.wikipedia.org/wiki/M51_(missile), https://missilethreat.csis.org/missile/m51/
SLBM data
load("slbm.dat")
library(DT)
datatable(slbm)library(ggplot2)
library(GGally)
#Here we use function from https://www.r-bloggers.com/multiple-regression-lines-in-ggpairs/
my_fn <- function(data, mapping, ...){
p <- ggplot(data = data, mapping = mapping) +
geom_point() +
geom_smooth(method=loess, fill="red", color="red", ...) +
geom_smooth(method=lm, fill="blue", color="blue", ...)
p
}
g = ggpairs(slbm,columns = 2:6, lower = list(continuous = my_fn))
g
Bayesian Linear Regression Model
library(rstanarm)## Loading required package: Rcpp## rstanarm (Version 2.17.3, packaged: 2018-02-17 05:11:16 UTC)## - Do not expect the default priors to remain the same in future rstanarm versions.## Thus, R scripts should specify priors explicitly, even if they are just the defaults.## - For execution on a local, multicore CPU with excess RAM we recommend calling## options(mc.cores = parallel::detectCores())## - Plotting theme set to bayesplot::theme_default().library(bayesplot)## This is bayesplot version 1.4.0## - Plotting theme set to bayesplot::theme_default()## - Online documentation at mc-stan.org/bayesplotoptions(mc.cores = parallel::detectCores())
fit.slbm.bs<-stan_glm(sqrt(R)~S+D+L+W+log(M)+log(P),
data=slbm,chains=4,iter=10000,seed=12345)
fit.slbm.bs## stan_glm
## family: gaussian [identity]
## formula: sqrt(R) ~ S + D + L + W + log(M) + log(P)
## observations: 26
## predictors: 7
## ------
## Median MAD_SD
## (Intercept) -40.6 106.7
## S 9.4 5.1
## D 37.4 19.8
## L 0.9 1.7
## W 3.3 5.6
## log(M) 6.5 15.0
## log(P) -7.3 6.0
## sigma 10.2 1.7
##
## Sample avg. posterior predictive distribution of y:
## Median MAD_SD
## mean_PPD 69.0 2.8
##
## ------
## For info on the priors used see help('prior_summary.stanreg').plot(fit.slbm.bs)
posterior_vs_prior(fit.slbm.bs)##
## Drawing from prior...
fit.slbm.bs2<-as.array(fit.slbm.bs)
mcmc_hist(fit.slbm.bs2)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
mcmc_trace(fit.slbm.bs2)
summary(fit.slbm.bs)##
## Model Info:
##
## function: stan_glm
## family: gaussian [identity]
## formula: sqrt(R) ~ S + D + L + W + log(M) + log(P)
## algorithm: sampling
## priors: see help('prior_summary')
## sample: 20000 (posterior sample size)
## observations: 26
## predictors: 7
##
## Estimates:
## mean sd 2.5% 25% 50% 75% 97.5%
## (Intercept) -40.0 109.1 -254.9 -112.3 -40.6 31.7 175.7
## S 9.4 5.4 -1.2 5.9 9.4 12.9 20.2
## D 37.4 20.3 -2.6 24.1 37.4 50.9 77.4
## L 0.9 1.8 -2.7 -0.3 0.9 2.0 4.4
## W 3.2 5.8 -8.5 -0.5 3.3 7.1 14.8
## log(M) 6.4 15.6 -24.3 -3.6 6.5 16.6 37.0
## log(P) -7.3 6.2 -19.5 -11.3 -7.3 -3.3 5.0
## sigma 10.4 1.8 7.6 9.2 10.2 11.4 14.6
## mean_PPD 69.0 2.9 63.2 67.1 69.0 70.9 74.7
## log-posterior -110.5 2.4 -116.3 -111.8 -110.1 -108.8 -107.1
##
## Diagnostics:
## mcse Rhat n_eff
## (Intercept) 1.2 1.0 8809
## S 0.0 1.0 12972
## D 0.2 1.0 9063
## L 0.0 1.0 12118
## W 0.1 1.0 11991
## log(M) 0.2 1.0 8316
## log(P) 0.1 1.0 13099
## sigma 0.0 1.0 10055
## mean_PPD 0.0 1.0 19069
## log-posterior 0.0 1.0 6234
##
## For each parameter, mcse is Monte Carlo standard error, 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).M51 Operational Range
M51 SLBM
slbm.pp <- posterior_predict(fit.slbm.bs,
newdata = data.frame(L=13,D=2.35,S=3,W=2,
M=53000,P=3000),seed=12345)
#Note. Payload is an expert estimate for the model
Range.km <- slbm.pp^2
Range.km <- Range.km[1:10000,]
quantile(Range.km,probs = c(0.1,0.5,0.9))## 10% 50% 90%
## 7968.799 11103.029 14849.838Range.km <- slbm.pp^2
Range.km <- Range.km[1:10000,]
ggplot(data=as.data.frame(Range.km), aes(Range.km)) +
geom_histogram(bins = 30,col="black",fill="green") +
geom_vline(xintercept = mean(Range.km), color = "red") +
geom_errorbarh(aes(y=-5, xmin=quantile(Range.km,0.1),
xmax=quantile(Range.km,0.9)),
data=as.data.frame(Range.km), col="#0094EA", size=3) +
ggtitle(label="Probability density of M51 range",
subtitle = "P(8000<Range.km<14800)=0.8")
M51 PEM flight test reliability
Here we use open source data on M51 flight test events and Bayesian Bernoulli model.
library(rstan)
M51_PEM_data <- c(TRUE,TRUE,TRUE,TRUE,FALSE,TRUE)
y <- as.integer(M51_PEM_data)
N <- length(y)
data_list <- list(N=N,y=y,a=2,b=1)
fit.ber <- stan(file = "bernouilly.stan",data = data_list,chains = 2,
iter = 5000,warmup = 1000,seed = 12345)
stan_hist(fit.ber,pars = "theta")
stan_trace(fit.ber)
RSTAN model
library(rstan)
print(fit.ber)## Inference for Stan model: bernouilly.
## 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
## theta 0.78 0.00 0.13 0.47 0.70 0.80 0.88 0.97 3169 1
## lp__ -5.32 0.02 0.78 -7.46 -5.49 -5.02 -4.83 -4.77 2036 1
##
## Samples were drawn using NUTS(diag_e) at Tue Dec 25 11:31:54 2018.
## 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).stan_model(file = "bernouilly.stan")## S4 class stanmodel 'bernouilly' coded as follows:
## data {
## int<lower=0> N;
## int<lower=0,upper=1> y[N];
## int a;
## int b;
## }
## parameters {
## real<lower=0,upper=1> theta;
## }
## model {
## theta ~ beta(a,b);
## for (n in 1:N)
## y[n] ~ bernoulli(theta);
## }
##
##
##
##
## ```
Conclusions
1.Applying Bayesian linear regression model for SLBM data (M,P,D,L,S,W) we can produce posterior distribution sample given M51 data (M=53000 kg,P=3000 kg,D=2.35 m,L=13 m,S=3,W=2) with 80% credible interval \(P(8000\le Range\le14800)=0.8\) and the mean \(Range=11000\).
2.M51 SLBM has operational Range corresponding to ICBM which could survive bolt-out-of-the blue attack and hit any possible target world wide.
3.Applying Bayesian Bernoulli model with \(beta(2,1)\) prior for the M51 flight test data we can produce posterior sample for the current M51 95% credible reliability interval as \(P(0.47\le \theta\le0.97)=0.95\) and the mean value \(\theta=0.80\). Compare these results with Trident II D5 DASO http://rpubs.com/alex-lev/377492.