# Bayesian data analysis - PyStan demos

Author: Aki Vehtari, Tuomas Sivula

Examples how to use Stan in Python with PyStan. This notebook assumes basic knowledge of Bayesian inference and MCMC. The Stan models are stored in separate .stan-files.

## Setting up the PyStan environment

We begin by importing the PyStan module as well at the matplotlib module for basic graphics facilities.

We also import some utilities by Michael Betancourt and introduced in PyStan workflow case study in module `stan_utility`.

``````

In [1]:

%matplotlib inline

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

import pystan

``````
``````

In [2]:

# add utilities directory to path
import os, sys
util_path = os.path.abspath(os.path.join(os.path.pardir, 'utilities_and_data'))
if util_path not in sys.path and os.path.exists(util_path):
sys.path.insert(0, util_path)

# import from utilities
import stan_utility

``````
``````

In [3]:

# edit default plot settings
plt.rc('font', size=12)

``````

## Bernoulli model

Toy data with sequence of failures (0) and successes (1). We would like to learn about the unknown probability of success.

``````

In [4]:

data = dict(N=10, y=[0,1,0,0,1,1,1,0,1,0])

``````

Bernoulli model with a Beta(1,1) (uniform) prior

``````

In [5]:

with open('bern.stan') as file:

``````
``````

// Bernoulli model
data {
int<lower=0> N;
int<lower=0,upper=1> y[N];
}
parameters {
real<lower=0,upper=1> theta;
}
model {
theta ~ beta(1,1);
y ~ bernoulli(theta);
}

``````

Given the Stan program we then use the `compile_model` method of `stan_utility` module to compile the Stan program into a C++ executable. This utility function automatically saves a cached version of the compiled model to the disk for possible future use.

``````

In [6]:

model = stan_utility.compile_model('bern.stan')

``````
``````

INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_4586f2dc76604848221fafe6413762a9 NOW.

``````

Getting the model again with the utility function uses the cached version automatically.

``````

In [7]:

del model
model = stan_utility.compile_model('bern.stan')

``````
``````

Using cached StanModel

``````

Sample form the posterior, show the summary and plot the histogram of the posterior draws. We recommend explicitly specifying the seed of Stan's random number generator, as we have done here, so that we can reproduce these exactly results in the future, at least when using the same machine, operating system, and interface. This is especially helpful for the more subtle pathologies that may not always be found, which results in seemingly stochastic behavior.

``````

In [8]:

fit = model.sampling(data=data, seed=194838)
print(fit)
samples = fit.extract(permuted=True)
plt.hist(samples['theta'], 50);

``````
``````

Inference for Stan model: anon_model_4586f2dc76604848221fafe6413762a9.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.

mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
theta    0.5  3.7e-3   0.14   0.24    0.4    0.5    0.6   0.76   1435    1.0
lp__   -8.85    0.02   0.72 -10.91  -9.05  -8.57  -8.37  -8.32   1704    1.0

Samples were drawn using NUTS at Tue Oct 23 12:58:08 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).

``````

## Binomial model

Instead of sequence of 0's and 1's, we can summarize the data with the number of experiments and the number successes:

``````

In [9]:

data = dict(N=10, y=8)

``````

And then we use Binomial model with Beta(1,1) prior for the probability of success.

``````

In [10]:

with open('binom.stan') as file:

``````
``````

// Binomial model with beta(1,1,) prior
data {
int<lower=0> N;
int<lower=0> y;
}
parameters {
real<lower=0,upper=1> theta;
}
model {
theta ~ beta(1,1);
y ~ binomial(N,theta);
}

``````

Sample from the posterior and plot the posterior. The histogram should look similar as in the Bernoulli case.

``````

In [11]:

model = stan_utility.compile_model('binom.stan')
fit = model.sampling(data=data, seed=194838)
samples = fit.extract(permuted=True)
plt.hist(samples['theta'], 50);

``````
``````

INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_a05259dcade68de9675d2ad9b1145f54 NOW.

``````

Re-run the model with a new data. Here the same compiled model object is reused.

``````

In [12]:

data = dict(N=10, y=10)
fit = model.sampling(data=data, seed=194838)
samples = fit.extract(permuted=True)
plt.hist(samples['theta'], 50);

``````
``````

``````

## Comparison of two groups with Binomial

An experiment was performed to estimate the effect of beta-blockers on mortality of cardiac patients. A group of patients were randomly assigned to treatment and control groups:

• out of 674 patients receiving the control, 39 died
• out of 680 receiving the treatment, 22 died
``````

In [13]:

data = dict(N1=674, y1=39, N2=680, y2=22)

``````

To analyse whether the treatment is useful, we can use Binomial model for both groups and compute odds-ratio:

``````

In [14]:

with open('binom2.stan') as file:

``````
``````

//  Comparison of two groups with Binomial
data {
int<lower=0> N1;
int<lower=0> y1;
int<lower=0> N2;
int<lower=0> y2;
}
parameters {
real<lower=0,upper=1> theta1;
real<lower=0,upper=1> theta2;
}
model {
theta1 ~ beta(1,1);
theta2 ~ beta(1,1);
y1 ~ binomial(N1,theta1);
y2 ~ binomial(N2,theta2);
}
generated quantities {
real oddsratio;
oddsratio = (theta2/(1-theta2))/(theta1/(1-theta1));
}

``````

Sample from the posterior and plot the posterior

``````

In [15]:

model = stan_utility.compile_model('binom2.stan')
fit = model.sampling(data=data, seed=194838)
samples = fit.extract(permuted=True)
plt.hist(samples['oddsratio'], 50);

``````
``````

INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_f7225537f3d9d3af861d4f45c7f8da01 NOW.

``````

## Linear Gaussian model

The following file has Kilpisjärvi summer month temperatures 1952-2013:

``````

In [16]:

data_path = os.path.abspath(
os.path.join(
os.path.pardir,
'utilities_and_data',
'kilpisjarvi-summer-temp.csv'
)
)
d = np.loadtxt(data_path, dtype=np.double, delimiter=';', skiprows=1)
x = d[:, 0]
y = d[:, 4]
N = len(x)
xpred = 2016

``````

Plot the data

``````

In [17]:

# make slightly wider figure
wide_figsize = plt.rcParams['figure.figsize'].copy()
wide_figsize[0] *= 1.4
plt.figure(figsize=wide_figsize)
plt.scatter(x, y)
plt.xlabel('year')
plt.ylabel('summer temperature at Kilpisjarvi');

``````
``````

``````

To analyse whether the average summer month temperature is rising, we use a linear model with Gaussian model for the unexplained variation.

### Gaussian linear model with adjustable priors

The folloing Stan code allows also setting hyperparameter values as data allowing easier way to use different priors in different analyses:

``````

In [18]:

with open('lin.stan') as file:

``````
``````

// Gaussian linear model with adjustable priors
data {
int<lower=0> N; // number of data points
vector[N] x; //
vector[N] y; //
real xpred; // input location for prediction
real pmualpha; // prior mean for alpha
real psalpha;  // prior std for alpha
real pmubeta;  // prior mean for beta
real psbeta;   // prior std for beta
}
parameters {
real alpha;
real beta;
real<lower=0> sigma;
}
transformed parameters {
vector[N] mu;
mu = alpha + beta*x;
}
model {
alpha ~ normal(pmualpha, psalpha);
beta ~ normal(pmubeta, psbeta);
y ~ normal(mu, sigma);
}
generated quantities {
real ypred;
vector[N] log_lik;
ypred = normal_rng(alpha + beta*xpred, sigma);
for (i in 1:N)
log_lik[i] = normal_lpdf(y[i] | mu[i], sigma);
}

``````

Create a list with data and priors:

``````

In [19]:

data = dict(
N = N,
x = x,
y = y,
xpred = xpred,
pmualpha = y.mean(),     # centered
psalpha  = 100,          # weakly informative prior
pmubeta  = 0,            # a priori increase and decrese as likely
psbeta   = (.1--.1)/6.0  # avg temp probably does not increase more than 1
# degree per 10 years
)

``````

Run Stan

``````

In [20]:

model = stan_utility.compile_model('lin.stan')
fit = model.sampling(data=data, seed=194838)
samples = fit.extract(permuted=True)

``````
``````

INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_6d95d5cbf9b9fb7d99a793b066a7bb72 NOW.
WARNING:pystan:51 of 4000 iterations saturated the maximum tree depth of 10 (1.275%)
WARNING:pystan:Run again with max_treedepth larger than 10 to avoid saturation

``````

Check the `n_eff` and `Rhat`

``````

In [21]:

print(fit)

``````
``````

Inference for Stan model: anon_model_6d95d5cbf9b9fb7d99a793b066a7bb72.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.

mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
alpha       -28.83    0.44  15.53 -58.72 -39.77 -28.66 -18.78   1.37   1251    1.0
beta          0.02  2.2e-4 7.8e-3 4.1e-3   0.01   0.02   0.02   0.03   1253    1.0
sigma         1.13  3.0e-3   0.11   0.94   1.05   1.12    1.2   1.35   1278    1.0
mu[1]         8.72  7.2e-3   0.27   8.18   8.54   8.72   8.91   9.26   1458    1.0
mu[2]         8.74  7.0e-3   0.27   8.22   8.56   8.74   8.93   9.27   1473    1.0
mu[3]         8.76  6.8e-3   0.26   8.25   8.59   8.76   8.94   9.27   1489    1.0
mu[4]         8.78  6.6e-3   0.26   8.28   8.61   8.78   8.95   9.28   1507    1.0
mu[5]          8.8  6.4e-3   0.25   8.31   8.63    8.8   8.97   9.29   1527    1.0
mu[6]         8.82  6.2e-3   0.24   8.34   8.66   8.82   8.99    9.3   1549    1.0
mu[7]         8.84  6.0e-3   0.24   8.37   8.68   8.84    9.0    9.3   1573    1.0
mu[8]         8.86  5.7e-3   0.23    8.4   8.71   8.86   9.02   9.31   1600    1.0
mu[9]         8.88  5.5e-3   0.22   8.43   8.73   8.88   9.03   9.32   1630    1.0
mu[10]         8.9  5.3e-3   0.22   8.47   8.75    8.9   9.05   9.32   1664    1.0
mu[11]        8.92  5.1e-3   0.21    8.5   8.77   8.92   9.06   9.33   1701    1.0
mu[12]        8.94  4.9e-3   0.21   8.53    8.8   8.94   9.08   9.34   1743    1.0
mu[13]        8.96  4.7e-3    0.2   8.56   8.82   8.96   9.09   9.34   1789    1.0
mu[14]        8.97  4.6e-3    0.2   8.59   8.84   8.98   9.11   9.36   1842    1.0
mu[15]        8.99  4.4e-3   0.19   8.61   8.87    9.0   9.12   9.36   1902    1.0
mu[16]        9.01  4.2e-3   0.19   8.64   8.89   9.01   9.14   9.37   1969    1.0
mu[17]        9.03  4.0e-3   0.18   8.68   8.91   9.03   9.16   9.38   2046    1.0
mu[18]        9.05  3.8e-3   0.18    8.7   8.93   9.05   9.17   9.39   2132    1.0
mu[19]        9.07  3.6e-3   0.17   8.73   8.96   9.07   9.19    9.4   2230    1.0
mu[20]        9.09  3.5e-3   0.17   8.76   8.98   9.09    9.2   9.41   2341    1.0
mu[21]        9.11  3.3e-3   0.16   8.78    9.0   9.11   9.22   9.43   2467    1.0
mu[22]        9.13  3.1e-3   0.16   8.81   9.02   9.13   9.24   9.44   2608    1.0
mu[23]        9.15  3.0e-3   0.16   8.84   9.04   9.15   9.26   9.45   2765    1.0
mu[24]        9.17  2.8e-3   0.15   8.87   9.06   9.17   9.27   9.46   2939    1.0
mu[25]        9.19  2.7e-3   0.15   8.89   9.09   9.18   9.29   9.48   3128    1.0
mu[26]        9.21  2.6e-3   0.15   8.92   9.11    9.2   9.31   9.49   3328    1.0
mu[27]        9.23  2.5e-3   0.15   8.94   9.13   9.22   9.33    9.5   3534    1.0
mu[28]        9.24  2.4e-3   0.15   8.96   9.15   9.24   9.34   9.52   3736    1.0
mu[29]        9.26  2.3e-3   0.14   8.99   9.17   9.26   9.36   9.54   3923    1.0
mu[30]        9.28  2.2e-3   0.14   9.01   9.18   9.28   9.38   9.56   4080    1.0
mu[31]         9.3  2.2e-3   0.14   9.03    9.2    9.3    9.4   9.58   4153    1.0
mu[32]        9.32  2.2e-3   0.14   9.05   9.22   9.32   9.42    9.6   4169    1.0
mu[33]        9.34  2.2e-3   0.14   9.06   9.24   9.34   9.44   9.62   4132    1.0
mu[34]        9.36  2.3e-3   0.15   9.08   9.26   9.36   9.46   9.65   4048    1.0
mu[35]        9.38  2.3e-3   0.15   9.09   9.28   9.38   9.48   9.67   3923    1.0
mu[36]         9.4  2.4e-3   0.15   9.11    9.3    9.4    9.5    9.7   3772    1.0
mu[37]        9.42  2.5e-3   0.15   9.12   9.32   9.42   9.52   9.72   3604    1.0
mu[38]        9.44  2.6e-3   0.15   9.14   9.33   9.44   9.54   9.75   3430    1.0
mu[39]        9.46  2.7e-3   0.16   9.15   9.35   9.45   9.56   9.77   3257    1.0
mu[40]        9.48  2.9e-3   0.16   9.17   9.37   9.48   9.58    9.8   3092    1.0
mu[41]        9.49  3.0e-3   0.16   9.18   9.39   9.49    9.6   9.82   2937    1.0
mu[42]        9.51  3.2e-3   0.17   9.19    9.4   9.51   9.62   9.85   2795    1.0
mu[43]        9.53  3.3e-3   0.17    9.2   9.42   9.53   9.65   9.87   2665    1.0
mu[44]        9.55  3.5e-3   0.18   9.21   9.43   9.55   9.67    9.9   2548    1.0
mu[45]        9.57  3.7e-3   0.18   9.21   9.45   9.57   9.69   9.93   2442    1.0
mu[46]        9.59  3.8e-3   0.19   9.22   9.47   9.59   9.71   9.96   2348    1.0
mu[47]        9.61  4.0e-3   0.19   9.23   9.48   9.61   9.74   9.99   2263    1.0
mu[48]        9.63  4.2e-3    0.2   9.24    9.5   9.63   9.76  10.02   2187    1.0
mu[49]        9.65  4.4e-3    0.2   9.24   9.51   9.65   9.78  10.05   2119    1.0
mu[50]        9.67  4.6e-3   0.21   9.25   9.53   9.67    9.8  10.08   2059    1.0
mu[51]        9.69  4.7e-3   0.21   9.26   9.54   9.69   9.83  10.11   2004    1.0
mu[52]        9.71  4.9e-3   0.22   9.27   9.56   9.71   9.85  10.15   1955    1.0
mu[53]        9.73  5.1e-3   0.22   9.28   9.57   9.73   9.87  10.18   1911    1.0
mu[54]        9.74  5.3e-3   0.23   9.28   9.59   9.75   9.89  10.21   1858    1.0
mu[55]        9.76  5.6e-3   0.24   9.29    9.6   9.76   9.92  10.24   1817    1.0
mu[56]        9.78  5.7e-3   0.24   9.29   9.62   9.78   9.94  10.27   1789    1.0
mu[57]         9.8  5.9e-3   0.25    9.3   9.63    9.8   9.97   10.3   1756    1.0
mu[58]        9.82  6.2e-3   0.26    9.3   9.65   9.82   9.99  10.33   1726    1.0
mu[59]        9.84  6.4e-3   0.26   9.31   9.66   9.84  10.01  10.36   1698    1.0
mu[60]        9.86  6.6e-3   0.27   9.32   9.68   9.86  10.04  10.39   1673    1.0
mu[61]        9.88  6.8e-3   0.28   9.33   9.69   9.88  10.06  10.42   1649    1.0
mu[62]         9.9  7.0e-3   0.28   9.33   9.71    9.9  10.08  10.46   1628    1.0
ypred         9.95    0.02   1.17   7.63   9.16   9.97   10.7  12.32   3579    1.0
log_lik[1]   -1.14  4.0e-3   0.13  -1.43  -1.21  -1.13  -1.05  -0.92   1030   1.01
log_lik[2]   -2.94    0.01   0.53  -4.15  -3.25  -2.88  -2.56  -2.05   1298   1.01
log_lik[3]   -1.23  3.9e-3   0.15  -1.58  -1.32  -1.21  -1.12  -0.98   1541    1.0
log_lik[4]   -1.25  4.4e-3   0.15  -1.59  -1.34  -1.24  -1.14  -0.99   1221   1.01
log_lik[5]   -1.26  4.4e-3   0.15   -1.6  -1.35  -1.25  -1.15   -1.0   1252   1.01
log_lik[6]   -1.56  5.5e-3   0.22  -2.04   -1.7  -1.54  -1.41   -1.2   1560    1.0
log_lik[7]   -1.08  3.2e-3    0.1   -1.3  -1.15  -1.08  -1.01  -0.88   1062   1.01
log_lik[8]   -1.08  2.7e-3    0.1  -1.29  -1.14  -1.08  -1.01  -0.89   1440    1.0
log_lik[9]   -2.87    0.01   0.46  -3.88  -3.14  -2.82  -2.55  -2.11   1314   1.01
log_lik[10]  -1.64  5.2e-3   0.21  -2.11  -1.77  -1.61  -1.48  -1.28   1701    1.0
log_lik[11]  -1.75  5.4e-3   0.23  -2.24  -1.89  -1.73  -1.59  -1.36   1733    1.0
log_lik[12]  -1.06  2.7e-3    0.1  -1.26  -1.12  -1.06  -0.99  -0.87   1268   1.01
log_lik[13]  -1.22  3.4e-3   0.12  -1.49   -1.3  -1.22  -1.14   -1.0   1341   1.01
log_lik[14]  -2.32  7.2e-3   0.31  -2.97  -2.51   -2.3   -2.1  -1.79   1802    1.0
log_lik[15]  -1.09  2.6e-3    0.1  -1.29  -1.15  -1.08  -1.02  -0.91   1476    1.0
log_lik[16]  -1.07  2.7e-3    0.1  -1.26  -1.13  -1.06   -1.0  -0.88   1292   1.01
log_lik[17]  -1.87  4.7e-3   0.21  -2.33  -2.01  -1.86  -1.72   -1.5   2053    1.0
log_lik[18]  -1.89  4.8e-3   0.22  -2.38  -2.03  -1.87  -1.74  -1.53   2024    1.0
log_lik[19]  -2.55  7.9e-3   0.33  -3.28  -2.75  -2.51  -2.32   -2.0   1708   1.01
log_lik[20]  -1.06  2.7e-3    0.1  -1.25  -1.12  -1.06   -1.0  -0.88   1302    1.0
log_lik[21]  -2.98  9.7e-3   0.39  -3.87  -3.22  -2.93   -2.7  -2.31   1667   1.01
log_lik[22]  -1.35  2.4e-3   0.12   -1.6  -1.43  -1.35  -1.27  -1.14   2388    1.0
log_lik[23]  -1.41  2.4e-3   0.12  -1.67  -1.49   -1.4  -1.32  -1.19   2626    1.0
log_lik[24]  -4.13    0.01   0.58  -5.39  -4.48  -4.09  -3.73  -3.11   1754    1.0
log_lik[25]  -1.44  2.2e-3   0.12  -1.69  -1.51  -1.43  -1.35  -1.22   2987    1.0
log_lik[26]  -1.31  2.1e-3   0.11  -1.52  -1.38   -1.3  -1.23  -1.11   2510    1.0
log_lik[27]  -1.07  2.5e-3   0.09  -1.26  -1.13  -1.07  -1.01   -0.9   1432    1.0
log_lik[28]  -1.22  2.1e-3    0.1  -1.42  -1.28  -1.21  -1.15  -1.03   2164    1.0
log_lik[29]  -1.76  2.7e-3   0.16  -2.11  -1.86  -1.75  -1.65  -1.49   3482    1.0
log_lik[30]  -2.18  4.3e-3   0.22  -2.68  -2.32  -2.17  -2.03   -1.8   2744    1.0
log_lik[31]  -2.08  3.8e-3   0.21  -2.53   -2.2  -2.06  -1.93  -1.73   2946    1.0
log_lik[32]  -1.64  2.2e-3   0.14  -1.95  -1.73  -1.63  -1.54   -1.4   3996    1.0
log_lik[33]   -1.4  1.9e-3   0.11  -1.63  -1.47  -1.39  -1.32   -1.2   3503    1.0
log_lik[34]  -1.09  2.4e-3   0.09  -1.28  -1.15  -1.09  -1.03  -0.91   1491    1.0
log_lik[35]  -1.05  2.6e-3   0.09  -1.24  -1.11  -1.05  -0.99  -0.87   1336    1.0
log_lik[36]  -2.82  7.5e-3   0.34  -3.55  -3.02   -2.8  -2.57  -2.23   2053    1.0
log_lik[37]  -1.36  2.0e-3   0.11  -1.59  -1.43  -1.35  -1.28  -1.15   3039    1.0
log_lik[38]  -1.06  2.6e-3   0.09  -1.24  -1.12  -1.05  -0.99  -0.88   1357    1.0
log_lik[39]  -1.33  2.2e-3   0.11  -1.57   -1.4  -1.33  -1.26  -1.12   2723    1.0
log_lik[40]  -1.09  2.5e-3    0.1  -1.28  -1.15  -1.09  -1.02  -0.91   1469    1.0
log_lik[41]  -1.14  2.4e-3    0.1  -1.35  -1.21  -1.14  -1.08  -0.95   1732    1.0
log_lik[42]  -1.12  2.5e-3    0.1  -1.32  -1.18  -1.11  -1.05  -0.93   1605    1.0
log_lik[43]  -1.05  2.6e-3   0.09  -1.23  -1.11  -1.05  -0.98  -0.87   1337    1.0
log_lik[44]  -1.34  2.7e-3   0.13  -1.62  -1.42  -1.33  -1.25  -1.11   2268    1.0
log_lik[45]  -1.09  2.6e-3    0.1  -1.29  -1.16  -1.09  -1.03  -0.91   1466    1.0
log_lik[46]  -1.38  2.9e-3   0.14  -1.69  -1.47  -1.37  -1.29  -1.14   2239    1.0
log_lik[47]  -1.07  2.6e-3    0.1  -1.26  -1.13  -1.07   -1.0  -0.88   1426    1.0
log_lik[48]  -1.21  2.8e-3   0.12  -1.47  -1.28   -1.2  -1.13   -1.0   1833    1.0
log_lik[49]  -1.22  2.9e-3   0.12  -1.49   -1.3  -1.21  -1.14   -1.0   1825    1.0
log_lik[50]  -1.05  2.6e-3    0.1  -1.25  -1.12  -1.05  -0.99  -0.87   1354    1.0
log_lik[51]  -2.24  7.3e-3   0.32  -2.93  -2.43  -2.21  -2.01  -1.69   1913    1.0
log_lik[52]  -1.45  4.0e-3   0.18  -1.85  -1.56  -1.44  -1.33  -1.16   1916    1.0
log_lik[53]  -1.11  2.9e-3   0.11  -1.34  -1.18   -1.1  -1.04  -0.92   1469    1.0
log_lik[54]  -1.51  4.6e-3    0.2  -1.94  -1.63  -1.48  -1.37  -1.18   1845    1.0
log_lik[55]  -1.22  3.4e-3   0.14  -1.54   -1.3  -1.21  -1.12  -0.98   1617    1.0
log_lik[56]  -1.17  3.3e-3   0.13  -1.46  -1.24  -1.16  -1.08  -0.95   1530    1.0
log_lik[57]  -1.46  5.0e-3    0.2  -1.93  -1.58  -1.44  -1.32  -1.13   1673    1.0
log_lik[58]  -1.06  2.7e-3    0.1  -1.27  -1.13  -1.06  -0.99  -0.87   1379    1.0
log_lik[59]   -1.5  5.5e-3   0.22   -2.0  -1.63  -1.47  -1.34  -1.14   1621    1.0
log_lik[60]  -1.42  5.2e-3   0.21   -1.9  -1.54  -1.39  -1.27  -1.09   1605    1.0
log_lik[61]  -1.72  7.2e-3   0.29  -2.37  -1.89  -1.69  -1.52  -1.26   1554    1.0
log_lik[62]  -1.65  6.8e-3   0.27  -2.27  -1.81  -1.62  -1.45  -1.21   1610    1.0
lp__        -38.06    0.04   1.28 -41.36 -38.61 -37.74 -37.14 -36.62   1078    1.0

