Stan: Bayesian Modeling Examples

Sir Panda

Stan: Bayesian Modeling Examples

statistics

Sample scripts and examples for Bayesian modeling using Stan.

Author

Sir Panda

Published

January 1, 2018

Keywords

Stan, Bayesian statistics, statistical modeling, MCMC, probabilistic programming, sensitivity analysis, uncertainty analysis, robust analysis

Please don’t mind this post. I use this to try out various highlighting styles for my code and formats.

Algorithm: Bootstrap Confidence Interval

Input: Data X = \{x_1, x_2, \ldots, x_n\}, statistic \theta(X), confidence level \alpha

Output: Confidence interval [\theta_L, \theta_U]

For b = 1, \ldots, B:

Draw bootstrap sample X_b^* by sampling n observations from X with replacement.
Compute \theta_b^* = \theta(X_b^*).

Sort \{\theta_1^*, \theta_2^*, \ldots, \theta_B^*\} in ascending order.
Set \theta_L = \text{quantile}(\theta^*, \alpha/2).
Set \theta_U = \text{quantile}(\theta^*, 1-\alpha/2).
Return [\theta_L, \theta_U].

\begin{algorithm} \caption{Metropolis–Hastings} \begin{algorithmic} \STATE initialize $\theta_0$ \FOR{$t = 1$ to $T$} \STATE propose $\theta^* \sim q(\cdot \mid \theta_{t-1})$ \ENDFOR \end{algorithmic} \end{algorithm}

Y^{\operatorname{Post}} = \beta_{0} + \beta_{1}^{\operatorname{Group}} + \beta_{2}^{\operatorname{Base}} + \beta_{3}^{\operatorname{Age}} + \beta_{4}^{\operatorname{Z}} + \beta_{5}^{\operatorname{R1}} + \beta_{6}^{\operatorname{R2}} + \epsilon

Y^{\operatorname{Post}} = \beta_{0} + \beta_{1}^{\operatorname{Group}} + \beta_{2}^{\operatorname{Base}} + \beta_{3}^{\operatorname{Age}} + \beta_{4}^{\operatorname{Z}} + \beta_{5}^{\operatorname{R1}} + \beta_{6}^{\operatorname{R2}} + \epsilon

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

Stan

data {
  int<lower=0> J;         // number of schools
  real y[J];              // estimated treatment effects
  real<lower=0> sigma[J]; // standard error of effect estimates
}
parameters {
  real mu;                // population treatment effect
  real<lower=0> tau;      // standard deviation in treatment effects
  vector[J] eta;          // unscaled deviation from mu by school
}
transformed parameters {
  vector[J] theta = mu + tau * eta;        // school treatment effects
}
model {
  target += normal_lpdf(eta | 0, 1);       // prior log-density
  target += normal_lpdf(y | theta, sigma); // log-likelihood
}

df1 <- read.csv("../../../ts.csv")
#> Error in `file()`:
#> ! cannot open the connection
df1$y <- ts(df1$Sales)
#> Error:
#> ! object 'df1' not found
df1$ds <- as.Date(df1$Time.Increment)
#> Error:
#> ! object 'df1' not found

splits <- initial_time_split(df1, prop = 0.5)
#> Error:
#> ! object 'df1' not found
train <- training(splits)
#> Error:
#> ! object 'splits' not found
test <- testing(splits)
#> Error:
#> ! object 'splits' not found
interactive <- TRUE
# Forecasting with auto.arima
library("forecast")
md <- auto.arima(train$y)
#> Error:
#> ! object 'train' not found
fc <- forecast(md, h = 12)
#> Error:
#> ! object 'md' not found

model_fit_arima_no_boost <- arima_reg() %>%
  set_engine(engine = "auto_arima") %>%
  fit(y ~ ds, data = training(splits))
#> Error:
#> ! object 'splits' not found

# Model 2: arima_boost ----
model_fit_arima_boosted <- arima_boost(
  min_n = 2,
  learn_rate = 0.015
) %>%
  set_engine(engine = "auto_arima_xgboost") %>%
  fit(y ~ ds + as.numeric(ds) + factor(month(ds, label = TRUE),
    ordered = F
  ),
  data = training(splits)
  )
