Bayesian Linear Regression for Heavy Artillery. K9 Thunder case
Alexander Levakov, Senior Research Fellow, Ph.D
September, 2020
K9 Thunder
1 About K9 Thunder
The K9 Thunder is a South Korean self-propelled 155 mm howitzer designed and developed by the Agency for Defense Development and Samsung Aerospace Industries for the Republic of Korea Armed Forces, and is now manufactured by Hanwha Defense. K9 howitzers operate in groups with the K10 automatic ammunition resupply vehicle variant. The entire K9 fleet operated by the ROK Armed Forces is now undergoing upgrades to K9A1 standard, and a further development of a K9A2 variant is in process.
For more see: https://en.wikipedia.org/wiki/K9_Thunder
https://www.army-technology.com/projects/thunderselfpropelled/
2 Data
All data for long range guns (LRG) is here: https://en.wikipedia.org/wiki/List_of_railway_artillery
Some interesting statistics is here:http://rpubs.com/alex-lev/357849
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, unionlibrary(DT)
lrg <- read.delim("lrg.txt")
datatable(lrg)#lrg3 Bayesian Linear Regression Model
Here we fit Bayesian linear regression model for data with Range_M < 70000. For more about Bayesian linear model see: https://rpubs.com/alex-lev/360401
library(readxl)
library(ggplot2)
library(tidyverse)## ── Attaching packages ───────────────────────────────────────────────── tidyverse 1.3.0 ──## ✔ tibble 3.0.3 ✔ purrr 0.3.3
## ✔ tidyr 1.1.1 ✔ 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()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 settinglibrary(rstanarm)## Loading required package: Rcpp## Registered S3 method overwritten by 'xts':
## method from
## as.zoo.xts zoo## rstanarm (Version 2.18.2, packaged: 2018-11-08 22:19:38 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(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':
##
## extractoptions(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
lrg2 <- lrg%>%filter(Range_M<70000)
fit.lrg.2.Bayes <- rstanarm::stan_glm(Range_M~Velocity_MPS+Barrel_M+Elevation,data = na.omit(lrg2),iter=5000,chain=4,seed=12345)
# Verification of the model
summary(fit.lrg.2.Bayes)##
## Model Info:
##
## function: stan_glm
## family: gaussian [identity]
## formula: Range_M ~ Velocity_MPS + Barrel_M + Elevation
## algorithm: sampling
## priors: see help('prior_summary')
## sample: 10000 (posterior sample size)
## observations: 24
## predictors: 4
##
## Estimates:
## mean sd 2.5% 25% 50% 75% 97.5%
## (Intercept) -44164.5 7280.1 -58485.9 -48901.6 -44161.3 -39423.9 -29574.6
## Velocity_MPS 59.1 7.9 43.3 53.9 59.1 64.2 74.5
## Barrel_M 736.2 196.3 345.7 610.0 733.8 864.5 1123.6
## Elevation 409.8 124.3 165.2 327.7 411.6 490.4 652.2
## sigma 5341.9 889.8 3922.5 4715.5 5232.0 5840.1 7427.3
## mean_PPD 31666.9 1561.7 28606.6 30647.8 31655.3 32672.8 34784.7
## log-posterior -256.4 1.7 -260.8 -257.3 -256.0 -255.1 -254.1
##
## Diagnostics:
## mcse Rhat n_eff
## (Intercept) 69.0 1.0 11125
## Velocity_MPS 0.1 1.0 7895
## Barrel_M 2.3 1.0 7388
## Elevation 1.4 1.0 8001
## sigma 11.3 1.0 6191
## mean_PPD 15.4 1.0 10317
## log-posterior 0.0 1.0 3741
##
## 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).plot(fit.lrg.2.Bayes)
stan_hist(fit.lrg.2.Bayes)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
stan_trace(fit.lrg.2.Bayes)
draws <- as.data.frame(fit.lrg.2.Bayes)
bayesplot::mcmc_pairs(draws)
hist(bayes_R2(fit.lrg.2.Bayes),col="green",main = "R2 distribution of the model",xlab = "R2",xlim = c(0.75,0.95),breaks = 30)
fit.lrg.2.Bayes %>%
rhat() %>%
mcmc_rhat() +
yaxis_text()
4 Predicting K9 Thunder Range from posterior sample
Here we apply open source data on K9 Thunder.
pp <- posterior_predict(fit.lrg.2.Bayes,newdata = data.frame(Velocity_MPS=928,Barrel_M=155*52/1000,Elevation=60,seed=12345))
qq <- quantile(pp,probs = c(0.025,0.5,0.975))
qq## 2.5% 50% 97.5%
## 29530.31 41153.21 53089.29\[P(29179.38<Range<53156.71)=0.95\]
hist(pp,breaks = 50,col = "red",main = "Posterior destribution for K9 Thunder Range",xlab = "Range")
abline(v = qq[2])
qqq <- as.vector(qq)
lines(c(qqq[1],qqq[3]),c(-10,-10),col="blue",lwd=5)
5 Conclusion
- Bayesian linear regression model
Range_M~Velocity_MPS+Barrel_M+Elevationwith normal priors is relevant with all variables in the data set for historical long range guns data (\(Range < 70 km\)). - Applying Bayesian linear model on real heavy artillery gun K9 Thunder we have got nice result for the range of fire relevant results in terms of 95% credible interval for posterior sample as \(P(29179.38<Range<53156.71)=0.95\).