『StanとRでベイズ統計モデリング』をpythonでやってみる

TL;DR

今後は以下の書籍を読んで勉強した内容をブログに書いていこうと思います(｀・ω・´)

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

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

選定理由

• 現在興味を持っている ベイズ統計モデリング の書籍では, 最も網羅的かつ実践的な内容に思えるため（他の書籍だとモデルやサンプリングの結果の評価が薄い点が気になる. 実際の業務で使う際には評価の観点がとても重要になると考えるため）
• pythonの本だとコードをコピペして理解した気になってしまうため, Stan/RStanという慣れていない言語の方が都合がよい
• 可視化や図が多くてわかりやすそう
• 表紙がかわいい(｀・ω・´)

やり方

• 30分で読む・ブログに書くが完結できる分量にする. もしそれ以上かかるようなら読む分量を分割する
• RStanのコードはできるだけpythonで置き換えて実行する. コードをコピペするだけだと理解できてないことが多いが, 別言語に置き換えるにはちゃんと理解してないとできない.
• 当たり前だが, 自分の考えと書籍や参考資料に書いてある内容を明確に分ける

参考文献

• StanとRでベイズ統計モデリング (Wonderful R)
• 「StanとRでベイズ統計モデリング」松浦健太郎 という本を書きました - StatModeling Memorandum

pymc4の現状整理

ちょっとブログ更新の間があいてたのですが、そろそろ再開したいと思います.

pymc4のソースコードリーディングがだいたい完了したので、現状のライブラリの現状を整理します.

現状整理

• 主な処理はtensorflow probabilityを内部的に使用
• 可視化はarvizを使用
• 他ライブラリを使いやすくするための薄いラッパー（2018/10/30現在: 56コミットしかない）
• メインコントリビューターがarviz等の整備に注力しているため、pymc4自体の開発は進んでない（逆に言えば他ライブラリの整備の優先度が高い）

今後やっていきたいこと

• [ ] pymc4にコントリビュートしていく（ほとんどの処理がtensorflow probabilityにあるため、主にREADMEやexampleの整備になる想定）
• [ ] tensorflow probabilityのソースコードリーディング（深く理解していくために必要）
• [ ] StanとRでベイズ統計モデリング (Wonderful R) を読み進めて、pymc4で実行してみる

参考資料

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

サンプリングのテスト追加と微修正[6717b3e, 007cd06, d896372] - pymc4のソースコード読んでみた

TL;DR

• 主にサンプリングのテスト追加

コミット

2018/07/06のコミットです.

• fix pylint error · pymc-devs/pymc4@6717b3e · GitHub
• add test for sampling · pymc-devs/pymc4@007cd06 · GitHub
• solve pycodestyle errors · pymc-devs/pymc4@d896372 · GitHub

以下が変更対象ファイルです.

• pymc4/model/base.py
  • 不要な変数を削除
• tests/test_sampling.py
  • テスト追加
  • pycodestyle対応

tests/test_sampling.py

正規分布の確率分布からのサンプリングをテストし, 期待通り平均0/標準偏差1になっていることを確認している.

def test_sample():
    model = pm.Model()
    
    @model.define
    def sample(cfg):
        mu = ed.Normal(0., 1., name="mu")
    trace = pm.sample(model) # => Acceptance rate: 0.9916
    assert 0. == pytest.approx(trace["mu"].mean(), 1) # => 0.014200708
    assert 1. == pytest.approx(trace["mu"].std(), 1) # => 0.9805069

pymc4/inference/sampling/sample.py

ついでに sample 関数の実装を確認する. 主に次の3つの処理を実行している.

1. サンプリング時の初期値を0.5で設定（乱数ではない）
2. tfp.mcmc （tensorflow_probabilityのMCMC実装）でサンプリング処理を定義
3. tf.Session 内でサンプリング処理を実行
   1. サンプリング後の状態を取得
   2. acceptされたサンプル数から%を計算