#> Error:
#> ! object 'splits' not found

# Model 3: ets ----
model_fit_ets <- exp_smoothing() %>%
  set_engine(engine = "ets") %>%
  fit(y ~ ds, data = training(splits))
#> Error:
#> ! object 'splits' not found

model_fit_lm <- linear_reg() %>%
  set_engine("lm") %>%
  fit(y ~ as.numeric(ds) + factor(month(ds, label = TRUE),
    ordered = FALSE
  ),
  data = training(splits)
  )
#> Error:
#> ! object 'splits' not found

# Model 4: prophet ----
model_fit_prophet <- prophet_reg() %>%
  set_engine(engine = "prophet") %>%
  fit(y ~ ds, data = training(splits))
#> Error:
#> ! object 'splits' not found
# Model 6: earth ----
model_spec_mars <- mars(mode = "regression") %>%
  set_engine("earth")


recipe_spec <- recipe(y ~ ds, data = training(splits)) %>%
  step_date(ds, features = "month", ordinal = FALSE) %>%
  step_mutate(date_num = as.numeric(ds)) %>%
  step_normalize(date_num) %>%
  step_rm(ds)
#> Error:
#> ! object 'splits' not found

wflw_fit_mars <- workflow() %>%
  add_recipe(recipe_spec) %>%
  add_model(model_spec_mars) %>%
  fit(training(splits))
#> Error:
#> ! object 'recipe_spec' not found

models_tbl <- modeltime_table(
  model_fit_arima_no_boost,
  model_fit_arima_boosted,
  model_fit_ets,
  model_fit_prophet,
  model_fit_lm,
  wflw_fit_mars
)
#> Error:
#> ! object 'model_fit_arima_no_boost' not found

calibration_tbl <- models_tbl %>%
  modeltime_calibrate(new_data = testing(splits))
#> Error:
#> ! object 'splits' not found

calibration_tbl %>%
  modeltime_forecast(
    new_data    = testing(splits),
    actual_data = df1
  ) %>%
  plot_modeltime_forecast(
    .legend_max_width = 25, # For mobile screens
    .interactive      = interactive
  )
#> Error:
#> ! object 'splits' not found
calibration_tbl %>%
  modeltime_accuracy() %>%
  table_modeltime_accuracy(
    .interactive = FALSE
  )
#> Error:
#> ! object 'calibration_tbl' not found

refit_tbl <- calibration_tbl %>%
  modeltime_refit(data = df1)
#> Error:
#> ! object 'calibration_tbl' not found

refit_tbl %>%
  modeltime_forecast(h = "4 years", actual_data = df1) %>%
  plot_modeltime_forecast(
    .legend_max_width = 10, # For mobile screens
    .interactive      = TRUE
  )
#> Error:
#> ! object 'refit_tbl' not found

plot(greybox::forecast(smooth::adam(df1$y, h = 12, holdout = TRUE)))
#> Error in `h()`:
#> ! error in evaluating the argument 'x' in selecting a method for function 'plot': object 'df1' not found

plot(greybox::forecast(smooth::es(df1$y, h = 12, holdout = TRUE, silent = FALSE)))
#> Error in `h()`:
#> ! error in evaluating the argument 'x' in selecting a method for function 'plot': object 'df1' not found

s1 <- bayesforecast::stan_naive(
  ts = df1$y, chains = 4,
  iter = 4000, cores = 8
)
#> Error:
#> ! object 'df1' not found

plot(s1) +
  theme_bw()
#> Error in `h()`:
#> ! error in evaluating the argument 'x' in selecting a method for function 'plot': object 's1' not found
check_residuals(s1) +
  theme_light()
#> Error:
#> ! object 's1' not found
autoplot(
  object = forecast(s1, h = 12, biasadj = TRUE, PI = TRUE),
  include = 100
) +
  theme_bw()
#> Error:
#> ! object 's1' not found
autoplot(
  object = forecast(s1, h = 52, biasadj = TRUE, PI = TRUE),
  include = 100
) +
  theme_bw()
