PyStanで単回帰

StanとRでベイズ統計モデリング (Wonderful R)

この本のp34〜の単回帰をPyStanで実装します. まずはデータセットをdata-salary.txtという名前で保存. 次に、以下を同じディレクトリに置き、pythonコードを実行します.

• data-salary.txt
• simlinreg.stan
• simlinreg.py

//simlinreg.stan
data {
    int N;
    real X[N];
    real Y[N];
}

parameters {
    real a;
    real b;
    real<lower=0> sigma;
}

model {
    for (n in 1:N) {
        Y[n] ~ normal(a + b*X[n], sigma);
    }
}

# simlinreg.py
import pystan
import pandas as pd
import matplotlib.pyplot as plt

# load data
d = pd.read_csv("data-salary.txt")

# set data
data = {
    'N': len(d.index),
    'X': d['X'],
    'Y': d['Y']
}

# exe
fit = pystan.stan(file='simlinreg.stan', data=data, iter=1000, chains=4)

# show result
print(fit)

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

mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
a     -115.5    3.24  73.82 -252.1 -165.0 -120.0 -65.74  36.57  520.0   1.01
b      21.82    0.07   1.65  18.35   20.7   21.9  22.92  24.86  519.0   1.01
sigma  86.18    0.65  15.63  62.71  75.16  83.98  94.59 122.52  576.0   1.01
lp__  -93.63    0.06    1.3 -97.03 -94.26 -93.26 -92.66 -92.15  492.0   1.01
    
Samples were drawn using NUTS at Sun Apr 23 17:32:48 2017.
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).
"""

# show result
fit.plot()
plt.show()

Rのlm関数による推定値がa=-119.7、b=21.9なのでいい感じに推定できてますね.