Demystifying SS-N-32
Alexander Levakov, Senior Research Fellow, Ph.D
December, 2018
SS-N-32
SS-N-32
The SS-N-32 “Bulava” is an intercontinental-range, submarine-launched, solid propellant ballistic missile. Alongside the SS-25 and the SS-27, both land-based ICBMs, the Bulava represents a core component of Russia’s future strategic nuclear force. Development of the program began in the 1990’s with official production contacts going into effect in the 2007-2008 timeframe. The Bulava was designed to be deployed onto Russia’s Borey-class ballistic missile submarines (SSBN’s), also referred to as Dolgorukiy-class or Project 955 submarines, which are capable of holding 12-16 missiles each.
For more see: https://missilethreat.csis.org/missile/ss-n-32-bulava/
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
We apply SLBM data frame without SS-N-32 data for the purpose of model accuracy i.e to be free of noise in our data.
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[-16,],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: 25
## predictors: 7
## ------
## Median MAD_SD
## (Intercept) -40.0 110.2
## S 9.4 5.5
## D 37.1 20.3
## L 0.9 1.8
## W 3.2 5.7
## log(M) 6.5 15.7
## log(P) -7.4 6.1
## sigma 10.5 1.8
##
## Sample avg. posterior predictive distribution of y:
## Median MAD_SD
## mean_PPD 67.9 3.0
##
## ------
## 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: 25
## predictors: 7
##
## Estimates:
## mean sd 2.5% 25% 50% 75% 97.5%
## (Intercept) -40.0 114.0 -265.1 -114.8 -40.0 33.7 189.6
## S 9.4 5.6 -1.6 5.7 9.4 13.1 20.6
## D 37.3 21.2 -4.7 23.6 37.1 50.9 79.8
## L 0.9 1.9 -2.8 -0.3 0.9 2.1 4.5
## W 3.1 6.0 -8.7 -0.7 3.2 7.0 15.1
## log(M) 6.5 16.3 -26.7 -4.1 6.5 17.1 38.7
## log(P) -7.4 6.5 -20.3 -11.5 -7.4 -3.3 5.5
## sigma 10.8 1.9 7.8 9.4 10.5 11.8 15.2
## mean_PPD 68.0 3.1 61.8 65.9 67.9 70.0 74.1
## log-posterior -107.5 2.4 -113.3 -108.8 -107.1 -105.7 -104.0
##
## Diagnostics:
## mcse Rhat n_eff
## (Intercept) 1.2 1.0 9010
## S 0.0 1.0 14324
## D 0.2 1.0 9847
## L 0.0 1.0 12107
## W 0.1 1.0 14134
## log(M) 0.2 1.0 8580
## log(P) 0.1 1.0 14045
## sigma 0.0 1.0 10090
## mean_PPD 0.0 1.0 18807
## log-posterior 0.0 1.0 5926
##
## 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).SS-N-32 Operational Range
Here we use open source data (see above) for the Bayesian linear regression model.
slbm.pp <- posterior_predict(fit.slbm.bs,
newdata = data.frame(L=12.1,D=2,S=3,W=2,
M=36800,P=1150),seed=12345)
Range.km <- slbm.pp^2
Range.km <- Range.km[1:10000,]
quantile(Range.km,probs = c(0.05,0.5,0.95))## 5% 50% 95%
## 5751.587 9275.903 13441.715Range.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 SS-N-32 range")
SS-N-32 PEM flight test reliability
Here we use open source data on SS-N-32 flight test events and Bayesian Bernoulli model. See http://rpubs.com/alex-lev/453699
library(rstan)
SSN32_PEM_data <- c(FALSE,TRUE,TRUE,TRUE,FALSE,FALSE,FALSE,TRUE,TRUE,TRUE,
FALSE,FALSE,FALSE,TRUE,TRUE,TRUE,FALSE,
TRUE,TRUE,FALSE,TRUE,TRUE,TRUE,FALSE,FALSE,
TRUE,FALSE,TRUE,TRUE,TRUE,TRUE,TRUE)
y <- as.integer(SSN32_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)## 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)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.63 0.00 0.08 0.46 0.58 0.63 0.69 0.78 3039 1
## lp__ -23.60 0.01 0.70 -25.57 -23.79 -23.32 -23.14 -23.09 3766 1
##
## Samples were drawn using NUTS(diag_e) at Sun Dec 23 17:25:47 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 SS-N-32 data (M=36800 kg,P=1150 kg,D=2 m,L=12.1 m,S=3,W=2) with 90% credible interval \(P(5752\le Range\le13441)=0.9\) and the mean \(Range=9276\).
2.SS-N-32 has operational Range corresponding to SLBM which could hit any possible target world wide.
3.Applying Bayesian Bernoulli model with \(beta(2,1)\) prior for the SS-N-32 flight test data we can produce posterior sample for the current SS-N-32 95% credible reliability interval as \(P(0.46\le \theta\le0.78)=0.95\) and the mean value \(\theta=0.63\). Compare these results with Trident II D5 DASO http://rpubs.com/alex-lev/377492.