#> Error:
#> ! object 's1' not found
ztable(summary(s1))
#> Error in `h()`:
#> ! error in evaluating the argument 'object' in selecting a method for function 'summary': object 's1' not found
(meanf(df1$y))
#> Error:
#> ! object 'df1' not found
(forecast(s1, h = 12))
#> Error:
#> ! object 's1' not found

autoplot(
  object = forecast(s1, h = 6, biasadj = TRUE, PI = TRUE),
  include = 100
) +
  theme_bw()
#> Error:
#> ! object 's1' not found

df <- as.data.frame(df1)
#> Error:
#> ! object 'df1' not found
m <- prophet(df1, growth = "linear", yearly.seasonality = "auto")
#> Error:
#> ! object 'df1' not found
future <- make_future_dataframe(m,
  periods = 60,
  freq = "months",
  include_history = TRUE
)
#> Error:
#> ! object 'm' not found
forecast <- predict(m, future)
#> Error:
#> ! object 'm' not found
plot(m, forecast) +
  theme_bw()
#> Error in `h()`:
#> ! error in evaluating the argument 'x' in selecting a method for function 'plot': object 'm' not found

plot(greybox::forecast(smooth::auto.adam(df1$y,
  h = 12, holdout = TRUE, ic = "AICc",
  regressors = "select"
)))
#> Error in `h()`:
#> ! error in evaluating the argument 'x' in selecting a method for function 'plot': object 'df1' not found
adamAutoARIMAAir <- auto.adam(df1$y, h = 50)
#> Error:
#> ! object 'df1' not found

plot(greybox::forecast(adam(df1$y,
  main = "Parametric prediction interval",
  h = 5, sim = 1000, lev = .99, interval = "prediction"
)))
#> Error in `h()`:
#> ! error in evaluating the argument 'x' in selecting a method for function 'plot': object 'df1' not found
plot(greybox::forecast(auto.adam(df1$y,
  h = 5, interval = "complete", nsim = 100,
  main = "Complete prediction interval"
)))
#> Error in `h()`:
#> ! error in evaluating the argument 'x' in selecting a method for function 'plot': object 'df1' not found

plot(greybox::forecast(adam(df1$y,
  c("CCN", "ANN", "AAN", "AAdN"),
  h = 10, holdout = TRUE,
  ic = "AICc"
)))
#> Error in `h()`:
#> ! error in evaluating the argument 'x' in selecting a method for function 'plot': object 'df1' not found

plot(greybox::forecast(adam(df1$y,
  model = "NNN", lags = 7,
  orders = c(0, 1, 1), constant = TRUE,
  h = 5, interval = "complete", nsim = 100
)))
#> Error in `h()`:
#> ! error in evaluating the argument 'x' in selecting a method for function 'plot': object 'df1' not found

plot(greybox::forecast((adam(df1$y,
  model = "NNN", lags = c(24, 24 * 7, 24 * 365),
  orders = list(ar = c(3, 2, 2, 2), i = c(2, 1, 1, 1), ma = c(3, 2, 2, 2), select = TRUE),
  initial = "backcasting"
))))
#> Error in `h()`:
#> ! error in evaluating the argument 'x' in selecting a method for function 'plot': object 'df1' not found

adamARIMA
#> Error:
#> ! object 'adamARIMA' not found
# Apply models
adamPoolBJ <- vector("list", 3)
adamPoolBJ[[1]] <- adam(df1$y, "ZZN",
  h = 10, holdout = TRUE,
  ic = "BICc"
)
#> Error:
#> ! object 'df1' not found
adamPoolBJ[[2]] <- adam(df1$y, "NNN",
  orders = list(ar = 3, i = 2, ma = 3, select = TRUE),
  h = 10, holdout = TRUE,
  ic = "BICc"
)
#> Error:
#> ! object 'df1' not found
adamPoolBJ[[3]] <- adam(df1$y, "MMN",
  h = 10, holdout = TRUE,
  ic = "BICc",
  regressors = "select"
)
#> Error:
#> ! object 'df1' not found

# Extract BICc values
adamsICs <- sapply(adamPoolBJ, BICc)
#> Error in `UseMethod()`:
#> ! no applicable method for 'logLik' applied to an object of class "NULL"