Samples were drawn using NUTS at Tue Oct 23 13:01:01 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).

``````

Check the treedepth, E-BFMI, and divergences

``````

In [22]:

stan_utility.check_treedepth(fit)
stan_utility.check_energy(fit)
stan_utility.check_div(fit)

``````
``````

51 of 4000 iterations saturated the maximum tree depth of 10 (1.275%)
Run again with max_depth set to a larger value to avoid saturation
0.0 of 4000 iterations ended with a divergence (0.0%)

``````

We get a warning that several iterations that saturated the maximum treedepth". The reason for this is a very strong posterior correlation between alpha and beta as shown in the next plot. This doesn't imvalidate the results, but leads to suboptimal performance. We'll later look at alternative which reduces the posterior correlation. The following figure shows the strong posterior correlation in the current case.

``````

In [23]:

samples = fit.extract(permuted=True)
# preview 500 posterior samples of (alpha, beta)
plt.scatter(samples['alpha'][:500], samples['beta'][:500], 5)
plt.xlabel('alpha')
plt.ylabel('beta');

``````
``````

``````

Compute the probability that the summer temperature is increasing.

``````

In [24]:

print('Pr(beta > 0) = {}'.format(np.mean(samples['beta'] > 0)))

``````
``````

Pr(beta > 0) = 0.98975

``````

Plot the data, the model fit and prediction for year 2016.

``````

In [25]:

# make slightly wider figure of 3 plots
figsize = plt.rcParams['figure.figsize'].copy()
figsize[0] *= 1.4  # width
figsize[1] *= 2.5  # height
fig, axes = plt.subplots(3, 1, figsize=figsize)

# plot 1: scatterplot and lines
color_scatter = 'C0'  # 'C0' for default color #0
color_line = 'C1'     # 'C1' for default color #1
# lighten color_line
1 - 0.1*(1 - np.array(mpl.colors.to_rgb(color_line)))
)
# plot
ax = axes[0]
ax.fill_between(
x,
np.percentile(samples['mu'], 5, axis=0),
np.percentile(samples['mu'], 95, axis=0),
)
ax.plot(
x,
np.percentile(samples['mu'], 50, axis=0),
color=color_line,
linewidth=1
)
ax.scatter(x, y, 5, color=color_scatter)
ax.set_xlabel('year')
ax.set_ylabel('summer temperature at Kilpisjarvi')
ax.set_xlim((1952, 2013))

# plot 2: histogram
ax = axes[1]
ax.hist(samples['beta'], 50)
ax.set_xlabel('beta')

# plot 3: histogram
ax = axes[2]
ax.hist(samples['ypred'], 50)
ax.set_xlabel('y-prediction for x={}'.format(xpred))

# make figure compact
fig.tight_layout();

``````
``````

``````

### Gaussian linear model with standardized data

In the above we used the unnormalized data and as x values are far away from zero, this will cause very strong posterior dependency between alpha and beta (did you use ShinyStan for the above model?). The strong posterior dependency can be removed by normalizing the data to have zero mean. The following Stan code makes it in Stan. In generated quantities we do correspnding transformation back to the original scale.

``````

In [26]:

with open('lin_std.stan') as file:

``````
``````

// Gaussian linear model with standardized data
data {
int<lower=0> N; // number of data points
vector[N] x; //
vector[N] y; //
real xpred; // input location for prediction
}
transformed data {
vector[N] x_std;
vector[N] y_std;
real xpred_std;
x_std = (x - mean(x)) / sd(x);
y_std = (y - mean(y)) / sd(y);
xpred_std = (xpred - mean(x)) / sd(x);
}
parameters {
real alpha;
real beta;
real<lower=0> sigma_std;
}
transformed parameters {
vector[N] mu_std;
mu_std = alpha + beta*x_std;
}
model {
alpha ~ normal(0, 1);
beta ~ normal(0, 1);
y_std ~ normal(mu_std, sigma_std);
}
generated quantities {
vector[N] mu;
real<lower=0> sigma;
real ypred;
vector[N] log_lik;
mu = mu_std*sd(y) + mean(y);
sigma = sigma_std*sd(y);
ypred = normal_rng((alpha + beta*xpred_std)*sd(y)+mean(y), sigma*sd(y));
for (i in 1:N)
log_lik[i] = normal_lpdf(y[i] | mu[i], sigma);
}

``````

Run Stan

``````

In [27]:

model = stan_utility.compile_model('lin_std.stan')
fit = model.sampling(data=data, seed=194838)
samples = fit.extract(permuted=True)

``````
``````

INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_161eb30961f09398a0825482beb4b070 NOW.

``````

Check the `n_eff` and `Rhat`

``````

In [28]:

print(fit)

``````
``````

Inference for Stan model: anon_model_161eb30961f09398a0825482beb4b070.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.

mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
alpha       4.0e-4  2.1e-3   0.12  -0.24  -0.08 6.8e-5   0.08   0.25   3614    1.0
beta          0.31  2.1e-3   0.12   0.08   0.23   0.31   0.39   0.55   3471    1.0
sigma_std     0.97  1.5e-3   0.09   0.81   0.91   0.97   1.03   1.17   3560    1.0
mu_std[1]    -0.52  3.9e-3   0.23  -0.97  -0.68  -0.53  -0.37  -0.06   3525    1.0
mu_std[2]    -0.51  3.8e-3   0.23  -0.94  -0.66  -0.51  -0.36  -0.05   3528    1.0
mu_std[3]    -0.49  3.7e-3   0.22  -0.92  -0.64  -0.49  -0.35  -0.04   3532    1.0
mu_std[4]    -0.47  3.7e-3   0.22  -0.89  -0.62  -0.48  -0.33  -0.03   3536    1.0
mu_std[5]    -0.46  3.6e-3   0.21  -0.86   -0.6  -0.46  -0.32  -0.03   3540    1.0
mu_std[6]    -0.44  3.5e-3   0.21  -0.83  -0.58  -0.44   -0.3  -0.02   3544    1.0
mu_std[7]    -0.42  3.4e-3    0.2   -0.8  -0.56  -0.42  -0.29  -0.02   3549    1.0
mu_std[8]     -0.4  3.3e-3    0.2  -0.78  -0.54  -0.41  -0.28-7.0e-3   3554    1.0
mu_std[9]    -0.39  3.2e-3   0.19  -0.75  -0.51  -0.39  -0.26 1.6e-3   3559    1.0
mu_std[10]   -0.37  3.1e-3   0.19  -0.72  -0.49  -0.37  -0.25 8.5e-3   3565    1.0
mu_std[11]   -0.35  3.0e-3   0.18  -0.69  -0.47  -0.35  -0.24   0.01   3572    1.0
mu_std[12]   -0.34  2.9e-3   0.18  -0.67  -0.45  -0.34  -0.22   0.02   3578    1.0
mu_std[13]   -0.32  2.9e-3   0.17  -0.64  -0.43  -0.32  -0.21   0.03   3586    1.0
mu_std[14]    -0.3  2.8e-3   0.17  -0.62  -0.41   -0.3  -0.19   0.04   3593    1.0
mu_std[15]   -0.28  2.7e-3   0.16  -0.59  -0.39  -0.29  -0.18   0.05   3601    1.0
mu_std[16]   -0.27  2.6e-3   0.16  -0.56  -0.37  -0.27  -0.16   0.06   3610    1.0
mu_std[17]   -0.25  2.6e-3   0.15  -0.54  -0.35  -0.25  -0.15   0.06   3619    1.0
mu_std[18]   -0.23  2.5e-3   0.15  -0.52  -0.33  -0.23  -0.13   0.07   3629    1.0
mu_std[19]   -0.21  2.4e-3   0.15  -0.49  -0.31  -0.22  -0.12   0.08   3638    1.0
mu_std[20]    -0.2  2.4e-3   0.14  -0.47  -0.29   -0.2  -0.11   0.09   3633    1.0
mu_std[21]   -0.18  2.3e-3   0.14  -0.45  -0.27  -0.18  -0.09    0.1   3628    1.0
mu_std[22]   -0.16  2.3e-3   0.14  -0.43  -0.25  -0.16  -0.07   0.11   3623    1.0
mu_std[23]   -0.15  2.2e-3   0.13  -0.41  -0.23  -0.15  -0.06   0.13   3619    1.0
mu_std[24]   -0.13  2.2e-3   0.13  -0.38  -0.21  -0.13  -0.04   0.14   3616    1.0
mu_std[25]   -0.11  2.1e-3   0.13  -0.36   -0.2  -0.11  -0.03   0.15   3613    1.0
mu_std[26]   -0.09  2.1e-3   0.13  -0.34  -0.18   -0.1-9.9e-3   0.16   3610    1.0
mu_std[27]   -0.08  2.1e-3   0.13  -0.32  -0.16  -0.08 7.0e-3   0.18   3609    1.0
mu_std[28]   -0.06  2.1e-3   0.12   -0.3  -0.14  -0.06   0.02   0.19   3608    1.0
mu_std[29]   -0.04  2.1e-3   0.12  -0.28  -0.12  -0.04   0.04    0.2   3609    1.0
mu_std[30]   -0.03  2.1e-3   0.12  -0.27  -0.11  -0.03   0.06   0.22   3610    1.0
mu_std[31] -8.2e-3  2.0e-3   0.12  -0.25  -0.09-9.1e-3   0.07   0.24   3613    1.0
mu_std[32]  9.0e-3  2.1e-3   0.12  -0.24  -0.07 9.2e-3   0.09   0.26   3616    1.0
mu_std[33]    0.03  2.1e-3   0.12  -0.22  -0.06   0.03   0.11   0.28   3620    1.0
mu_std[34]    0.04  2.1e-3   0.13   -0.2  -0.04   0.04   0.13    0.3   3625    1.0
mu_std[35]    0.06  2.1e-3   0.13  -0.19  -0.02   0.06   0.14   0.31   3631    1.0
mu_std[36]    0.08  2.1e-3   0.13  -0.17-10.0e-3   0.08   0.16   0.33   3637    1.0
mu_std[37]     0.1  2.2e-3   0.13  -0.16 5.6e-3    0.1   0.18   0.36   3644    1.0
mu_std[38]    0.11  2.2e-3   0.13  -0.15   0.02   0.11    0.2   0.38   3651    1.0
mu_std[39]    0.13  2.2e-3   0.14  -0.14   0.04   0.13   0.22    0.4   3658    1.0
mu_std[40]    0.15  2.3e-3   0.14  -0.13   0.05   0.15   0.24   0.43   3665    1.0
mu_std[41]    0.16  2.3e-3   0.14  -0.11   0.07   0.16   0.26   0.45   3672    1.0
mu_std[42]    0.18  2.4e-3   0.14   -0.1   0.08   0.18   0.28   0.48   3673    1.0
mu_std[43]     0.2  2.5e-3   0.15   -0.1    0.1    0.2    0.3    0.5   3664    1.0
mu_std[44]    0.22  2.5e-3   0.15  -0.09   0.11   0.21   0.32   0.53   3655    1.0
mu_std[45]    0.23  2.6e-3   0.16  -0.08   0.13   0.23   0.34   0.55   3646    1.0
mu_std[46]    0.25  2.7e-3   0.16  -0.07   0.14   0.25   0.36   0.57   3638    1.0
mu_std[47]    0.27  2.7e-3   0.17  -0.06   0.16   0.26   0.38    0.6   3629    1.0
mu_std[48]    0.28  2.8e-3   0.17  -0.05   0.17   0.28    0.4   0.62   3622    1.0
mu_std[49]     0.3  2.9e-3   0.17  -0.04   0.18    0.3   0.42   0.65   3614    1.0
mu_std[50]    0.32  3.0e-3   0.18  -0.03    0.2   0.32   0.44   0.68   3607    1.0
mu_std[51]    0.34  3.1e-3   0.18  -0.02   0.21   0.33   0.46    0.7   3601    1.0
mu_std[52]    0.35  3.2e-3   0.19-9.1e-3   0.22   0.35   0.48   0.73   3594    1.0
mu_std[53]    0.37  3.2e-3   0.19-4.7e-3   0.24   0.37    0.5   0.75   3588    1.0
mu_std[54]    0.39  3.3e-3    0.2 1.3e-3   0.25   0.39   0.52   0.78   3583    1.0
mu_std[55]     0.4  3.4e-3    0.2 6.7e-3   0.26    0.4   0.54   0.81   3577    1.0
mu_std[56]    0.42  3.5e-3   0.21   0.01   0.28   0.42   0.56   0.83   3572    1.0
mu_std[57]    0.44  3.6e-3   0.22   0.02   0.29   0.44   0.58   0.86   3568    1.0
mu_std[58]    0.46  3.7e-3   0.22   0.03   0.31   0.45    0.6   0.89   3563    1.0
mu_std[59]    0.47  3.8e-3   0.23   0.03   0.32   0.47   0.63   0.92   3559    1.0
mu_std[60]    0.49  3.9e-3   0.23   0.04   0.34   0.49   0.65   0.94   3555    1.0
mu_std[61]    0.51  4.0e-3   0.24   0.05   0.35   0.51   0.67   0.97   3552    1.0
mu_std[62]    0.53  4.1e-3   0.24   0.05   0.36   0.52   0.69    1.0   3548    1.0
mu[1]          8.7  4.6e-3   0.27   8.18   8.52    8.7   8.88   9.25   3525    1.0
mu[2]         8.72  4.5e-3   0.26   8.22   8.55   8.72    8.9   9.26   3528    1.0
mu[3]         8.74  4.3e-3   0.26   8.25   8.57   8.74   8.91   9.27   3532    1.0
mu[4]         8.76  4.2e-3   0.25   8.28    8.6   8.76   8.93   9.28   3536    1.0
mu[5]         8.78  4.1e-3   0.25   8.32   8.62   8.78   8.95   9.28   3540    1.0
mu[6]          8.8  4.0e-3   0.24   8.35   8.64    8.8   8.96   9.29   3544    1.0
mu[7]         8.82  3.9e-3   0.23   8.38   8.67   8.82   8.98    9.3   3549    1.0
mu[8]         8.84  3.8e-3   0.23   8.41   8.69   8.84   8.99    9.3   3554    1.0
mu[9]         8.86  3.7e-3   0.22   8.44   8.72   8.86   9.01   9.31   3559    1.0
mu[10]        8.88  3.6e-3   0.21   8.48   8.74   8.88   9.02   9.32   3565    1.0
mu[11]         8.9  3.5e-3   0.21   8.51   8.76    8.9   9.04   9.33   3572    1.0
mu[12]        8.92  3.4e-3    0.2   8.54   8.79   8.92   9.06   9.34   3578    1.0
mu[13]        8.94  3.3e-3    0.2   8.57   8.81   8.94   9.07   9.35   3586    1.0
mu[14]        8.96  3.2e-3   0.19    8.6   8.84   8.96   9.09   9.36   3593    1.0
mu[15]        8.98  3.1e-3   0.19   8.63   8.86   8.98   9.11   9.37   3601    1.0
mu[16]         9.0  3.0e-3   0.18   8.66   8.88    9.0   9.12   9.38   3610    1.0
mu[17]        9.02  3.0e-3   0.18   8.69   8.91   9.02   9.14   9.39   3619    1.0
mu[18]        9.04  2.9e-3   0.17   8.71   8.93   9.04   9.16    9.4   3629    1.0
mu[19]        9.06  2.8e-3   0.17   8.74   8.95   9.06   9.17   9.41   3638    1.0
mu[20]        9.08  2.7e-3   0.17   8.77   8.97   9.08   9.19   9.42   3633    1.0
mu[21]         9.1  2.7e-3   0.16   8.79    9.0    9.1   9.21   9.43   3628    1.0
mu[22]        9.12  2.6e-3   0.16   8.82   9.02   9.12   9.23   9.45   3623    1.0
mu[23]        9.14  2.6e-3   0.15   8.84   9.04   9.14   9.25   9.46   3619    1.0
mu[24]        9.16  2.5e-3   0.15   8.87   9.06   9.16   9.26   9.47   3616    1.0
mu[25]        9.18  2.5e-3   0.15   8.89   9.09   9.18   9.28   9.49   3613    1.0
mu[26]         9.2  2.5e-3   0.15   8.92   9.11    9.2    9.3    9.5   3610    1.0
mu[27]        9.22  2.4e-3   0.15   8.94   9.13   9.22   9.32   9.52   3609    1.0
mu[28]        9.24  2.4e-3   0.14   8.96   9.15   9.24   9.34   9.53   3608    1.0
mu[29]        9.26  2.4e-3   0.14   8.98   9.17   9.26   9.36   9.55   3609    1.0
mu[30]        9.28  2.4e-3   0.14    9.0   9.19   9.28   9.38   9.57   3610    1.0
mu[31]         9.3  2.4e-3   0.14   9.02   9.21    9.3    9.4   9.59   3613    1.0
mu[32]        9.32  2.4e-3   0.14   9.04   9.23   9.32   9.42   9.61   3616    1.0
mu[33]        9.34  2.4e-3   0.14   9.06   9.25   9.34   9.44   9.63   3620    1.0
mu[34]        9.36  2.4e-3   0.15   9.08   9.27   9.36   9.46   9.66   3625    1.0
mu[35]        9.38  2.4e-3   0.15    9.1   9.28   9.38   9.48   9.68   3631    1.0
mu[36]         9.4  2.5e-3   0.15   9.11    9.3    9.4    9.5    9.7   3637    1.0
mu[37]        9.42  2.5e-3   0.15   9.13   9.32   9.42   9.52   9.73   3644    1.0
mu[38]        9.44  2.5e-3   0.15   9.14   9.34   9.44   9.54   9.75   3651    1.0
mu[39]        9.46  2.6e-3   0.16   9.15   9.36   9.46   9.57   9.78   3658    1.0
mu[40]        9.48  2.6e-3   0.16   9.17   9.38   9.48   9.59   9.81   3665    1.0
mu[41]         9.5  2.7e-3   0.16   9.18   9.39    9.5   9.61   9.84   3672    1.0
mu[42]        9.52  2.8e-3   0.17   9.19   9.41   9.52   9.64   9.87   3673    1.0
mu[43]        9.54  2.8e-3   0.17    9.2   9.43   9.54   9.66   9.89   3664    1.0
mu[44]        9.56  2.9e-3   0.18   9.21   9.44   9.56   9.68   9.92   3655    1.0
mu[45]        9.58  3.0e-3   0.18   9.23   9.46   9.58    9.7   9.95   3646    1.0
mu[46]         9.6  3.1e-3   0.19   9.24   9.48    9.6   9.73   9.98   3638    1.0
mu[47]        9.62  3.2e-3   0.19   9.25   9.49   9.62   9.75   10.0   3629    1.0
mu[48]        9.64  3.3e-3    0.2   9.26   9.51   9.64   9.77  10.04   3622    1.0
mu[49]        9.66  3.4e-3    0.2   9.27   9.53   9.66    9.8  10.07   3614    1.0
mu[50]        9.68  3.5e-3   0.21   9.28   9.54   9.68   9.82   10.1   3607    1.0
mu[51]         9.7  3.6e-3   0.21   9.29   9.56    9.7   9.85  10.13   3601    1.0
mu[52]        9.72  3.7e-3   0.22    9.3   9.57   9.72   9.87  10.16   3594    1.0
mu[53]        9.74  3.8e-3   0.23   9.31   9.59   9.74   9.89  10.19   3588    1.0
mu[54]        9.76  3.9e-3   0.23   9.31    9.6   9.76   9.92  10.22   3583    1.0
mu[55]        9.78  4.0e-3   0.24   9.32   9.62   9.78   9.94  10.25   3577    1.0
mu[56]         9.8  4.1e-3   0.24   9.33   9.64    9.8   9.96  10.28   3572    1.0
mu[57]        9.82  4.2e-3   0.25   9.34   9.65   9.82   9.99  10.31   3568    1.0
mu[58]        9.84  4.3e-3   0.26   9.34   9.67   9.84  10.01  10.34   3563    1.0
mu[59]        9.86  4.4e-3   0.26   9.35   9.68   9.86  10.04  10.38   3559    1.0
mu[60]        9.88  4.5e-3   0.27   9.36    9.7   9.88  10.06  10.41   3555    1.0
mu[61]         9.9  4.6e-3   0.28   9.37   9.72    9.9  10.09  10.44   3552    1.0
mu[62]        9.92  4.7e-3   0.28   9.37   9.73   9.92  10.11  10.48   3548    1.0
sigma         1.13  1.8e-3    0.1   0.94   1.06   1.12   1.19   1.36   3560    1.0
ypred        10.01    0.02   1.32   7.43   9.15  10.01  10.89  12.66   4046    1.0
log_lik[1]   -1.13  2.3e-3   0.13  -1.42  -1.21  -1.12  -1.04  -0.92   3230    1.0
log_lik[2]   -2.97  8.7e-3   0.53   -4.1  -3.31  -2.94  -2.59  -2.06   3757    1.0
log_lik[3]   -1.24  2.6e-3   0.15  -1.57  -1.32  -1.22  -1.13  -0.99   3323    1.0
log_lik[4]   -1.24  2.6e-3   0.15  -1.59  -1.33  -1.22  -1.13  -0.99   3546    1.0
log_lik[5]   -1.25  2.5e-3   0.15   -1.6  -1.34  -1.23  -1.14   -1.0   3573    1.0
log_lik[6]   -1.55  3.5e-3   0.21  -2.05  -1.68  -1.52  -1.39  -1.19   3762    1.0
log_lik[7]   -1.08  1.9e-3    0.1  -1.29  -1.14  -1.07   -1.0  -0.88   3122    1.0
log_lik[8]   -1.08  1.8e-3    0.1  -1.29  -1.15  -1.08  -1.01  -0.89   3001    1.0
log_lik[9]    -2.9  7.4e-3   0.45  -3.86  -3.18  -2.86  -2.57  -2.12   3743    1.0
log_lik[10]  -1.65  3.4e-3   0.21  -2.12  -1.78  -1.63   -1.5  -1.29   3805    1.0
log_lik[11]  -1.74  3.6e-3   0.22  -2.25  -1.87  -1.71  -1.58  -1.36   3813    1.0
log_lik[12]  -1.06  1.7e-3    0.1  -1.26  -1.12  -1.05  -0.99  -0.87   3150    1.0
log_lik[13]  -1.22  2.1e-3   0.12  -1.49  -1.29  -1.21  -1.13   -1.0   3558    1.0
log_lik[14]   -2.3  4.9e-3    0.3  -2.96  -2.49  -2.27  -2.09   -1.8   3756    1.0
log_lik[15]  -1.09  1.7e-3    0.1  -1.29  -1.15  -1.09  -1.02   -0.9   3182    1.0
log_lik[16]  -1.06  1.7e-3    0.1  -1.26  -1.13  -1.06   -1.0  -0.88   3283    1.0
log_lik[17]  -1.86  3.4e-3   0.21  -2.33  -1.99  -1.84  -1.72  -1.51   3826    1.0
log_lik[18]   -1.9  3.4e-3   0.21  -2.36  -2.04  -1.89  -1.75  -1.52   3786    1.0
log_lik[19]  -2.56  5.3e-3   0.32  -3.26  -2.76  -2.54  -2.33  -1.99   3741    1.0
log_lik[20]  -1.06  1.6e-3   0.09  -1.25  -1.12  -1.05  -0.99  -0.88   3328    1.0
log_lik[21]  -2.99  6.4e-3   0.39  -3.85  -3.24  -2.97  -2.71  -2.31   3723    1.0
log_lik[22]  -1.35  1.9e-3   0.12   -1.6  -1.43  -1.35  -1.27  -1.15   3739    1.0
log_lik[23]  -1.41  2.0e-3   0.12  -1.67  -1.49  -1.41  -1.33   -1.2   3763    1.0
log_lik[24]  -4.12  9.6e-3   0.58  -5.41  -4.48  -4.08  -3.72  -3.14   3590    1.0
log_lik[25]  -1.43  2.0e-3   0.12  -1.69  -1.51  -1.42  -1.35  -1.22   3674    1.0
log_lik[26]   -1.3  1.7e-3   0.11  -1.53  -1.37   -1.3  -1.23  -1.11   3625    1.0
log_lik[27]  -1.07  1.6e-3   0.09  -1.27  -1.14  -1.07  -1.01  -0.89   3349    1.0
log_lik[28]  -1.22  1.6e-3    0.1  -1.42  -1.28  -1.21  -1.15  -1.04   3582    1.0
log_lik[29]  -1.76  2.6e-3   0.16   -2.1  -1.86  -1.75  -1.65  -1.48   3787    1.0
log_lik[30]  -2.19  3.7e-3   0.23  -2.68  -2.33  -2.17  -2.02   -1.8   3683    1.0
log_lik[31]  -2.08  3.4e-3   0.21  -2.53  -2.21  -2.06  -1.93  -1.72   3675    1.0
log_lik[32]  -1.65  2.3e-3   0.14  -1.95  -1.73  -1.64  -1.55  -1.39   3671    1.0
log_lik[33]   -1.4  1.8e-3   0.11  -1.64  -1.47   -1.4  -1.32   -1.2   3658    1.0
log_lik[34]  -1.09  1.6e-3   0.09  -1.28  -1.15  -1.09  -1.03  -0.91   3414    1.0
log_lik[35]  -1.05  1.6e-3   0.09  -1.24  -1.11  -1.05  -0.99  -0.87   3355    1.0
log_lik[36]  -2.83  5.8e-3   0.35  -3.61  -3.04   -2.8  -2.58  -2.24   3634    1.0
log_lik[37]  -1.35  1.8e-3   0.11  -1.59  -1.43  -1.35  -1.27  -1.15   3705    1.0
log_lik[38]  -1.05  1.6e-3   0.09  -1.25  -1.12  -1.05  -0.99  -0.87   3354    1.0
log_lik[39]  -1.33  1.9e-3   0.11  -1.57   -1.4  -1.32  -1.25  -1.12   3694    1.0
log_lik[40]  -1.09  1.6e-3    0.1  -1.28  -1.15  -1.08  -1.02   -0.9   3410    1.0
log_lik[41]  -1.15  1.7e-3    0.1  -1.36  -1.21  -1.14  -1.08  -0.97   3360    1.0
log_lik[42]  -1.12  1.7e-3    0.1  -1.32  -1.18  -1.11  -1.05  -0.94   3311    1.0
log_lik[43]  -1.05  1.6e-3   0.09  -1.24  -1.11  -1.04  -0.98  -0.87   3260    1.0
log_lik[44]  -1.35  2.2e-3   0.13  -1.63  -1.43  -1.34  -1.26  -1.12   3472    1.0
log_lik[45]  -1.09  1.7e-3    0.1   -1.3  -1.15  -1.08  -1.02   -0.9   3390    1.0
log_lik[46]  -1.37  2.3e-3   0.14  -1.67  -1.46  -1.37  -1.27  -1.13   3754    1.0
log_lik[47]  -1.07  1.7e-3    0.1  -1.27  -1.13  -1.06   -1.0  -0.89   3167    1.0
log_lik[48]  -1.22  2.1e-3   0.12  -1.48  -1.29  -1.21  -1.13  -1.01   3294    1.0
log_lik[49]  -1.23  2.1e-3   0.12   -1.5   -1.3  -1.22  -1.14  -1.02   3289    1.0
log_lik[50]  -1.05  1.7e-3    0.1  -1.25  -1.11  -1.05  -0.99  -0.87   3090    1.0
log_lik[51]  -2.21  5.0e-3   0.31  -2.86  -2.41  -2.19  -1.99  -1.68   3751    1.0
log_lik[52]  -1.44  2.8e-3   0.17  -1.81  -1.55  -1.42  -1.31  -1.14   3744    1.0
log_lik[53]  -1.11  1.9e-3   0.11  -1.34  -1.17   -1.1  -1.03   -0.9   3360    1.0
log_lik[54]  -1.49  3.1e-3   0.19  -1.91  -1.61  -1.47  -1.35  -1.16   3745    1.0
log_lik[55]  -1.21  2.3e-3   0.14  -1.51  -1.29   -1.2  -1.11  -0.97   3571    1.0
log_lik[56]  -1.16  2.2e-3   0.13  -1.43  -1.23  -1.15  -1.07  -0.94   3470    1.0
log_lik[57]  -1.48  3.6e-3   0.21  -1.96   -1.6  -1.46  -1.33  -1.15   3393    1.0
log_lik[58]  -1.06  1.9e-3    0.1  -1.27  -1.12  -1.06  -0.99  -0.87   2818    1.0
log_lik[59]  -1.52  3.9e-3   0.23  -2.04  -1.65  -1.49  -1.35  -1.16   3397    1.0
log_lik[60]   -1.4  3.3e-3    0.2  -1.85  -1.52  -1.38  -1.26  -1.07   3692    1.0
log_lik[61]  -1.75  5.0e-3   0.29  -2.43  -1.93  -1.71  -1.54  -1.29   3468    1.0
log_lik[62]  -1.62  4.3e-3   0.26  -2.21  -1.79   -1.6  -1.43  -1.19   3726    1.0
lp__        -28.78    0.03   1.24  -31.9 -29.35 -28.44 -27.88 -27.39   1691    1.0

