PyStan入門してみた(1) - カレーちゃんブログ

犬４匹本の輪読会などでどっぷりStanにはまっていて、Pythonから使えるPyStanに入門してみました。StanはRから使う人が多くて情報もRの方が多いんだけど、慣れているPythonから使えた方が楽なのでということで。
このシリーズではPyStanならではの使い方をなるべく書いていきます。Stanの使い方はアヒル本を見ていただくのが圧倒的にわかりやすいでしょうし、私も勉強中です。

Stanとは
Pystanのインストール
簡単な線形回帰を動かしてみる
まとめ

Stanとは

高速なMCMCサンプラー
モデルを書いて、データを入れればサンプリングしてくれる
R、Python、Julia, などから使える
詳しくは、公式や公式（日本語）を参照

Pystanのインストール

Pystanのインストールには、以下の手順を踏みますが、大概の人は１行目は終わっていると思うのでpip install pystanで終わりでしょうか。

Python3.3以上、Cython0.22以上、NumPy1.7以上をインストール
Pystanをインストール

簡単な線形回帰を動かしてみる

データ作成

公式のGetting Startedとかそのあとのページが良くできているので、それを見つつ、とても簡単なサンプルを作って、動かしてみます。
サンプルは、Yが2×X+誤差（誤差は平均0、分散5の正規分布）となるようなデータにしてみます。

import numpy as np
import pandas as pd
import pystan
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline

X = np.array([i * 2 for i in range(1,21)])
np.random.seed(1000)
epsilon = np.random.normal(0,5,20)
Y = 2 * X + epsilon
plt.scatter(X, Y)

f:id:currypurin:20171001071539p:plain

.stanファイル

.stanファイルでモデルを記述して（今回はliner.stanというファイル名で.stanファイルを作成しています。）

data { 
  int<lower=0> N; // 学習データの数
  real X[N]; // 単回帰モデルの入力x
  real Y[N]; // 単回帰モデルの出力y
}
parameters { // parameterセクション
  real alpha; // 単回帰モデルの切片
  real beta; // 単回帰モデルのxの係数
  real<lower=0> sigma; // 単回帰モデルからの誤差の標準偏差
}
model { // モデルを宣言するmodelセクション
  for (n in 1:N)  {
    Y[n]  ~  normal(alpha + beta * X[n], sigma); // 線形単回帰モデル
  }
}

pystan.stan

PyStanの場合は、辞書型のデータをpystan.stanに渡します。

stan_data = {'N': 20, 'X':X, 'Y': Y}
fit = pystan.stan(file='./liner.stan',data=stan_data, iter=3000, chains=3, thin=1)

.stanのコンパイルとサンプリングを別にしたい場合には、以下のように書くとよいです。
こうすると、時間がかかるコンパイルが一度で済むので、基本はこの方法を使うとよいみたい。

sm = pystan.StanModel(file='liner.stan')
fit = sm.sampling(data=stan_data, iter=3000, chains=3,thin=1)

> fit
Inference for Stan model: anon_model_fc2a1126cc783f976d050e6d29c933e3.
3 chains, each with iter=3000; warmup=1500; thin=1; 
post-warmup draws per chain=1500, total post-warmup draws=4500.

        mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
alpha  -0.51    0.04   1.72  -3.93  -1.65  -0.48   0.63   2.83   1906    1.0
beta    2.01  1.7e-3   0.07   1.87   1.96   2.01   2.06   2.15   1850    1.0
sigma   3.59    0.02   0.68   2.55   3.11   3.49   3.94   5.23   1752    1.0
lp__  -33.49    0.04   1.34 -37.05 -34.08 -33.12 -32.52 -31.97   1366    1.0

Samples were drawn using NUTS at Sun Oct  1 07:12:22 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).

のようにMCMCが終わったことがわかります。
fit.plot()としてやることで、以下のような可視化ができます。 f:id:currypurin:20171001074406p:plain

可視化

自分でサンプリングデータをいじりたい場合は、samples = fit.extract(permuted=True)とすると、全てのチェインを統合したものがsamplesに入るので、後はpandasを使って可視化する。この辺りは後日詳しくまとめたいですが、とりあえず以下のような感じ。

samples = fit.extract(permuted=True)
df = pd.DataFrame(samples)

plt.plot(X,Y,"o")
plt.plot(X,samples["alpha"].mean() + X_base * samples["beta"].mean())
plt.xlim(0,42)
plt.xlabel("X")
plt.ylabel("Y")
plt.title("Pystanでの線形回帰")

sns.violinplot(data=df.iloc[:,0:3])
plt.title("alpha,beta,sigmaのバイオリンプロット")
plt.ylim(-2,8)

f:id:currypurin:20171001075931p:plain f:id:currypurin:20171001075945p:plain