def sample(model,  # pylint: disable-msg=too-many-arguments
           num_results=5000,
           num_burnin_steps=3000,
           step_size=.4,
           num_leapfrog_steps=3,
           numpy=True):
    initial_state = []
    for name, (_, shape, _) in model.unobserved.items():
        initial_state.append(.5 * tf.ones(shape, name="init_{}".format(name)))

    states, kernel_results = tfp.mcmc.sample_chain(
        num_results=num_results,
        num_burnin_steps=num_burnin_steps,
        current_state=initial_state,
        kernel=tfp.mcmc.HamiltonianMonteCarlo(
            target_log_prob_fn=model.target_log_prob_fn(),
            step_size=step_size,
            num_leapfrog_steps=num_leapfrog_steps))

    if numpy:
        with tf.Session() as sess:
            states, is_accepted_ = sess.run([states, kernel_results.is_accepted])
            accepted = np.sum(is_accepted_)
            print("Acceptance rate: {}".format(accepted / num_results))
    return dict(zip(model.unobserved.keys(), states))

各々の処理を見ていく.

---

1. サンプリング時の初期値を0.5で設定（乱数ではない）

unobserved の数分の0.5の行列のリストを生成している.

    initial_state = []
    for name, (_, shape, _) in model.unobserved.items():
        initial_state.append(.5 * tf.ones(shape, name="init_{}".format(name)))

1. tfp.mcmc （tensorflow_probabilityのMCMC実装）でサンプリング処理を定義

    states, kernel_results = tfp.mcmc.sample_chain(
        num_results=num_results, # => 結果として必要なサンプリング数
        num_burnin_steps=num_burnin_steps, # => 最初の方の削除対象の不要なサンプリング数
        current_state=initial_state, # => 初期値（すべて0.5の行列）
        kernel=tfp.mcmc.HamiltonianMonteCarlo( # => サンプリングにはハミルトンモンテカルロ法を指定
            target_log_prob_fn=model.target_log_prob_fn(),
            step_size=step_size,
            num_leapfrog_steps=num_leapfrog_steps))

1. tf.Session 内でサンプリング処理を実行
   1. サンプリング後の状態を取得
   2. acceptされたサンプル数から%を計算

この辺りの細かい挙動は tfp.mcmc 内の実装を実際に読んでいく必要がありそう.

    if numpy:
        with tf.Session() as sess:
            states, is_accepted_ = sess.run([states, kernel_results.is_accepted])
            accepted = np.sum(is_accepted_)
            print("Acceptance rate: {}".format(accepted / num_results))
    return dict(zip(model.unobserved.keys(), states))

TODO

• [ ] tfp.mcmc 辺りの実装を確認する

参考資料

• Get mcmc sampling to work by sharanry · Pull Request #9 · pymc-devs/pymc4 · GitHub

対数確率計算用インターセプタ（実装はちょっと怪しい）[f223e4e] - pymc4のソースコード読んでみた

TL;DR

• 対数確率計算用インターセプタの実装を読んだが, なんだか実装が怪しい気がする.

コミット

2018/07/01のコミットです.

• make target_log_prob_fn to use class based interceptor · pymc-devs/pymc4@f223e4e · GitHub

以下が変更対象ファイルです.

• pymc4/init.py
  • import追加: from .inference.sampling.sample import sample
• pymc4/model/base.py
• pymc4/util/interceptors.py

CollectLogProb - pymc4/util/interceptors.py

CollectLogProb クラスの親クラスは SetState なのでそちらから確認する.

• 引数の states はdict(key=rvの名前: value=rvの実体インスタンス)
  • observed
  • unobserved

class SetState(Interceptor):
    def __init__(self, state):
        self.state = state

    def before(self, f, *args, **kwargs):
        if kwargs['name'] in self.state:
            kwargs['value'] = self.state[kwargs['name']]
        return f, args, kwargs

以前の記事 edward2のinterception処理 を読み返すと, 次のような関数の引数が **kwargs に対応する. 上記の SetState では name がkeyword引数で渡される. 例) name=“y”.

ed.Poisson(rate=1.5, name="y")

擬似コードで説明してみる.

# keyword引数
kwargs = {"name": "y"}

# すべてのRV変数の辞書
state = {"x": "foo", "y": "bar"}

if kwargs['name'] in state: # kwargs['name'] == "y"
    kwargs['value'] = state[kwargs['name']]

print(kwargs) # => {"name": "y", "value": "bar"}

次に CollectLogProb の方を見てみる. 単にRV変数すべての対数確率を計算して足してるだけのはずです.

