ST308 Bayesian Inference

Class 8


Philipp Sterzinger

p.sterzinger@lse.ac.uk

21 March 2025

Recap

Normal Linear Model

Assume that \[ y_i = \beta_0 + x_i^\top \beta + \varepsilon_i, \]

for intercept \(\beta_0\), and regression coefficients \(\beta \in \Re^p\) and where, conditional on \(X\), the errors are indepent

\[ \varepsilon_i \mid_{x_i} \sim \mathcal{N}(0, \sigma^2) \]

In frequentist analysis, we condition on the covariates, and treat \(\beta_0, \beta, \sigma^2\) as fixed

Bayesian Normal Linear Model

We now assume that

\[ y_i = \beta_0 + x_i^\top \beta + \varepsilon_i, \]

and

\[ \begin{aligned} \varepsilon_i \mid_{x_i} &\overset{\textrm{i.i.d.}}{\sim} \mathcal{N}(0, \sigma^2) \\ \beta_0 &\sim \pi_{\beta_0}(\beta_0) \\ \beta &\sim \pi_{\beta}(\beta) \\ \sigma^2 &\sim \pi_{\sigma^2}(\sigma^2) \end{aligned} \]

so that

\[ y \mid_{X, \beta_0, \beta, \sigma^2} \sim \mathcal{N}(\beta_0 + X \beta, \sigma^2 I_n) \]

Bayesian Normal Linear Model

Oftentimes, we are unable to get a closed form expression for the posterior, so that we have to resort to MCMC for inference.

We will use rstanarm for this

Stan

Stan is a probabilistic programming language for statistical inference written in C++. Wikipedia

rstanarm

rstanarm is an R package that emulates other R model-fitting functions but uses Stan (via the rstan package) for the back-end estimation. rstanrm documentation

Tutorial

rstanarm

install.packages("rstanarm")

To get started, and for future reference, check out the vignettes of rstanarm

Today we consider the bayesian linear model, but like with glm, you can estimate all kinds of GLMs such as bayesian logistic regression

stan_glm

The workhorse for estimation purposes is stan_glm

Similar in functionality and syntax to glm but allows bayesian approach

It is good practice to check the documentation of the functions we are using

library(rstanarm) 

?stan_glm

stan_glm

Key inputs

  • formula: An R formula as you would use in lm, glm, glmnet... such as y ~ 1
  • data: A R data.frame that contains our data, i.e. a response y and some covariates X
  • family: The family of GLM that we want to use, gaussian for the normal linear model, binomial for logistic regression. This also determines the default link function
    • For each family, we need to specify the link function, for the linear model, we need the identity link, for logistic regression we need the logit
    • For the linear normal model, we would thus use family = gaussian(link = "identity")
  • prior: The prior on the regression coefficients \(\beta\). The default priors depend on the model family, you can check these with default_prior_coef(family).
  • prior_intercept: The prior on the intercept parameter \(\beta_0\). Check with default_prior_intercept(family)
  • prior_aux: The prior on the auxiliary parameter of the GLM, e.g. \(\sigma^2\) for the gaussian(link = "identity"). Check the vignette for the model that you want to use

stan_glm

Output

A stanreg object

Useful functions:

  • summary
  • prior_summary
  • posterior_interval
  • posterior_predict
  • launch_shinystan

To look at the posterior sample that you got from a stan_glm call, do as.matrix(your_stanreg_object)

Example

\[ y \sim \mathcal{N}(\mu, \sigma^2) \]

\[ \mu \sim \mathcal{N}(0, 10), \quad \sigma^2 \sim \mathcal{N}^{+}(0, 10) \]

Code 
## Simulate data 
set.seed(1)
N=100
mu = 5
sigma = 4
obs = rnorm(N,mu,sigma)
toy_dat = data.frame(obs)

toy <- stan_glm(obs ~ 1, data = toy_dat,
                prior_intercept = normal(0, 10), prior_aux = normal(0,10))  

SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 3.3e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.33 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.01 seconds (Warm-up)
Chain 1:                0.016 seconds (Sampling)
Chain 1:                0.026 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 4e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 0.01 seconds (Warm-up)
Chain 2:                0.014 seconds (Sampling)
Chain 2:                0.024 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 4e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.009 seconds (Warm-up)
Chain 3:                0.015 seconds (Sampling)
Chain 3:                0.024 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 3e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 0.01 seconds (Warm-up)
Chain 4:                0.016 seconds (Sampling)
Chain 4:                0.026 seconds (Total)
Chain 4:

Example

Code 
summary(toy) 

Model Info:
 function:     stan_glm
 family:       gaussian [identity]
 formula:      obs ~ 1
 algorithm:    sampling
 sample:       4000 (posterior sample size)
 priors:       see help('prior_summary')
 observations: 100
 predictors:   1