# Calculate weights
adamsICWeights <- adamsICs - min(adamsICs)
#> Error:
#> ! object 'adamsICs' not found
adamsICWeights[] <- exp(-0.5 * adamsICWeights) /
  sum(exp(-0.5 * adamsICWeights))
#> Error:
#> ! object 'adamsICWeights' not found
names(adamsICWeights) <- c("ETS", "ARIMA", "ETSX")
#> Error:
#> ! object 'adamsICWeights' not found
round(adamsICWeights, 3)
#> Error:
#> ! object 'adamsICWeights' not found

adamPoolBJForecasts <- vector("list", 3)
# Produce forecasts from the three models
for (i in 1:3) {
  adamPoolBJForecasts[[i]] <- forecast(adamPoolBJ[[i]],
    h = 10, interval = "pred"
  )
}
#> Error in `forecast.NULL()`:
#> ! argument "new_data" is missing, with no default
# Produce combined conditional means and prediction intervals
finalForecast <- cbind(
  sapply(
    adamPoolBJForecasts,
    "[[", "mean"
  ) %*% adamsICWeights,
  sapply(
    adamPoolBJForecasts,
    "[[", "lower"
  ) %*% adamsICWeights,
  sapply(
    adamPoolBJForecasts,
    "[[", "upper"
  ) %*% adamsICWeights
)
#> Error:
#> ! object 'adamsICWeights' not found
# Give the appropriate names
colnames(finalForecast) <- c(
  "Mean", "Lower bound (2.5%)",
  "Upper bound (97.5%)"
)
#> Error:
#> ! object 'finalForecast' not found
# Transform the table in the ts format (for convenience)
finalForecast <- ts(finalForecast,
  start = start(adamPoolBJForecasts[[i]]$mean)
)
#> Error:
#> ! object 'finalForecast' not found
finalForecast
#> Error:
#> ! object 'finalForecast' not found

graphmaker(df1$y, finalForecast[, 1],
  lower = finalForecast[, 2], upper = finalForecast[, 3],
  level = 0.95
)
#> Error:
#> ! object 'finalForecast' not found

plot(forecast(forecastHybrid::hybridModel(df1$y,
  models = "aen",
  weights = "equal",
  cvHorizon = 8, num.cores = 4
),
h = 5
))
#> Error in `h()`:
#> ! error in evaluating the argument 'x' in selecting a method for function 'plot': object 'df1' not found

plot(forecast(forecastHybrid::hybridModel(df1$y,
  models = "fnst",
  weights = "equal",
  errorMethod = "RMSE"
)))
#> Error in `h()`:
#> ! error in evaluating the argument 'x' in selecting a method for function 'plot': object 'df1' not found

oesModel <- oes(df1$y, model = "YYY", occurrence = "auto")
#> Error:
#> ! object 'df1' not found

y0 <- stan_sarima(df1$y,
  refresh = 0, verbose = FALSE,
  open_progress = FALSE
)
#> Error:
#> ! object 'df1' not found
autoplot(forecast(y0, 5, .99), "red") +
  ggplot2::theme_bw()
#> Error:
#> ! object 'y0' not found

df1$month <- lubridate::month(df1$ds)
#> Error:
#> ! object 'df1' not found
gam1 <- mgcv::gam(Sales ~ s(month, bs = "cr", k = 12),
  data = df1, family = gaussian,
  correlation = SARIMA(form = ~month, p = 1),
  method = "REML"
) |>
  mgcv::plot.gam(lwd = 3, lty = 1, col = "#d46c5b")
#> Error:
#> ! object 'df1' not found

ts_plot(df1,
  title = "US Monthly Natural Gas Consumption",
  Ytitle = "Billion Cubic Feet"
)
#> Error:
#> ! object 'df1' not found

months <- c("2022-01-01", "2022-03-01", "2022-06-01")
ubereats <- c(7327.55, 4653.53, 4833.21)
doordash <- c(1304.54, 2000.35, 1643.58)
grubhub <- c(1199.85, 941.68, 623.27)
total <- c(7222.85, 7464.68, 7100.06)

df <- data.frame(months, ubereats, doordash, grubhub, total)
df$months <- as.Date(df$months)