Samples were drawn using NUTS at Tue Oct 23 13:01:59 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).

``````

We get now much better effective sample sizes `n_eff`.

Check the treedepth, E-BFMI, and divergences

``````

In [29]:

stan_utility.check_treedepth(fit)
stan_utility.check_energy(fit)
stan_utility.check_div(fit)

``````
``````

0 of 4000 iterations saturated the maximum tree depth of 10 (0.0%)
0.0 of 4000 iterations ended with a divergence (0.0%)

``````

Everything is fine now. The next figure shows that with the standardized data there is not much posterior correlation:

``````

In [30]:

samples = fit.extract(permuted=True)
plt.figure()
# preview 2000 samples
plt.scatter(samples['alpha'][:2000], samples['beta'][:2000], 5)
plt.xlabel('alpha')
plt.ylabel('beta');

``````
``````

``````

Compute the probability that the summer temperature is increasing.

``````

In [31]:

print('Pr(beta > 0) = {}'.format(np.mean(samples['beta'] > 0)))

``````
``````

Pr(beta > 0) = 0.9915

``````

Plot the data, the model fit and prediction for year 2016.

``````

In [32]:

# make slightly wider figure of 3 plots
figsize = plt.rcParams['figure.figsize'].copy()
figsize[0] *= 1.4  # width
figsize[1] *= 2.5  # height
fig, axes = plt.subplots(3, 1, figsize=figsize)

