サンプリングのテスト追加と微修正[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