2018-10-15

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

pymc4 python

f:id:yukinagae:20171122095115p:plain

TL;DR

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

コミット

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

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

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つの処理を実行している.

サンプリング時の初期値を0.5で設定（乱数ではない）
tfp.mcmc （tensorflow_probabilityのMCMC実装）でサンプリング処理を定義
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))

各々の処理を見ていく.

サンプリング時の初期値を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)))

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))

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

2018-10-12

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

pymc4 python

f:id:yukinagae:20171122095115p:plain

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 処理の挙動が怪しい気がするので確認したい

参考資料

2018-10-11

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

pymc4 python

f:id:yukinagae:20171122095115p:plain

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

2018-10-10

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

pymc4 python

f:id:yukinagae:20171122095115p:plain

TL;DR

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

コミット

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

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

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

参考資料

2018-10-09

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

pymc4 python

f:id:yukinagae:20171122095115p:plain

TL;DR

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

コミット

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

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

2018-10-05

target_log_prob_fn: MCMCサンプリングの実装 [66f95ca, 8aaa0ff, cded7c5] - pymc4のソースコード読んでみた

pymc4 python

f:id:yukinagae:20171122095115p:plain

TL;DR

前回の log_prob_fn と同様に複数の未観測変数の対数確率の合計を計算している様子なので細かい点は省略

コミット

2018/06/21から2018/06/23の間のコミットです.

target_log_prob_fn

前回の log_prob_fn と同様に複数の未観測変数の対数確率の合計を計算しているように見えます. ただ今回の target_log_prob_fn では値ではなく関数を返している点が異なります.

コメントを読むと Unnormalized target density as a function of unobserved states. とあるため、未観測変数の正規化していない対数確率（正規化していないので正確には確率とは呼ばないはず）を計算しているはずです.

def target_log_prob_fn(self, *args, **kwargs):
    """
    Unnormalized target density as a function of unobserved states.
    """

    def log_joint_fn(*args, **kwargs):
        states = dict(zip(self.unobserved.keys(), args))
        states.update(self.observed)
        log_probs = []

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

        with ed.interception(interceptor):
            self._f(self._cfg)
        log_prob = sum(log_probs)
        return log_prob
    
    return log_joint_fn

TODO

[ ] states.update(self.observed) で観測変数も対象にしている理由が不明
[ ] log_prob_fn との使い分けが不明

参考資料

2018-10-04

log_prob_fn: MCMCサンプリングの実装 [66f95ca, 8aaa0ff, cded7c5] - pymc4のソースコード読んでみた

pymc4 python

f:id:yukinagae:20171122095115p:plain

TL;DR

複数の未観測変数の対数確率の合計を計算しています. これがMCMCサンプリングのコアの部分です.

コミット

2018/06/21から2018/06/23の間のコミットです.

以下ファイルが修正されています.

pymc4/model/base.py: 対数確率を計算するメソッドを2つ実装
- log_prob_fn <= 今回はここを読みます.
- target_log_prob_fn <= 次回
pymc4/inference_sampling_sample.py
- iteritems() の代わりに keys() を使うことで未使用の変数を減らしているだけです.

Model - model/base.py

log_prob_fn

複数の未観測変数の対数確率の合計を計算しています.

注意) x 引数は関数内で使用されていないので, 後のコミットで削除されています.

def log_prob_fn(self, x,  *args, **kwargs):
    logp = 0
    for i in self.unobserved.keys():
        logp += self.unobserved[i].rv.distribution.log_prob(value=kwargs.get(i))

    return logp

前提条件として, 次のように掛け算（x）は値のlogを取ることで足し算（+）に変換できます.

変数1の確率関数 x 変数2の確率関数 x 変数3の確率関数 x … x 変数nの確率関数 ↓
変数1の対数確率関数 + 変数2の対数確率関数 + 変数3の対数確率関数 + … x 変数nの対数確率関数