# 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