Estimates:
              mean   sd   10%   50%   90%
(Intercept) 5.4    0.4  5.0   5.4   5.9  
sigma       3.6    0.3  3.3   3.6   4.0  

Fit Diagnostics:
           mean   sd   10%   50%   90%
mean_PPD 5.4    0.5  4.8   5.4   6.1  

The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).

MCMC diagnostics
              mcse Rhat n_eff
(Intercept)   0.0  1.0  2839 
sigma         0.0  1.0  2767 
mean_PPD      0.0  1.0  3277 
log-posterior 0.0  1.0  1685 

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).
Code 
# By hand summary 
coefs <- as.matrix(toy) 
means <- colMeans(coefs) 
sds <- apply(coefs, 2, sd) 
quants <- apply(coefs, 2, function(x){quantile(x, probs = c(.1, .5, .9))}) 

rbind(means, sds, quants)
      (Intercept)     sigma
means   5.4315202 3.6465385
sds     0.3514498 0.2639027
10%     4.9770134 3.3257507
50%     5.4347150 3.6308075
90%     5.8755196 4.0009727

Example

Code 
prior_summary(toy) 
Priors for model 'toy' 
------
Intercept (after predictors centered)
 ~ normal(location = 0, scale = 10)

Auxiliary (sigma)
 ~ half-normal(location = 0, scale = 10)
------
See help('prior_summary.stanreg') for more details
Code 
ci95 <- posterior_interval(toy, prob = 0.95)
round(ci95,2)
            2.5% 97.5%
(Intercept) 4.73  6.12
sigma       3.17  4.19

Example

Code 
library(MASS) 
library(pracma)

X <- coefs 

density_est <- kde2d(X[,1], X[,2], n = 100)
f <- function(x0, y0) {
  interp2(density_est$x, density_est$y, density_est$z, x0, y0, method = "linear")
}


vx=seq(4, 7,length=201)
vy=seq(3, 5,length=201)
z=outer(vx,vy, f)

xhist <- hist(X[,1], plot=FALSE)
yhist <- hist(X[,2], plot=FALSE)

top <- max(c(xhist$density, yhist$density))
nf <- layout(matrix(c(2,0,1,3),2,2,byrow=TRUE), c(3,1), c(1,3), TRUE)
par(mar=c(3,3,1,1))
image(vx,vy,z,col=rev(heat.colors(101)))
contour(vx,vy,z,col="blue",add=TRUE)
points(X,cex=.2)
par(mar=c(0,3,1,1))
barplot(xhist$density, axes=FALSE, ylim=c(0, top), space=0,col="light green")
lines((density(X[,1])$x-xhist$breaks[1])/diff(xhist$breaks)[1],
density(X[,1])$y,col="red")
par(mar=c(3,0,1,1))
barplot(yhist$density, axes=FALSE, xlim=c(0, top), space=0, 
horiz=TRUE,col="light green")
lines(density(X[,2])$y,(density(X[,2])$x-yhist$breaks[1])/
diff(yhist$breaks)[1],col="red")

Activity 1

Repeat the above exercise for the following model \[ y_i\sim \text{Poisson}(\lambda), \;\;i=1,...,N, \] with priors \[ \lambda \sim \text{Normal}_{>0}(0,100), \]

Data can be simulated in the following way:

Code 
set.seed(1)
N = 100
lambda = 2
y = rpois(N, exp(lambda))
mean(y)
[1] 7.45

Activity 1

Code 
toy_dat = data.frame(y)

toy <- stan_glm(y~1,data=toy_dat,family=poisson,
                prior_intercept = normal(0, 100))

SAMPLING FOR MODEL 'count' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 1.5e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.012 seconds (Warm-up)
Chain 1:                0.021 seconds (Sampling)
Chain 1:                0.033 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'count' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 1.1e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.11 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 0.013 seconds (Warm-up)
Chain 2:                0.02 seconds (Sampling)
Chain 2:                0.033 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'count' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 4e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.014 seconds (Warm-up)
Chain 3:                0.023 seconds (Sampling)
Chain 3:                0.037 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'count' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 5e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 0.012 seconds (Warm-up)
Chain 4:                0.02 seconds (Sampling)
Chain 4:                0.032 seconds (Total)
Chain 4:
Code 
prior_summary(toy)
Priors for model 'toy' 
------
Intercept (after predictors centered)
 ~ normal(location = 0, scale = 100)
------
See help('prior_summary.stanreg') for more details
Code 
toy
stan_glm
 family:       poisson [log]
 formula:      y ~ 1
 observations: 100
 predictors:   1
------
            Median MAD_SD
(Intercept) 2.0    0.0   

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

Activity 1