ts_plot(df,
  title = "US Monthly Natural Gas Consumption",
  Ytitle = "Billion Cubic Feet"
)

df$total <- ts(df$total)
plot(greybox::forecast(adam(df$total,
  h = 5, holdout = TRUE, ic = "AICc",
  regressors = "select"
)))
#> Error in `h()`:
#> ! error in evaluating the argument 'x' in selecting a method for function 'plot': The number of in-sample observations is not positive. Cannot do anything.

m <- prophet(df, growth = "linear", yearly.seasonality = "auto")
#> Error in `fit.prophet()`:
#> ! Dataframe must have columns 'ds' and 'y' with the dates and values respectively.
future <- make_future_dataframe(m,
  periods = 60,
  freq = "months",
  include_history = TRUE
)
#> Error:
#> ! object 'm' not found
forecast <- predict(m, future)
#> Error:
#> ! object 'm' not found
plot(m, forecast) +
  theme_bw()
#> Error in `h()`:
#> ! error in evaluating the argument 'x' in selecting a method for function 'plot': object 'm' not found

# Plotting actual vs. fitted and forecasted
test_forecast(actual = df1$y, forecast.obj = fc, test = test$y)
#> Error:
#> ! object 'fc' not found
plot(fc)
#> Error in `h()`:
#> ! error in evaluating the argument 'x' in selecting a method for function 'plot': object 'fc' not found

import delimited "../../../ts.csv", clear
encode timeincrement, generate(t)
parallel initialize 8, f
tsset t
tsline sales
ac sales
dfuller sales, regress trend
ac d.sales
ac sales
arima sales, arima(1,1,1)
estat ic
arima sales, arima(1,0,0)
estat ic
tsappend, add(30)
browse
predict fsales, y dynamic(m(2022m6))
label variable fsales "forecasted sales"
tsline sales fsales
graph export "forecastedsales.png", replace
#> . import delimited "../../../ts.csv", clear
#> file ../../../ts.csv not found
#> r(601);
#> 
#> r(601);

Python

from pandas import DataFrame
import statsmodels.api as sm
import matplotlib as plot

Stock_Market = {'Year': [2017,2017,2017,2017,2017,
                         2017,2017,2017,2017,2017,
                         2017,2017,2016,2016,2016,
                         2016,2016,2016,2016,2016,
                         2016,2016,2016,2016],
                'Month': [12, 11,10,9,8,7,6,5,4,
                          3,2,1,12,11,10,9,8,7,6,
                          5,4,3,2,1],
                'Interest_Rate':[2.75,2.5,2.5,2.5,2.5,2.5,
                                 2.5,2.25,2.25,2.25,2,2,2,1.75,1.75,
                                 1.75,1.75,1.75,1.75,1.75,1.75,1.75,1.75,1.75],
                'Unemployment_Rate':[5.3,5.3,5.3,5.3,5.4,5.6,
                                     5.5,5.5,5.5,5.6,5.7,5.9,6,5.9,5.8,6.1,
                                     6.2,6.1,6.1,6.1,5.9,6.2,6.2,6.1],
                'Stock_Index_Price': [1464,1394,1357,1293,1256,1254,1234,1195,1159,1167,1130,
                                      1075,1047,965,943,958,971,949,884,866,876,822,704,719]
                }

df = DataFrame(Stock_Market,columns=['Year','Month','Interest_Rate',
                                     'Unemployment_Rate','Stock_Index_Price'])

X = df[['Interest_Rate','Unemployment_Rate']]

# here we have 2 variables for the multiple linear regression. If you just want to use one variable for simple linear regression, then use X = df['Interest_Rate'] for example

Y = df['Stock_Index_Price']

X = sm.add_constant(X) # adding a constant

model = sm.OLS(Y, X).fit()
predictions = model.predict(X)