# plot 1: scatterplot and lines
color_scatter = 'C0'  # 'C0' for default color #0
color_line = 'C1'     # 'C1' for default color #1
# lighten color_line
1 - 0.1*(1 - np.array(mpl.colors.to_rgb(color_line)))
)
# plot
ax = axes[0]
ax.fill_between(
x,
np.percentile(samples['mu'], 5, axis=0),
np.percentile(samples['mu'], 95, axis=0),
)
ax.plot(
x,
np.percentile(samples['mu'], 50, axis=0),
color=color_line,
linewidth=1
)
ax.scatter(x, y, 5, color=color_scatter)
ax.set_xlabel('year')
ax.set_ylabel('summer temperature at Kilpisjarvi')
ax.set_xlim((1952, 2013))

# plot 2: histogram
ax = axes[1]
ax.hist(samples['beta'], 50)
ax.set_xlabel('beta')

# plot 3: histogram
ax = axes[2]
ax.hist(samples['ypred'], 50)
ax.set_xlabel('y-prediction for x={}'.format(xpred))

# make figure compact
fig.tight_layout();

``````
``````

``````

## Linear Student's t model.

The temperatures used in the above analyses are averages over three months, which makes it more likely that they are normally distributed, but there can be extreme events in the feather and we can check whether more robust Student's t observation model woul give different results.

``````

In [33]:

with open('lin_t.stan') as file:

``````
``````

// Linear student-t model
data {
int<lower=0> N; // number of data points
vector[N] x; //
vector[N] y; //
real xpred; // input location for prediction
}
parameters {
real alpha;
real beta;
real<lower=0> sigma;
real<lower=1, upper=80> nu;
}
transformed parameters {
vector[N] mu;
mu = alpha + beta*x;
}
model {
nu ~ gamma(2, 0.1); // Juarez and Steel(2010)
y ~ student_t(nu, mu, sigma);
}
generated quantities {
real ypred;
vector[N] log_lik;
ypred = normal_rng(alpha + beta*xpred, sigma);
for (i in 1:N)
log_lik[i] = student_t_lpdf(y[i] | nu, mu[i], sigma);
}

``````
``````

In [34]:

model = stan_utility.compile_model('lin_t.stan')
fit = model.sampling(data=data, seed=194838)
samples = fit.extract(permuted=True)

``````
``````

INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_985e43a8fe28088e3317d00edd324bab NOW.
WARNING:pystan:145 of 4000 iterations saturated the maximum tree depth of 10 (3.625%)
WARNING:pystan:Run again with max_treedepth larger than 10 to avoid saturation

``````

Check the `n_eff` and `Rhat`

``````

In [35]:

print(fit)

``````
``````

Inference for Stan model: anon_model_985e43a8fe28088e3317d00edd324bab.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.

mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
alpha       -33.48    0.47  15.83 -64.89 -44.09 -33.68 -23.02  -2.18   1127    1.0
beta          0.02  2.4e-4 8.0e-3 5.8e-3   0.02   0.02   0.03   0.04   1128    1.0
sigma         1.09  2.7e-3   0.11   0.89   1.01   1.08   1.16   1.33   1707    1.0
nu            24.2    0.38  13.79   6.26  13.96  21.04  31.31  58.93   1347    1.0
mu[1]         8.65  7.9e-3   0.29    8.1   8.46   8.65   8.84   9.23   1357    1.0
mu[2]         8.67  7.6e-3   0.28   8.13   8.48   8.67   8.86   9.24   1372    1.0
mu[3]         8.69  7.4e-3   0.28   8.16   8.51   8.69   8.88   9.25   1389    1.0
mu[4]         8.72  7.2e-3   0.27    8.2   8.54   8.71   8.89   9.25   1408    1.0
mu[5]         8.74  6.9e-3   0.26   8.23   8.57   8.73   8.91   9.26   1428    1.0
mu[6]         8.76  6.7e-3   0.26   8.27   8.59   8.76   8.93   9.27   1451    1.0
mu[7]         8.78  6.5e-3   0.25    8.3   8.62   8.78   8.95   9.28   1476    1.0
mu[8]          8.8  6.3e-3   0.24   8.33   8.64    8.8   8.96   9.29   1503    1.0
mu[9]         8.82  6.0e-3   0.24   8.37   8.67   8.82   8.98    9.3   1534    1.0
mu[10]        8.85  5.8e-3   0.23    8.4    8.7   8.84    9.0   9.31   1567    1.0
mu[11]        8.87  5.6e-3   0.22   8.44   8.72   8.86   9.02   9.32   1605    1.0
mu[12]        8.89  5.4e-3   0.22   8.47   8.75   8.88   9.03   9.33   1647    1.0
mu[13]        8.91  5.2e-3   0.21    8.5   8.77   8.91   9.05   9.33   1695    1.0
mu[14]        8.93  5.0e-3   0.21   8.54    8.8   8.93   9.07   9.34   1748    1.0
mu[15]        8.95  4.7e-3    0.2   8.57   8.82   8.95   9.09   9.36   1808    1.0
mu[16]        8.98  4.5e-3    0.2   8.59   8.84   8.97    9.1   9.37   1876    1.0
mu[17]         9.0  4.3e-3   0.19   8.62   8.87   8.99   9.12   9.37   1953    1.0
mu[18]        9.02  4.1e-3   0.19   8.65   8.89   9.02   9.14   9.39   2040    1.0
mu[19]        9.04  3.9e-3   0.18   8.69   8.92   9.04   9.16    9.4   2140    1.0
mu[20]        9.06  3.7e-3   0.18   8.72   8.94   9.06   9.18   9.41   2257    1.0
mu[21]        9.08  3.5e-3   0.17   8.75   8.97   9.08    9.2   9.43   2391    1.0
mu[22]        9.11  3.3e-3   0.17   8.77   8.99    9.1   9.22   9.44   2538    1.0
mu[23]        9.13  3.2e-3   0.17    8.8   9.01   9.13   9.24   9.46   2706    1.0
mu[24]        9.15  3.0e-3   0.16   8.83   9.04   9.15   9.26   9.47   2893    1.0
mu[25]        9.17  2.8e-3   0.16   8.86   9.06   9.17   9.28   9.48   3101    1.0
mu[26]        9.19  2.7e-3   0.16   8.88   9.08   9.19    9.3    9.5   3327    1.0
mu[27]        9.21  2.6e-3   0.15   8.91   9.11   9.21   9.32   9.51   3565    1.0
mu[28]        9.23  2.5e-3   0.15   8.93   9.13   9.23   9.34   9.53   3772    1.0
mu[29]        9.26  2.4e-3   0.15   8.96   9.15   9.26   9.36   9.55   3909    1.0
mu[30]        9.28  2.3e-3   0.15   8.99   9.18   9.28   9.38   9.56   4021    1.0
mu[31]         9.3  2.3e-3   0.15   9.01    9.2    9.3    9.4   9.59   4097    1.0
mu[32]        9.32  2.3e-3   0.15   9.04   9.22   9.32   9.42   9.61   4131    1.0
mu[33]        9.34  2.3e-3   0.15   9.06   9.24   9.34   9.44   9.63   4121    1.0
mu[34]        9.36  2.3e-3   0.15   9.08   9.26   9.36   9.46   9.66   4066    1.0
mu[35]        9.39  2.4e-3   0.15   9.09   9.28   9.39   9.49   9.68   3972    1.0
mu[36]        9.41  2.4e-3   0.15   9.12    9.3   9.41   9.51   9.71   3847    1.0
mu[37]        9.43  2.5e-3   0.15   9.13   9.32   9.43   9.53   9.73   3702    1.0
mu[38]        9.45  2.6e-3   0.15   9.15   9.35   9.45   9.55   9.76   3545    1.0
mu[39]        9.47  2.7e-3   0.16   9.17   9.37   9.47   9.58   9.79   3383    1.0
mu[40]        9.49  2.8e-3   0.16   9.18   9.39   9.49    9.6   9.82   3225    1.0
mu[41]        9.52  3.0e-3   0.16    9.2    9.4   9.51   9.62   9.84   3073    1.0
mu[42]        9.54  3.1e-3   0.17   9.21   9.42   9.54   9.65   9.87   2931    1.0
mu[43]        9.56  3.2e-3   0.17   9.23   9.44   9.56   9.67    9.9   2799    1.0
mu[44]        9.58  3.4e-3   0.18   9.23   9.46   9.58    9.7   9.93   2679    1.0
mu[45]         9.6  3.5e-3   0.18   9.25   9.48    9.6   9.72   9.96   2568    1.0
mu[46]        9.62  3.7e-3   0.18   9.26    9.5   9.62   9.74   9.99   2458    1.0
mu[47]        9.64  3.9e-3   0.19   9.27   9.52   9.64   9.77  10.02   2359    1.0
mu[48]        9.67  4.1e-3   0.19   9.28   9.54   9.66   9.79  10.05   2271    1.0
mu[49]        9.69  4.3e-3    0.2   9.29   9.56   9.68   9.82  10.08   2168    1.0
mu[50]        9.71  4.5e-3   0.21    9.3   9.58   9.71   9.84  10.12   2076    1.0
mu[51]        9.73  4.7e-3   0.21   9.31    9.6   9.73   9.87  10.15   1994    1.0
mu[52]        9.75  4.9e-3   0.22   9.33   9.61   9.75   9.89  10.18   1922    1.0
mu[53]        9.77  5.2e-3   0.22   9.33   9.63   9.77   9.92  10.21   1857    1.0
mu[54]         9.8  5.4e-3   0.23   9.35   9.65   9.79   9.94  10.24   1800    1.0
mu[55]        9.82  5.6e-3   0.23   9.36   9.67   9.81   9.97  10.28   1748    1.0
mu[56]        9.84  5.8e-3   0.24   9.37   9.69   9.83   10.0  10.31   1702    1.0
mu[57]        9.86  6.1e-3   0.25   9.38    9.7   9.86  10.02  10.35   1663    1.0
mu[58]        9.88  6.3e-3   0.25   9.39   9.72   9.88  10.05  10.39   1628    1.0
mu[59]         9.9  6.5e-3   0.26   9.39   9.74    9.9  10.07  10.42   1596    1.0
mu[60]        9.93  6.7e-3   0.27    9.4   9.75   9.92   10.1  10.45   1567    1.0
mu[61]        9.95  7.0e-3   0.27   9.41   9.77   9.94  10.13  10.48   1541    1.0
mu[62]        9.97  7.2e-3   0.28   9.42   9.79   9.97  10.15  10.52   1517    1.0
ypred        10.02    0.02   1.12   7.74   9.28  10.01  10.76  12.22   3986    1.0
log_lik[1]    -1.1  3.7e-3   0.14  -1.41  -1.18  -1.09  -1.01  -0.87   1401    1.0
log_lik[2]   -3.07    0.01   0.53  -4.19   -3.4  -3.03   -2.7  -2.14   1398    1.0
log_lik[3]   -1.27  4.9e-3   0.19   -1.7  -1.38  -1.25  -1.14  -0.97   1481    1.0
log_lik[4]   -1.21  4.6e-3   0.17   -1.6  -1.31   -1.2   -1.1  -0.94   1363    1.0
log_lik[5]   -1.22  4.6e-3   0.17  -1.61  -1.32  -1.21  -1.11  -0.94   1379    1.0
log_lik[6]   -1.54  6.1e-3   0.24  -2.08  -1.68  -1.51  -1.37  -1.14   1485    1.0
log_lik[7]   -1.05  2.8e-3   0.11  -1.29  -1.12  -1.05  -0.98  -0.84   1630    1.0
log_lik[8]   -1.08  2.6e-3   0.12  -1.32  -1.15  -1.07   -1.0  -0.86   1986    1.0
log_lik[9]   -2.98    0.01   0.45  -3.94  -3.26  -2.95  -2.66  -2.17   1499    1.0
log_lik[10]  -1.73  6.7e-3   0.25  -2.29  -1.88  -1.71  -1.55   -1.3   1419    1.0
log_lik[11]  -1.74  5.8e-3   0.24   -2.3  -1.89  -1.72  -1.57  -1.33   1750    1.0
log_lik[12]  -1.04  2.5e-3   0.11  -1.25   -1.1  -1.03  -0.96  -0.84   1778    1.0
log_lik[13]   -1.2  3.5e-3   0.14   -1.5  -1.28  -1.19   -1.1  -0.96   1571    1.0
log_lik[14]  -2.31  6.9e-3   0.31  -3.02  -2.51  -2.29  -2.08  -1.78   2101    1.0
log_lik[15]  -1.08  2.4e-3   0.11  -1.31  -1.16  -1.08  -1.01  -0.87   2175    1.0
log_lik[16]  -1.04  2.6e-3   0.11  -1.26  -1.11  -1.04  -0.97  -0.84   1720    1.0
log_lik[17]  -1.88  4.9e-3   0.23   -2.4  -2.03  -1.86  -1.71  -1.48   2221    1.0
log_lik[18]  -1.98  6.0e-3   0.24  -2.52  -2.13  -1.96  -1.81  -1.56   1669    1.0
log_lik[19]  -2.64  8.0e-3   0.34  -3.34  -2.85  -2.62   -2.4  -2.04   1753    1.0
log_lik[20]  -1.04  2.5e-3    0.1  -1.25   -1.1  -1.03  -0.97  -0.84   1758    1.0
log_lik[21]  -3.04  8.8e-3   0.38  -3.84  -3.29  -3.01  -2.77  -2.35   1910    1.0
log_lik[22]  -1.39  2.8e-3   0.14  -1.68  -1.47  -1.38  -1.29  -1.14   2477    1.0
log_lik[23]  -1.45  2.8e-3   0.14  -1.75  -1.54  -1.44  -1.35   -1.2   2510    1.0
log_lik[24]  -3.97  9.5e-3    0.5  -5.05  -4.28  -3.93  -3.62  -3.09   2781    1.0
log_lik[25]  -1.44  2.5e-3   0.14  -1.73  -1.53  -1.44  -1.35   -1.2   3002    1.0
log_lik[26]  -1.31  2.3e-3   0.12  -1.55  -1.38   -1.3  -1.22  -1.09   2730    1.0
log_lik[27]  -1.06  2.2e-3    0.1  -1.26  -1.13  -1.06  -0.99  -0.86   2217    1.0
log_lik[28]  -1.22  1.9e-3   0.11  -1.45   -1.3  -1.22  -1.15  -1.02   3207    1.0
log_lik[29]  -1.82  3.3e-3   0.18   -2.2  -1.93   -1.8  -1.69  -1.51   2880    1.0
log_lik[30]  -2.23  4.2e-3   0.24  -2.75  -2.38  -2.21  -2.06  -1.82   3234    1.0
log_lik[31]  -2.13  3.9e-3   0.22   -2.6  -2.27  -2.11  -1.97  -1.75   3298    1.0
log_lik[32]  -1.68  2.6e-3   0.16  -2.02  -1.78  -1.67  -1.57  -1.41   3679    1.0
log_lik[33]  -1.42  2.1e-3   0.13  -1.69   -1.5  -1.41  -1.33   -1.2   3750    1.0
log_lik[34]  -1.07  2.2e-3    0.1  -1.27  -1.14  -1.07   -1.0  -0.87   2215    1.0
log_lik[35]  -1.03  2.4e-3    0.1  -1.23  -1.09  -1.02  -0.96  -0.83   1872    1.0
log_lik[36]  -2.86  6.1e-3   0.33  -3.58  -3.07  -2.84  -2.63  -2.28   2929    1.0
log_lik[37]  -1.36  2.0e-3   0.12  -1.62  -1.44  -1.36  -1.28  -1.15   3597    1.0
log_lik[38]  -1.03  2.4e-3    0.1  -1.24   -1.1  -1.03  -0.96  -0.84   1873    1.0
log_lik[39]  -1.33  2.1e-3   0.12  -1.59  -1.41  -1.33  -1.25  -1.11   3261    1.0
log_lik[40]  -1.06  2.3e-3    0.1  -1.27  -1.13  -1.06   -1.0  -0.86   2019    1.0
log_lik[41]  -1.14  2.2e-3   0.11  -1.37  -1.21  -1.14  -1.07  -0.95   2408    1.0
log_lik[42]  -1.11  2.2e-3   0.11  -1.34  -1.18  -1.11  -1.04  -0.91   2274    1.0
log_lik[43]  -1.03  2.4e-3    0.1  -1.24  -1.09  -1.02  -0.96  -0.84   1888    1.0
log_lik[44]  -1.38  2.9e-3   0.14  -1.69  -1.46  -1.36  -1.28  -1.12   2386    1.0
log_lik[45]  -1.07  2.5e-3   0.11  -1.28  -1.13  -1.06   -1.0  -0.86   1894    1.0
log_lik[46]  -1.37  2.9e-3   0.15  -1.69  -1.46  -1.36  -1.27  -1.11   2479    1.0
log_lik[47]  -1.05  2.3e-3   0.11  -1.27  -1.12  -1.05  -0.98  -0.86   2047    1.0
log_lik[48]  -1.23  2.9e-3   0.13  -1.52  -1.31  -1.22  -1.14   -1.0   2066    1.0
log_lik[49]  -1.24  3.1e-3   0.14  -1.55  -1.33  -1.23  -1.15  -1.01   1991    1.0
log_lik[50]  -1.03  2.4e-3    0.1  -1.24   -1.1  -1.03  -0.96  -0.83   1903    1.0
log_lik[51]  -2.23  6.5e-3    0.3  -2.87  -2.41  -2.21  -2.01  -1.69   2184    1.0
log_lik[52]  -1.43  4.0e-3   0.18  -1.83  -1.54  -1.42  -1.31  -1.13   2046    1.0
log_lik[53]  -1.08  2.8e-3   0.12  -1.32  -1.15  -1.08   -1.0  -0.86   1741    1.0
log_lik[54]  -1.48  4.5e-3    0.2  -1.93   -1.6  -1.47  -1.34  -1.14   1957    1.0
log_lik[55]  -1.19  3.5e-3   0.15   -1.5  -1.27  -1.18  -1.09  -0.94   1759    1.0
log_lik[56]  -1.13  3.3e-3   0.13  -1.42  -1.21  -1.12  -1.04  -0.89   1701    1.0
log_lik[57]  -1.54  6.0e-3   0.23  -2.07  -1.67  -1.51  -1.38  -1.15   1494    1.0
log_lik[58]  -1.04  2.5e-3   0.11  -1.27  -1.11  -1.04  -0.97  -0.84   1959    1.0
log_lik[59]  -1.58  6.7e-3   0.25  -2.16  -1.73  -1.55  -1.41  -1.17   1442    1.0
log_lik[60]  -1.38  5.1e-3   0.21  -1.85  -1.51  -1.37  -1.24  -1.04   1689    1.0
log_lik[61]  -1.84  8.5e-3   0.32  -2.56  -2.03   -1.8  -1.62  -1.31   1395    1.0
log_lik[62]  -1.61  6.6e-3   0.27  -2.23  -1.77  -1.58  -1.42  -1.16   1678    1.0
lp__        -56.57    0.04   1.52 -60.52 -57.31  -56.2 -55.45 -54.71   1281    1.0

