読者です 読者をやめる 読者になる 読者になる

kikuの四苦Hack

ガジェットやプログラミングで遊ぶ

MacでPyStan環境構築

こういう本を買ってみたのでpythonでやってみます.

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

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

# pipでPyStanをインストール
pip install pystan

# cythonとnumpyがインストールされる (ない場合).

本家を参考にテストコードを動かします.
作業ディレクトリに移動して以下のファイルを作成->pythonコードを実行.

//8schools.stan
data {
  int<lower=0> J; // number of schools 
  real y[J]; // estimated treatment effects
  real<lower=0> sigma[J]; // s.e. of effect estimates 
}
parameters {
  real mu; 
  real<lower=0> tau;
  real eta[J];
}
transformed parameters {
  real theta[J];
  for (j in 1:J)
    theta[j] <- mu + tau * eta[j];
}
model {
  eta ~ normal(0, 1);
  y ~ normal(theta, sigma);
}
# pystan_test.py
import pystan
import matplotlib.pyplot as plt

schools_dat = {'J': 8,
               'y': [28,  8, -3,  7, -1,  1, 18, 12],
               'sigma': [15, 10, 16, 11,  9, 11, 10, 18]}

fit = pystan.stan(file='8schools.stan', data=schools_dat, iter=1000, chains=4)

# 結果出力
print(fit)

"""
Inference for Stan model: anon_model_286b3180dfa752c4cfedaf0241add0e4.
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
mu          7.9    0.18    5.1  -2.06   4.67   7.75  11.09  18.55  787.0   1.01
tau        6.35    0.18   5.33   0.18   2.26   5.03   9.11  20.19  915.0    1.0
eta[0]     0.38    0.02   0.91  -1.53  -0.18   0.39   0.98   2.07 1734.0    1.0
eta[1]     0.04    0.02    0.9  -1.77  -0.52   0.04    0.6   1.93 1550.0    1.0
eta[2]    -0.22    0.02   0.95  -2.03  -0.86  -0.23   0.42   1.69 1962.0    1.0
eta[3]    -0.04    0.02   0.89  -1.76  -0.64  -0.06   0.57   1.73 1809.0    1.0
eta[4]    -0.35    0.02   0.88  -2.03  -0.98  -0.37    0.2   1.43 1755.0    1.0
eta[5]    -0.21    0.02   0.92  -1.99  -0.84  -0.18   0.39   1.55 1431.0    1.0
eta[6]     0.33    0.02    0.9  -1.42  -0.25   0.32   0.92   2.09 1467.0    1.0
eta[7]     0.03    0.02   0.99  -1.92  -0.63   0.04   0.69   1.96 1686.0    1.0
theta[0]  11.16    0.23   8.13  -2.48   5.69  10.19  15.32  31.45 1263.0    1.0
theta[1]   8.05    0.14   6.41  -5.07   4.11   7.84  11.96  20.81 2000.0    1.0
theta[2]   6.25    0.22   7.97  -11.8   1.73   6.63  11.09  21.32 1263.0    1.0
theta[3]   7.64    0.15   6.73  -5.77    3.4   7.64  11.51   21.6 2000.0    1.0
theta[4]   5.21    0.16   6.56  -9.48   1.43   5.68   9.55  16.98 1719.0    1.0
theta[5]   6.05    0.15   6.73  -9.31   2.11   6.41  10.36  18.54 1903.0    1.0
theta[6]  10.47    0.15   6.84  -1.75   5.91   9.74  14.38  26.55 2000.0    1.0
theta[7]   8.58    0.23   8.12  -6.02   3.61   8.05   12.9  26.96 1295.0   1.01
lp__      -5.12     0.1   2.74 -11.25  -6.89  -4.82  -3.16  -0.43  702.0    1.0

Samples were drawn using NUTS at Sun Apr 23 16:19:23 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).
"""

# 結果描画
fit.plot()
plt.show()

f:id:kiku11:20170423162505p:plain
結果の見方がようわからん...^^;