class CollectLogProb(SetState):
    def __init__(self, states):
        super().__init__(states)
        self.log_probs = []

    def before(self, f, *args, **kwargs):
        if kwargs['name'] not in self.state:
            raise RuntimeError(kwargs.get('name'), 'All RV should be present in state dict')
        return super().before(f, *args, **kwargs)

    def after(self, rv, *args, **kwargs):
        name = kwargs.get("name")
        for name in self.state:
            value = self.state[name]
            if kwargs.get("name") == name:
                kwargs["value"] = value
        log_prob = tf.reduce_sum(rv.distribution.log_prob(rv.value))
        self.log_probs.append(log_prob)
        return rv

    @property
    def result(self):
        return self.log_probs

しかし, after の処理は少し怪しい気がする.

    def after(self, rv, *args, **kwargs):
        name = kwargs.get("name")
        for name in self.state:
            value = self.state[name]
            if kwargs.get("name") == name:
                kwargs["value"] = value
        log_prob = tf.reduce_sum(rv.distribution.log_prob(rv.value))
        self.log_probs.append(log_prob)
        return rv

例えば, kwargs.get("name") の処理がおかしい気がする. name の値を辿ってみる.

name = kwargs.get("name") # 仮にname変数の値を `test` とする.
for name in self.state: # self.stateの中に `test` というkeyが存在するかチェック（ifでよいのでは？）
    value = self.state[name]
    if kwargs.get("name") == name: # name変数の値は `test` で、かつ `kwargs.get("name")` と同じ値なのでif文は不要では?
        kwargs["value"] = value

いろいろ気になるが, 次のように書き換えられる気がする. log_prob の処理もif文の中に含めないと, state に存在しないrv変数の対数確率も計算してしまう.

    def after(self, rv, *args, **kwargs):
        name = kwargs.get("name")
        if name in self.state:
            value = self.state[name]
            kwargs["value"] = value # TODO: このvalueをどこで使用してるかよくわかっていない
            log_prob = tf.reduce_sum(rv.distribution.log_prob(rv.value))
            self.log_probs.append(log_prob)
        return rv

TODO

• [ ] Modelクラスの self.observed + self.unobserved = self.varaibles のように見えるので, states = self.variables と単にしてよいのではないかという疑問
• [ ] CollectLogProb の after 処理の挙動が怪しい気がするので確認したい

参考資料

• Get mcmc sampling to work by sharanry · Pull Request #9 · pymc-devs/pymc4 · GitHub
• edward2のinterception処理 e334115, d07338e, 93bc07b - pymc4のソースコード読んでみた - オーストラリアで勉強してきたデータサイエンティストの口語自由詩

対数確率関数の計算をインターセプター処理 [f223e4e] - pymc4のソースコード読んでみた

TL;DR

• 変数に対する処理は基本的にインターセプターで対応する方針のため, 対数確率も同様に対応.

コミット

2018/07/01のコミットです.

• make target_log_prob_fn to use class based interceptor · pymc-devs/pymc4@f223e4e · GitHub

以下が変更対象ファイルです.

• pymc4/init.py
  • import追加: from .inference.sampling.sample import sample
• pymc4/model/base.py
• pymc4/util/interceptors.py

target_log_prob_fn - pymc4/model/base.py

states 変数に unobserved と observed のkey:valueを格納し, interceptors.CollectLogProb で states のインターセプト処理をしている. この CollectLogProb 内で対数確率の計算を行っている. その結果がリストで返るため, 最後に sum で合計を計算している.

    def target_log_prob_fn(self, *args, **kwargs):  # pylint: disable=unused-argument
        """
        Pass the states of the RVs as args in alphabetical order of the RVs.
        Compatible as `target_log_prob_fn` for tfp samplers.
        """

        def log_joint_fn(*args, **kwargs):  # pylint: disable=unused-argument
            states = dict(zip(self.unobserved.keys(), args))
            states.update(self.observed)
            log_probs = []
            interceptor = interceptors.CollectLogProb(states)
            with ed.interception(interceptor):
                self._f(self._cfg)

            log_prob = sum(interceptor.log_probs)
            return log_prob
        return log_joint_fn

CollectLogProb - pymc4/util/interceptors.py

• 次回読みます.

参考資料

• Get mcmc sampling to work by sharanry · Pull Request #9 · pymc-devs/pymc4 · GitHub

pycodestyle追加ともろもろ[a7cef9b, bd381b1, 89edc5c, 9f46878] - pymc4のソースコード読んでみた

TL;DR

• pycodestyle でコードチェックの追加
• 不要なテストの削除

コミット

2018/07/01のコミットです.