Samples were drawn using NUTS at Tue Oct 23 13:03:15 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).

``````

Without standardization we'are again get smaller `n_eff`'s, but large enough for practical purposes.

Check the treedepth, E-BFMI, and divergences

``````

In [36]:

stan_utility.check_treedepth(fit)
stan_utility.check_energy(fit)
stan_utility.check_div(fit)

``````
``````

145 of 4000 iterations saturated the maximum tree depth of 10 (3.625%)
Run again with max_depth set to a larger value to avoid saturation
0.0 of 4000 iterations ended with a divergence (0.0%)

``````

We see again many iterations saturating the maximum tree depth, which is harmful for the efficiency but doesn't invalidate the results.

Compute the probability that the summer temperature is increasing.

``````

In [37]:

samples = fit.extract(permuted=True)
print('Pr(beta > 0) = {}'.format(np.mean(samples['beta'] > 0)))

``````
``````

Pr(beta > 0) = 0.99225

``````

We get similar probability as with Gaussian obervation model.

Plot the data, the model fit, and the marginal posteriors for sigma and nu.

``````

In [38]:

# make slightly wider figure of 4 plots
figsize = plt.rcParams['figure.figsize'].copy()
figsize[0] *= 1.4  # width
figsize[1] *= 3  # height
fig, axes = plt.subplots(4, 1, figsize=figsize)

# plot 1: scatterplot and lines
color_scatter = 'C0'  # 'C0' for default color #0
color_line = 'C1'     # 'C1' for default color #1
# lighten color_line
1 - 0.1*(1 - np.array(mpl.colors.to_rgb(color_line)))
)
# plot
ax = axes[0]
ax.fill_between(
x,
np.percentile(samples['mu'], 5, axis=0),
np.percentile(samples['mu'], 95, axis=0),
)
ax.plot(
x,
np.percentile(samples['mu'], 50, axis=0),
color=color_line,
linewidth=1
)
ax.scatter(x, y, 5, color=color_scatter)
ax.set_xlabel('year')
ax.set_ylabel('summer temperature at Kilpisjarvi')
ax.set_xlim((1952, 2013))

# plot 2: histogram
ax = axes[1]
ax.hist(samples['beta'], 50)
ax.set_xlabel('beta')

# plot 3: histogram
ax = axes[2]
ax.hist(samples['sigma'], 50)
ax.set_xlabel('sigma')

# plot 4: histogram
ax = axes[3]
ax.hist(samples['nu'], 50)
ax.set_xlabel('nu')

# make figure compact
fig.tight_layout();

``````
``````

``````

The posterior of nu reveals that Gaussian model is likely to be ok, as Student's t with `nu > 20` is very close to Gaussian.

# Comparison of k groups with hierarchical models

Let's compare the temperatures in three summer months.

``````

In [39]:

data_path = os.path.abspath(
os.path.join(
os.path.pardir,
'utilities_and_data',
'kilpisjarvi-summer-temp.csv'
)
)
d = np.loadtxt(data_path, dtype=np.double, delimiter=';', skiprows=1)
# Is there difference between different summer months?
x = np.tile(np.arange(1, 4), d.shape[0]) # summer months are numbered from 1 to 3
y = d[:, 1:4].ravel()
N = len(x)
data = dict(
N = N,
K = 3,  # 3 groups
x = x,  # group indicators
y = y   # observations
)

``````

## Common variance (ANOVA) model

``````

In [40]:

with open('grp_aov.stan') as file:

``````
``````

// Comparison of k groups with common variance (ANOVA)
data {
int<lower=0> N; // number of data points
int<lower=0> K; // number of groups
int<lower=1,upper=K> x[N]; // group indicator
vector[N] y; //
}
parameters {
vector[K] mu;        // group means
real<lower=0> sigma; // common std
}
model {
y ~ normal(mu[x], sigma);
}

``````

Fit the model

``````

In [41]:

model = stan_utility.compile_model('grp_aov.stan')
fit = model.sampling(data=data, seed=194838)
samples = fit.extract(permuted=True)

``````
``````

INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_7564a28d1f0313ff396284979a4361d6 NOW.

``````

Check the `n_eff` and `Rhat`

``````

In [42]:

print(fit)

``````
``````

Inference for Stan model: anon_model_7564a28d1f0313ff396284979a4361d6.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.

mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
mu[1]   7.54  3.0e-3   0.19   7.16   7.41   7.54   7.66   7.91   3977    1.0
mu[2]  10.96  2.9e-3   0.19  10.59  10.83  10.96  11.09  11.35   4410    1.0
mu[3]   9.44  3.0e-3    0.2   9.05    9.3   9.44   9.57   9.82   4309    1.0
sigma   1.53  1.2e-3   0.08   1.38   1.47   1.53   1.58    1.7   4468    1.0
lp__  -170.8    0.03   1.41 -174.5 -171.5 -170.5 -169.8 -169.1   2066    1.0

Samples were drawn using NUTS at Tue Oct 23 13:04:16 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).

``````

Check the treedepth, E-BFMI, and divergences

``````

In [43]:

stan_utility.check_treedepth(fit)
stan_utility.check_energy(fit)
stan_utility.check_div(fit)

``````
``````

0 of 4000 iterations saturated the maximum tree depth of 10 (0.0%)
0.0 of 4000 iterations ended with a divergence (0.0%)

``````

Plot group mean distributions and matrix of probabilities that one mu is larger than other.

``````

In [44]:

# matrix of probabilities that one mu is larger than other
mu = fit.extract(permuted=True)['mu']
ps = np.zeros((3, 3))
for k1 in range(3):
for k2 in range(k1+1, 3):
ps[k1, k2] = np.mean(mu[:, k1] > mu[:, k2])
ps[k2, k1] = 1 - ps[k1, k2]
print("Matrix of probabilities that one mu is larger than other:")
print(ps)

# plot violins
plt.violinplot(mu, showmeans=True, showextrema=False)
plt.xticks(np.arange(1, 4))
plt.xlabel('group')
plt.ylabel('mu');

``````
``````

Matrix of probabilities that one mu is larger than other:
[[0. 0. 0.]
[1. 0. 1.]
[1. 0. 0.]]

``````

## Common variance and hierarchical prior for mean.

Results do not differ much from the previous, because there is only few groups and quite much data per group, but this works as an example of a hierarchical model.

``````

In [45]:

with open('grp_prior_mean.stan') as file:

``````
``````

// Comparison of k groups with common variance and
// hierarchical prior for the mean
data {
int<lower=0> N; // number of data points
int<lower=0> K; // number of groups
int<lower=1,upper=K> x[N]; // group indicator
vector[N] y; //
}
parameters {
real mu0;             // prior mean
real<lower=0> sigma0; // prior std
vector[K] mu;         // group means
real<lower=0> sigma;  // common std
}
model {
mu0 ~ normal(10,10);      // weakly informative prior
sigma0 ~ cauchy(0,4);     // weakly informative prior
mu ~ normal(mu0, sigma0); // population prior with unknown parameters
sigma ~ cauchy(0,4);      // weakly informative prior
y ~ normal(mu[x], sigma);
}

``````

Fit the model

``````

In [46]:

model = stan_utility.compile_model('grp_prior_mean.stan')
fit = model.sampling(data=data, seed=194838)
samples = fit.extract(permuted=True)

``````
``````

INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_323d20ff1d7ed8fcfbd94da1ca27c19c NOW.

``````

Plot group mean distributions and matrix of probabilities that one mu is larger than other.

``````

In [47]:

# matrix of probabilities that one mu is larger than other
mu = fit.extract(permuted=True)['mu']
ps = np.zeros((3, 3))
for k1 in range(3):
for k2 in range(k1+1, 3):
ps[k1, k2] = np.mean(mu[:, k1] > mu[:, k2])
ps[k2, k1] = 1 - ps[k1, k2]
print("Matrix of probabilities that one mu is larger than other:")
print(ps)

# plot violins
plt.violinplot(mu, showmeans=True, showextrema=False)
plt.xticks(np.arange(1, 4))
plt.xlabel('group')
plt.ylabel('mu');

``````
``````

Matrix of probabilities that one mu is larger than other:
[[0. 0. 0.]
[1. 0. 1.]
[1. 0. 0.]]

``````

## Unequal variance and hierarchical prior for mean and variance

Results do not differ much from the previous, because there is only few groups and quite much data per group, but this works as an example of a hierarchical model.

``````

In [48]:

with open('grp_prior_mean_var.stan') as file:

``````
``````

// Comparison of k groups with unequal variance and
// hierarchical priors for the mean and the variance
data {
int<lower=0> N; // number of data points
int<lower=0> K; // number of groups
int<lower=1,upper=K> x[N]; // group indicator
vector[N] y; //
}
parameters {
real mu0;                 // prior mean
real<lower=0> musigma0;   // prior std
vector[K] mu;             // group means
real lsigma0;             // prior mean
real<lower=0> lsigma0s;   // prior std
vector<lower=0>[K] sigma; // group stds
}
model {
mu0 ~ normal(10, 10);       // weakly informative prior
musigma0 ~ cauchy(0,10);    // weakly informative prior
mu ~ normal(mu0, musigma0); // population prior with unknown parameters
lsigma0 ~ normal(0,1);      // weakly informative prior
lsigma0s ~ normal(0,1);     // weakly informative prior
sigma ~ cauchy(lsigma0, lsigma0s); // population prior with unknown parameters
y ~ normal(mu[x], sigma[x]);
}

``````

Fit the model

``````

In [49]:

model = stan_utility.compile_model('grp_prior_mean_var.stan')
fit = model.sampling(data=data, seed=194838)
samples = fit.extract(permuted=True)

``````
``````

INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_2f157eda9830d6c0deb5559f82c524bb NOW.

``````

Plot group mean distributions and matrix of probabilities that one mu is larger than other:

``````

In [50]:

samples = fit.extract(permuted=True)
print("std(mu0): {:.2f}".format(np.std(samples['mu0'])))
mu = samples['mu']

# matrix of probabilities that one mu is larger than other
ps = np.zeros((3, 3))
for k1 in range(3):
for k2 in range(k1+1, 3):
ps[k1, k2] = np.mean(mu[:, k1] > mu[:, k2])
ps[k2, k1] = 1 - ps[k1, k2]
print("Matrix of probabilities that one mu is larger than other:")
print(ps)

# plot violins
plt.violinplot(mu, showmeans=True, showextrema=False)
plt.xticks(np.arange(1, 4))
plt.xlabel('group')
plt.ylabel('mu');

``````
``````

std(mu0): 2.38
Matrix of probabilities that one mu is larger than other:
[[0. 0. 0.]
[1. 0. 1.]
[1. 0. 0.]]

``````