Code 
ci95 <- posterior_interval(toy, prob = 0.95)
round(ci95,2)
            2.5% 97.5%
(Intercept) 1.93  2.08
Code 
#launch_shinystan(toy, ppd = FALSE)

Activity 2

Work with Boston dataset from the MASS library and fit a Bayesian linear regression model for the response variable ‘medv’ (median value home price) based on the independent variables.

The data can be loaded with the code below. A linear regression model (non-Bayesian) is also fit below for reference.

Code 
library(MASS)
#summary(Boston)
Code 
lreg <- lm(medv ~ ., data = Boston)
summary(lreg)

Call:
lm(formula = medv ~ ., data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.595  -2.730  -0.518   1.777  26.199 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.646e+01  5.103e+00   7.144 3.28e-12 ***
crim        -1.080e-01  3.286e-02  -3.287 0.001087 ** 
zn           4.642e-02  1.373e-02   3.382 0.000778 ***
indus        2.056e-02  6.150e-02   0.334 0.738288    
chas         2.687e+00  8.616e-01   3.118 0.001925 ** 
nox         -1.777e+01  3.820e+00  -4.651 4.25e-06 ***
rm           3.810e+00  4.179e-01   9.116  < 2e-16 ***
age          6.922e-04  1.321e-02   0.052 0.958229    
dis         -1.476e+00  1.995e-01  -7.398 6.01e-13 ***
rad          3.060e-01  6.635e-02   4.613 5.07e-06 ***
tax         -1.233e-02  3.760e-03  -3.280 0.001112 ** 
ptratio     -9.527e-01  1.308e-01  -7.283 1.31e-12 ***
black        9.312e-03  2.686e-03   3.467 0.000573 ***
lstat       -5.248e-01  5.072e-02 -10.347  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.745 on 492 degrees of freedom
Multiple R-squared:  0.7406,    Adjusted R-squared:  0.7338 
F-statistic: 108.1 on 13 and 492 DF,  p-value: < 2.2e-16

Activity 2

The model and priors are given below

\[ y\sim N(\beta_0 + X\beta,\sigma^2) \] \[ \sigma \sim N(0,10^2),\;\;\;\sigma>0 \]

\[ \beta_i \sim N(0,10^2), \;\;\;i=0,\dots,p + 1. \]

Activity 2

Code 
blreg <- stan_glm(medv ~ ., data = Boston, prior = normal(0, 10),
                prior_intercept = normal(0, 10), prior_aux = normal(0,10))

SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 1.4e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.644 seconds (Warm-up)
Chain 1:                0.3 seconds (Sampling)
Chain 1:                0.944 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 9e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 0.564 seconds (Warm-up)
Chain 2:                0.304 seconds (Sampling)
Chain 2:                0.868 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 6e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.729 seconds (Warm-up)
Chain 3:                0.27 seconds (Sampling)
Chain 3:                0.999 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 6e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 0.605 seconds (Warm-up)
Chain 4:                0.322 seconds (Sampling)
Chain 4:                0.927 seconds (Total)
Chain 4:
Code 
prior_summary(blreg)
Priors for model 'blreg' 
------
Intercept (after predictors centered)
 ~ normal(location = 0, scale = 10)

Coefficients
 ~ normal(location = [0,0,0,...], scale = [10,10,10,...])

Auxiliary (sigma)
 ~ half-normal(location = 0, scale = 10)
------
See help('prior_summary.stanreg') for more details

Activity 2

Code 
blreg
stan_glm
 family:       gaussian [identity]
 formula:      medv ~ .
 observations: 506
 predictors:   14
------
            Median MAD_SD
(Intercept)  34.9    5.0 
crim         -0.1    0.0 
zn            0.0    0.0 
indus         0.0    0.1 
chas          2.7    0.9 
nox         -15.5    3.6 
rm            3.8    0.4 
age           0.0    0.0 
dis          -1.4    0.2 
rad           0.3    0.1 
tax           0.0    0.0 
ptratio      -0.9    0.1 
black         0.0    0.0 
lstat        -0.5    0.1 

Auxiliary parameter(s):
      Median MAD_SD
sigma 4.8    0.2   

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg
Code 
ci95 <- posterior_interval(blreg, prob = 0.95)
round(ci95, 2)
              2.5% 97.5%
(Intercept)  25.30 44.54
crim         -0.17 -0.04
zn            0.02  0.07
indus        -0.11  0.13
chas          0.96  4.35
nox         -22.21 -8.76
rm            2.98  4.63
age          -0.03  0.02
dis          -1.82 -1.05
rad           0.17  0.43
tax          -0.02 -0.01
ptratio      -1.17 -0.69
black         0.00  0.01
lstat        -0.63 -0.42
sigma         4.47  5.07

Philipp Sterzinger - ST308 Week 9