print_model = model.summary()
print(print_model)
#>                             OLS Regression Results                            
#> ==============================================================================
#> Dep. Variable:      Stock_Index_Price   R-squared:                       0.898
#> Model:                            OLS   Adj. R-squared:                  0.888
#> Method:                 Least Squares   F-statistic:                     92.07
#> Date:                Fri, 14 Aug 2026   Prob (F-statistic):           4.04e-11
#> Time:                        12:21:34   Log-Likelihood:                -134.61
#> No. Observations:                  24   AIC:                             275.2
#> Df Residuals:                      21   BIC:                             278.8
#> Df Model:                           2                                         
#> Covariance Type:            nonrobust                                         
#> =====================================================================================
#>                         coef    std err          t      P>|t|      [0.025      0.975]
#> -------------------------------------------------------------------------------------
#> const              1798.4040    899.248      2.000      0.059     -71.685    3668.493
#> Interest_Rate       345.5401    111.367      3.103      0.005     113.940     577.140
#> Unemployment_Rate  -250.1466    117.950     -2.121      0.046    -495.437      -4.856
#> ==============================================================================
#> Omnibus:                        2.691   Durbin-Watson:                   0.530
#> Prob(Omnibus):                  0.260   Jarque-Bera (JB):                1.551
#> Skew:                          -0.612   Prob(JB):                        0.461
#> Kurtosis:                       3.226   Cond. No.                         394.
#> ==============================================================================
#> 
#> Notes:
#> [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

R

fit <- glm(mpg ~ cyl + disp, mtcars,
  family = gaussian()
)
# show the theoretical model
equatiomatic::extract_eq(fit)
#> Error in `loadNamespace()`:
#> ! there is no package called 'equatiomatic'

Stata

sysuse auto2, clear
parallel initialize 8, f
mfp: glm price mpg
twoway (fpfitci price mpg, estcmd(glm) fcolor(dkorange%20) alcolor(%40))  || scatter price mpg, mcolor(dkorange) scale(0.75)
graph export "mfp.png", replace
#> . sysuse auto2, clear
#> (1978 automobile data)
#> 
#> . parallel initialize 8, f
#> -initialize- invalid parallel subcommand (help parallel)
#> r(198);
#> 
#> r(198);

clear
set obs 100
/**
Generate variables from a normal distribution like before
*/
generate x = rnormal(0, 1)
generate y = rnormal(0, 1)
/**
We set up our model here

*/
parallel initialize 8, f
bayesmh y x, likelihood(normal({var})) prior({var}, normal(0, 10)) ///
prior({y:}, normal(0, 10)) rseed(1031) saving(coutput_pred, replace) mcmcsize(1000)
/**
We use the bayespredict command to make predictions from the model
*/
bayespredict (mean:@mean({_resid})) (var:@variance({_resid})), ///
rseed(1031) saving(coutput_pred, replace)
/**
Then we calculate the posterior predictive P-values
*/
bayesstats ppvalues {mean} using coutput_pred
#> . clear
#> 
#> . set obs 100
#> Number of observations (_N) was 0, now 100.
#> 
#> . /**
#> > Generate variables from a normal distribution like before
#> > */
#> . generate x = rnormal(0, 1)
#> 
#> . generate y = rnormal(0, 1)
#> 
#> . /**
#> > We set up our model here
#> > 
#> > */
#> . parallel initialize 8, f
#> -initialize- invalid parallel subcommand (help parallel)
#> r(198);
#> 
#> r(198);

library(lme4)
library(simr)

# Toy model
fm = lmer(y ~ x + (x | g), data = simdata)

# Extend sample size of `g`
fm_extended_g = extend(fm, along = 'g', n = 4)
# 4 levels of g
pwcurve_4g = powerCurve(fm_extended_g, fixed('x'), along = 'g', breaks = 4,
                        nsim = 50, seed = 123,
                        # No progress bar
                        progress = FALSE)
# 6 levels of g
# Create a destination object using any of the power curves above.
all_pwcurve = pwcurve_4g
# Combine results
all_pwcurve$ps = c(pwcurve_4g$ps[1])

# Combine the different numbers of levels.
all_pwcurve$xval = c(pwcurve_4g$nlevels)

print(all_pwcurve)
#> Power for predictor 'x', (95% confidence interval),
#> by number of levels in g:
#>       4: 44.00% (29.99, 58.75) - 40 rows
#> 
#> (0 errors, 1 warning, 40 messages)
#> 
#> Time elapsed: 0 h 0 m 7 s

plot(all_pwcurve, xlab = 'Levels of g')