犬4匹本の輪読会などでどっぷりStanにはまっていて、Pythonから使えるPyStanに入門してみました。StanはRから使う人が多くて情報もRの方が多いんだけど、慣れているPythonから使えた方が楽なのでということで。
このシリーズではPyStanならではの使い方をなるべく書いていきます。Stanの使い方はアヒル本を見ていただくのが圧倒的にわかりやすいでしょうし、私も勉強中です。
Stanとは
Pystanのインストール
Pystanのインストールには、以下の手順を踏みますが、大概の人は1行目は終わっていると思うので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)
.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()
としてやることで、以下のような可視化ができます。
可視化
自分でサンプリングデータをいじりたい場合は、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)
傾きであるbeta=2は狭い範囲で予想できていることがわかります。alphaとsigmaは範囲が広い結果となりました。
まとめ
とりあえず動かすことができました。
次回は、可視化について書く予定。