Theano素人がPyMC3で頑張るための小技

PythonのMCMC(Markov Chain Monte Carlo)ライブラリであるPyMC3は，試された方はご存知の通り，数値処理ライブラリTheanoをベースとして作られている．この組み合わせのおかげで，複雑な問題も扱える柔軟性と十分な計算処理能力持つと言われているが，いかんせんTheanoは素人には相当難しい．シンボルの扱いを先に定義し，後から数値をバインドして計算するコンセプトが理解を難しくさせている．ここでは，PyMC3 simulation作業で問題が生じたときにTheanoの機能に触らずににデバッグを行うための小技を紹介する．

---

PyMC3 simulation作業の流れ

一般的なsimulationの作業の流れは次のようなものである．

1. データセットの用意（ファイル読み込み）．
2. シミュレーションモデルの記述．（事前分布はこれとこれ，尤度関数はこれ，など．）
3. 初期値設定．MCMCプロセスによる数値計算．
4. 計算結果の考察．グラフ化など．

注意が必要なのが「統計モデルの記述」の部分である．単純なモデルであれば，Tutorialを参考に簡単にモデルを作ることができると思うが，モデルが複雑になると計算しない，あるいは，計算した結果が変だ，という事態に陥ってしまう．

デバッグのために変数の値を出力したい

変数の中間値を出力することは，デバッグの大きな助けになることは間違いないが，これが Theano & PyMC3 環境下では簡単にいかない．

例として以下のcodeを取り上げる．ロジスティック回帰のsimulationである．

import numpy as np
import pymc3 as pm

# １．データセットの用意
n_betl = np.array([59, 60, 62, 56, 63, 59, 62, 60])
y_betl = np.array([6, 13, 18, 28, 52, 53, 61, 60])
x1 = np.array([1.6907, 1.7242, 1.7552, 1.7842,
    1.8113, 1.8369, 1.8610, 1.8839])

# ２．統計モデルの記述
mymodel1 = pm.Model()
def invlogit(x):
    return pm.exp(x) / (1 + pm.exp(x))

with mymodel1:  # model definition
    mytheta = pm.Normal('theta', mu=0, sd=32, shape=2)   # 2-dim array
    p = invlogit(mytheta[0] + mytheta[1] * x1)
    y_obs = pm.Binomial('y_obs', n=n_betl, p=p, observed=y_betl)

# ３．MCMC計算
with mymodel1:  # running simulation
    start = pm.find_MAP()
    step = pm.NUTS()
    trace1 = pm.sample(10000, step, start=start)
    
# ４．計算結果出力 ... summaryとtraceplot
pm.summary(trace1)
pm.traceplot(trace1)

やりたいのは，MCMC計算の各ステップにおける変数値出力であることが多いが，どこにそのステートメントを入れるか？直観的にMCMC計算の部分に入れたいが，計算命令は，

    trace1 = pm.sample(10000, step, start=start)

の一行で行われるので，挿入することは難しい．ならばその前の「統計モデルの記述」のところでなんとかできないかと思いTheanoのドキュメントをさがし，以下のような行を追加してみた．
（**"theano.printing.Print()"**の箇所になります．）

import theano
（中略）
with mymodel1:  # model definition
    mytheta = pm.Normal('theta', mu=0, sd=32, shape=2)   # 2-dim array

	# 私は，MCMC計算途中のlogit_pを出力させたい．．．
    logit_p = mytheta[0] + mytheta[1] * x1
    my_printed = theano.printing.Print('logit_p = ')(logit_p)

    p = invlogit(mytheta[0] + mytheta[1] * x1)
    y_obs = pm.Binomial('y_obs', n=n_betl, p=p, observed=y_betl)

結果としては，初期値のオールゼロが出力された．

logit_p =  __str__ = [ 0.  0.  0.  0.  0.  0.  0.  0.]

ということでデバッグとしては意味をなさないものである．やはり統計モデル記述のところで「出力」を行ってもだめなようだ．

PyMC(MCMC計算)のtraceを使う

別の方法としてPyMCのtraceに，任意の変数を追加してロギングする方法に思いあたった．統計モデルの主要な変数は，trace dataとして各stepで記録されている．これに自分の見たい変数を追加できないかやってみた．

with mymodel1:  # model definition
    mytheta = pm.Normal('theta', mu=0, sd=32, shape=2)   # 2-dim array
    p = invlogit(mytheta[0] + mytheta[1] * x1)

    # for debug
	logit_p = pm.Deterministic('logit_p', (mytheta[0] + mytheta[1] * x1))
 
    y_obs = pm.Binomial('y_obs', n=n_betl, p=p, observed=y_betl)
   

上記の通り，pm.Deterministic() のクラスメソッドを用いた．本来，pm.Deterministicは関連する他の変数から（確定的に求められる）変数を導出するために用いるが，これを使うことにより，引数として渡した**'ligit_p'** の名前を持った変数がtrace記録の対象になる．よって計算後に，traceからデバッグ情報としてこれを参照することができるようになる．もちろんmymodel1 を構成する他の変数を参照することもできるし，適当に追加の計算式を入れて，その結果を出すことも可能である．

>>> mylog = trace1['logit_p']
>>> mylog[-5:]  ... 最後の５ステップのみを表示
array([[-2.60380812, -1.54381452, -0.56292492,  0.35468148,  1.21216885,
         2.02219381,  2.78475637,  3.50934901],
       [-2.44719254, -1.38360671, -0.39939295,  0.52132315,  1.38171646,
         2.19448653,  2.95963336,  3.68668158],
       [-2.73650878, -1.64559829, -0.63609903,  0.30827125,  1.19076899,
         2.02441999,  2.80922426,  3.55495113],
       [-2.71951223, -1.62859273, -0.61908514,  0.32529294,  1.20779796,
         2.04145585,  2.82626659,  3.57199962],
       [-2.51597252, -1.42300483, -0.41160189,  0.53454925,  1.41871117,
         2.25393425,  3.04021847,  3.78735161]])

また，traceなのでtraceplot()でplotすることもできる．

Simulationで求めたかった回帰パラメータthetaの他に，logit_pも出力されている．（logit_pが，size=8のベクトルになっていたのは意外でしたが...）

PyMC3でやっていけそうか

traceに出す変数を増やすことで，より多くのデバッグ情報を得ることができた．MCMCのモデリング作業では，各ライブラリは，あまり参考になるエラーメッセージを出してくれないので，結構「気疲れする」（さらには「気が狂いそうになる」）状況に陥る．今回の方法で少し状況が改善され，頑張っていけそうな気がしてきた．

今回は，結局Theanoのところまで「降りる」ことなく対処する方法となったが，PyMC3に限らず, Theanoは機械学習のための有効なツールであるとのことなので，こちらの方も難解だからと諦めることなく少しずつ調べていきたい．またPythonの一般的なデバッグ手法についてもtipsを増やしていきたい．

参考文献 (web site)

• PyMC3 Tutorial https://pymc-devs.github.io/pymc3/getting_started/
• Theano Tutorial http://deeplearning.net/software/theano/tutorial/
• Qiita, Theano の 基本メモ http://qiita.com/mokemokechicken/items/3fbf6af714c1f66f99e9
• Qiita, BUGS Example "Beetle" を解く http://qiita.com/TomokIshii/items/b38dfd13aec5300a86dd