• remove some tests · pymc-devs/pymc4@a7cef9b · GitHub
• fix pycodestyle errors · pymc-devs/pymc4@bd381b1 · GitHub
• remove test_interceptors · pymc-devs/pymc4@89edc5c · GitHub
• add pycodestyle to dev requirements · pymc-devs/pymc4@9f46878 · GitHub

以下が各コミットの概要です.

• a7cef9b
  • クラスのインスタンス化のテストなど, やり過ぎなテストを削除
• bd381b1
  • pycodestyle でコードスタイルをチェックした対応
• 89edc5c
  • a7cef9b と同様にクラスのインスタンス化のテストを削除
• 9f46878
  • pycodestyle をrequirements-dev.txtに追加

pycodestyle - 9f46878

pycodestyleとは?

pycodestyle: GitHub - PyCQA/pycodestyle: Simple Python style checker in one Python file

一言で言うと, pep8がpycodestyleという名前に変わった だけらしいです.

see: pep8 が pycodestyle に変わった話

以下のようにコードチェックをできるようです.

例）複数importする際には, , 区切りは不可

$ pycodestyle --show-source --show-pep8 testsuite/E40.py
testsuite/E40.py:2:10: E401 multiple imports on one line
import os, sys
         ^
    Imports should usually be on separate lines.

    Okay: import os\nimport sys
    E401: import sys, os

参考資料

• Get mcmc sampling to work by sharanry · Pull Request #9 · pymc-devs/pymc4 · GitHub
• GitHub - PyCQA/pycodestyle: Simple Python style checker in one Python file
• pep8 が pycodestyle に変わった話

（細かい修正なのであまり重要ではない） [c27e97f, a8e3dae, 4f5382f, e06d946, ca9f334] - pymc4のソースコード読んでみた

TL;DR

• 細かい修正が多いのであまり重要ではないです.

コミット

2018/06/26から2018/06/30の間のコミットです.

• minor fixes · pymc-devs/pymc4@c27e97f · GitHub
• add some tests · pymc-devs/pymc4@a8e3dae · GitHub
• add tests and fix model.configure() · pymc-devs/pymc4@4f5382f · GitHub
• fix lint errors · pymc-devs/pymc4@e06d946 · GitHub
• add tests · pymc-devs/pymc4@ca9f334 · GitHub

minor fixes - c27e97f

• pymc4/inference_sampling_sample.py
  • dictのfor文の処理で不要な変数を減らしているだけです.
• pymc4/model/base.py
  • target_log_prob_fnメソッド: この辺りは変更は多いがやっていることは変わらないため省略
  • unobservedプロパティ: 順序付きdictの初期化の処理を修正
    • unobserved = {} とdictを初期化するとpython3.5以下では順序無しのdictになる点を修正
    • see: Get mcmc sampling to work by sharanry · Pull Request #9 · pymc-devs/pymc4 · GitHub

add some tests - a8e3dae

• 実際にはテストの追加ではなくリファクタリングをしているだけなので省略

add tests and fix model.configure() - 4f5382f

次のようなテストが複数追加されている.

• Modelクラスのインスタンス化
• @model.define で ed.Normal(0., 1., name='normal') を生成したModelインスタンスにセット
• assertでModelインスタンスのプロパティをチェック

    model = pm.Model()
    @model.define
    def simple(cfg):
        ed.Normal(0., 1., name='normal')
    
    assert len(model.variables) == 1
    assert len(model.unobserved) == 1
    assert "normal" in model.variables

fix lint errors - e06d946

lint対応なので省略: #pylint: disable=unused-argument でlintをdisableしている箇所もある.

add tests - ca9f334

テストが複数追加されている.

例えば, 次のテストは（対数）確率関数の値を近似的（approx）にチェックしている.

def test_model_log_prob_fn():
    model = pm.Model()
    @model.define
    def simple(cfg):
        mu = ed.Normal(0., 1., name="mu")
    
    log_prob_fn = model.target_log_prob_fn()
    with tf.Session():
        assert -0.91893853 == pytest.approx(log_prob_fn(0).eval(), 0.00001)

念のためにscipy.statsのdistributionモジュールで値をチェックしたみた.

from scipy.stats import norm
print(norm(0, 1).logpdf(0)) # => -0.9189385332046727

参考資料

• Get mcmc sampling to work by sharanry · Pull Request #9 · pymc-devs/pymc4 · GitHub