StanのPythonバインディングであるPyStanが公開されて久しいですが、検索してもあんまり情報がヒットしません。ちょっと寂しいと思ったので、インストールやtraceplotの出力なども含めて、以下の本の5.1節「重回帰」の一部を実行してみました(ステマです)。
StanとRでベイズ統計モデリング (Wonderful R)
- 作者:健太郎, 松浦
- 発売日: 2016/10/25
- メディア: 単行本
本自体の紹介は以前の記事を読んでいただければと思います。
インストール
Windows 7 64bit、Python 3系でのインストール手順を説明します。
- AnacondaでPythonなどをインストール
- コマンドプロンプトから
pip install pystan
でpystan
をインストール - Visual C++ Build Toolsをインストール
以上でPyStanが動くようになります。その他のOSにつきましては公式ページを読んでください。ほぼ同様の手順であり、難しくありません。
なおこの記事で使用したバージョンは次の通りです。
- Python 3.5.2 | Anaconda 4.2.0 (64-bit)
- PyStan 2.12.0.0
モデルファイル(Stanファイル)の説明
本の5.1節では、「大学の授業の出席率」を応答変数とし、過去の「アルバイトが好きかどうかというアンケート結果」と過去の「テストの点数」を説明変数として重回帰を行います。サポートページであるGitHubリポジトリに各種ファイルを置いてあります。
データやStanファイルの説明は本に任せて、この記事ではPythonからの使い方に注力します。
Stanファイルのコンパイルとサンプリング
pystan
を使ったStanファイルのコンパイルとサンプリングは以下のように行います。
import pandas as pd from pystan import StanModel import matplotlib.pyplot as plt import pickle d = pd.read_csv('input/data-attendance-1.txt') d.Score /= 200 data = d.to_dict('list') data.update({'N':len(d)}) stanmodel = StanModel(file='model/model5-3.stan') # NUTS (No U-Turn Sampler) fit_nuts = stanmodel.sampling(data=data, n_jobs=1) mcmc_sample = fit_nuts.extract() mu_est = mcmc_sample['mu'] # ADVI (Automatic Differentiation Variational Inference) fit_vb = stanmodel.vb(data=data) vb_sample = pd.read_csv(fit_vb['args']['sample_file'].decode('utf-8'), comment='#') vb_sample = vb_sample.drop([0,1]) mu_est = vb_sample.filter(regex='mu\.\d+') with open('output/model_and_result.pkl', 'wb') as f: pickle.dump(stanmodel, f) pickle.dump(fit_nuts, f)
- 6~9行目: データを読み込んで
pandas
のDataFrame
を作り、そこからStanに渡せるようにdictionary
に変換しています。7行目ではデータのScore
をスケーリングしています。推定したいパラメータのスケールを同じぐらいに揃えるのが目的です(収束が改善する場合があります)。 - 11行目: Stanファイルのコンパイルだけを行っています。初期値を変えたり
iter
を変えたりしてサンプリングをいろいろ試したい場合に便利です。 - 14行目: NUTSというアルゴリズムでサンプリングを実行しています。デフォルトだとコア数の分だけ並列で実行されます。しかし、Windowsの場合は無限コンパイル地獄に落ちる不具合があるので、
n_jobs=1
で並列を止める必要があります。このページを参照。 - 15~16行目: 結果からMCMCサンプルを入手する一例です。
- 19~22行目: NUTSではなく、自動変分ベイズ(ADVI)で推定する場合です。結果の形式がまだ洗練されておらず、ファイルから読み込んでからゴミの行を除去する必要があります(簡単ですが)。
- 24~26行目: コンパイルされたモデルと結果を
pickle
で保存しておきます。
実行速度は基本的にC++コンパイラに依存するだけなので、Rの場合(少し古いgcc)と大きな差はなさそうです。
収束の要約とtraceplotの保存
本の4.4.8項「収束診断をファイルへ出力する」に相当します。pystan
のデフォルトのプロット関数を使ってfit_nuts.plot()
のようにするとtraceplotを描くことができますが、4つのchainを線種を変えて描くのではなくて1つのchainにまとめて描いてしまったり、表示がつぶれたりする不満点があります。代わりの関数としては、PyMCの開発チームが作成中のmcmcplotlibに期待していますが、例も使い方も載ってないので躊躇してしまいます。ここでは自作でRの{coda}
パッケージ風のtraceplotを、複数ページのpdfで作成することにしました。
import pickle import sys import numpy as np import math import matplotlib.pyplot as plt import seaborn as sns from matplotlib.backends.backend_pdf import PdfPages with open('output/model_and_result.pkl', 'rb') as f: stanmodel = pickle.load(f) fit_nuts = pickle.load(f) ## summary with open('output/fit-summary.txt', 'w') as f: f.write(str(fit_nuts)) ## traceplot palette = sns.color_palette() ms = fit_nuts.extract(permuted=False, inc_warmup=True) iter_from = fit_nuts.sim['warmup'] iter_range = np.arange(iter_from, ms.shape[0]) paraname = fit_nuts.sim['fnames_oi'] num_pages = math.ceil(len(paraname)/4) with PdfPages('output/fit-traceplot.pdf') as pdf: for pg in range(num_pages): plt.figure() for pos in range(4): pi = pg*4 + pos if pi >= len(paraname): break plt.subplot(4, 2, 2*pos+1) plt.tight_layout() [plt.plot(iter_range + 1, ms[iter_range,ci,pi], color=palette[ci]) for ci in range(ms.shape[1])] plt.title(paraname[pi]) plt.subplot(4, 2, 2*(pos+1)) plt.tight_layout() [sns.kdeplot(ms[iter_range,ci,pi], color=palette[ci]) for ci in range(ms.shape[1])] plt.title(paraname[pi]) pdf.savefig() plt.close()
- 9~11行目:
pickle
で保存した結果を読み込んでいます。 - 14~15行目: MCMCサンプルの要約(Rhatなどの収束診断を含む)をファイルに出力しています。
- 18~40行目: traceplotを描いています。
- 18行目: chainごとの色を
seaborn
パッケージのデフォルトカラーパレットで決めています。このままですとchainの数が6を超えたらエラーが出ます(色が巡回するようにすればOKと思います)。 - 19行目:
extract
関数はMCMCサンプルを取り出す関数ですが、デフォルトだと各chainを混ぜて、iterationの順番も混ぜて出力します。開発チームの想定と異なる使われ方を防ぐためだそうです(ここ参照)。traceplotをプロットするにはchainごとにiterationの順番を保持したままのMCMCサンプルが欲しいので、このように引数を指定しています。 - 20~21行目: ここではRの
{rstan}
パッケージと同様に、warmup期間を除いた部分のtraceplotと密度関数を描くようにしました。もしwarmup期間も含めてプロットしたい場合にはiter_from = 0
としてください。 - 22~23行目: 複数ページのpdfを作るためにパラメータ数を取得しています。1ページに4つのパラメータを描きます。
- 25~40行目: 残りは
matplotlib
の話なのでここでは省略します。複数ページのpdfの作成方法は公式ページを参考にしました。
- 18行目: chainごとの色を
結果の図は以下になります(pdfの1ページ目です)。
実測値と予測値のプロット
本の5.1.8項「図によるモデルのチェック」に相当します。図5.3を描きます。横軸に実測値(データ)・縦軸に予測分布の中央値と予測区間をプロットしたグラフです。
このようなグラフを描くにはMCMCサンプルを集計したり加工したりする必要があります。岩波DS1のGitHubリポジトリ内に、sinhrksさん(pandas
のコミッタ)に書いてもらったコードがありますので、よくそれらを参考にしています(example1,example2,example3)。
Pythonコードは以下になりました。
import pickle import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns with open('output/model_and_result.pkl', 'rb') as f: stanmodel = pickle.load(f) fit_nuts = pickle.load(f) d_ori = pd.read_csv('input/data-attendance-1.txt') ms = fit_nuts.extract() quantile = [10, 50, 90] colname = ['p' + str(x) for x in quantile] d_qua = pd.DataFrame(np.percentile(ms['y_pred'], q=quantile, axis=0).T, columns=colname) d = pd.concat([d_ori, d_qua], axis=1) d0 = d[d.A == 0] d1 = d[d.A == 1] palette = sns.color_palette() fig = plt.figure() ax = fig.add_subplot(1,1,1) ax.plot([0,0.5], [0,0.5], 'k--', alpha=0.7) ax.errorbar(d0.Y, d0.p50, yerr=[d0.p50 - d0.p10, d0.p90 - d0.p50], fmt='o', ecolor='gray', ms=10, mfc=palette[0], alpha=0.8, marker='o') ax.errorbar(d1.Y, d1.p50, yerr=[d1.p50 - d1.p10, d1.p90 - d1.p50], fmt='o', ecolor='gray', ms=10, mfc=palette[1], alpha=0.8, marker='^') ax.set_aspect('equal') ax.set_xlim(0, 0.5) ax.set_ylim(0, 0.5) ax.set_xlabel('Observed') ax.set_ylabel('Predicted') fig.savefig('output/fig5-3.png')
残りはmatplotlib
の話なのでここでは省略します。結果の図は以下になります。
MCMCサンプルの散布図行列
同じく本の5.1.8項「図によるモデルのチェック」の、図5.5に相当する図を描きます。Pythonで散布図行列を描くにはseaborn
パッケージのPairGrid
を使うのが便利です。ただし、上三角行列の部分はカスタマイズする必要があります。ここを参考にしました。
import pickle import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns import math from matplotlib.patches import Ellipse from scipy import stats with open('output/model_and_result.pkl', 'rb') as f: stanmodel = pickle.load(f) fit_nuts = pickle.load(f) d_ori = pd.read_csv('input/data-attendance-1.txt') ms = fit_nuts.extract() def corrfunc(x, y, **kws): r, _ = stats.spearmanr(x, y) ax = plt.gca() ax.axis('off') ellcolor = plt.cm.RdBu(0.5*(r+1)) txtcolor = 'black' if math.fabs(r) < 0.5 else 'white' ax.add_artist(Ellipse(xy=[.5, .5], width=math.sqrt(1+r), height=math.sqrt(1-r), angle=45, facecolor=ellcolor, edgecolor='none', transform=ax.transAxes)) ax.text(.5, .5, '{:.0f}'.format(r*100), color=txtcolor, fontsize=28, horizontalalignment='center', verticalalignment='center', transform=ax.transAxes) d = pd.DataFrame({'b1':ms['b1'], 'b2':ms['b2'], 'b3':ms['b3'], 'sigma':ms['sigma'], 'mu1':ms['mu'][:,0], 'mu50':ms['mu'][:,49], 'lp__':ms['lp__']}, columns=['b1', 'b2', 'b3', 'sigma', 'mu1', 'mu50', 'lp__']) sns.set(font_scale=2) g = sns.PairGrid(d) g = g.map_lower(sns.kdeplot, cmap='Blues_d') g = g.map_diag(sns.distplot, kde=False) g = g.map_upper(corrfunc) g.fig.subplots_adjust(wspace=0.05, hspace=0.05) for ax in g.axes.flatten(): for t in ax.get_xticklabels(): _ = t.set(rotation=40) g.savefig('output/fig5-5.png')
- 17~26行目: 上三角行列を作成しています。Spearmanの順位相関係数に応じて、楕円の太さ細さと色を変えています。数値は相関係数×100を表示しています。
- 36行目: subplotの間の空白を狭くしています。空いていると少し見づらかったので。
- 37~39行目: X軸ラベルを少し傾けています。
PairGrid
なので少しややこしいです。
結果の図は以下になります。