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の対数確率関数

参考資料

2018-10-03

PyMC4のInstallation failsというissueに対応するPR送った

TL;DR

次のissueにあるように, 現状だと依存性の解決の部分でfailしてインストールできないのでとりあえずforkして dependency-resolution というbranchで修正してみた.

see: Installation fails · Issue #23 · pymc-devs/pymc4 · GitHub

修正点

これを

tf-nightly==1.9.0.dev20180607
tfp-nightly==0.3.0.dev20180725
tb-nightly==1.9.0a20180613

こうするだけの簡単なお仕事.

tf-nightly
tfp-nightly
tb-nightly

とりあえず本体にPR送った(｀・ω・´)

Specified the latest versions of tf/tfp/tfb-nightly by yukinagae · Pull Request #24 · pymc-devs/pymc4 · GitHub

おまけ

pip コマンドよくわかってないオプションがあったので軽く調べた.

pip install --user git+https://github.com/pymc-devs/pymc4.git#egg=pymc4

—user
- ユーザのローカルディレクトリ（ ~/.local/ など）に対象のパッケージをインストールする
egg=[プロジェクト名を明示的に指定]

また, 次のようにbranch名も指定できる * @[branch名]

pip install --user git+https://github.com/yukinagae/pymc4@dependecny-resolution.git#egg=pymc4

参考資料

pip install — pip 18.0 documentation

2018-10-02

Add model.target_log_prob_fn() sampling [a703c21] - pymc4のソースコード読んでみた

pymc4 python

f:id:yukinagae:20171122095115p:plain

TL;DR

ターゲットとなる unobserved（未観測なRandomVariableインスタンス） の対数確率の合計を返すメソッドを実装しています.

コミット

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

Add model.target_log_prob_fn() · pymc-devs/pymc4@a703c21 · GitHub

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

pymc4/model/base.py: Model クラスに次の２つのメソッドが追加
- target_log_prob_fn(self, args, kwargs)
- unobserved(self) *@property
pymc4/util/interceptors.py
- VariableDescription クラスに rv（RandomVariableのインスタンス） のプロパティが追加されただけです. この追加プロパティは target_log_prob_fn で使用されています（後述）.

target_log_prob_fn

unobserved（未観測なrv） の log_probability（対数確率） の合計を返します.

※ちなみにちょっとした罠ですが, 変数名の i というのはこのコードの実装では index ではなく, dictの key なので, 単に k などとした方がよいです.

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

unobserved

self.variables （RandomVariableインスタンスのリスト）から unobserved を抽出しているだけです（＝ observed ではないものを抽出）.

ちなみに, わざわざ collections.OrderedDict で順序を保持したdictionary配列にしている理由はこの時点ではよくわかりません.

    @property
    def unobserved(self):
        unobserved = {}
        for i in self.variables:
            if self.variables[i] not in self.observed.values():
                unobserved[i] = self.variables[i]
         unobserved = collections.OrderedDict(unobserved)
        return unobserved

参考資料

Get mcmc sampling to work by sharanry · Pull Request #9 · pymc-devs/pymc4 · GitHub
確率質量関数 - Wikipedia: probability mass function（PMF）- 離散値用の確率関数
確率密度関数 - Wikipedia: probability density function（PDF）- 連続値用の確率変数

2018-10-01

python3.5削除 + pep8対応 [bb4de21] - pymc4のソースコード読んでみた

pymc4 python

f:id:yukinagae:20171122095115p:plain

TL;DR

python3.5は対応せず、3.6以上対応の方針
pep8のコードスタイルの修正なので、特に重要な点はなさそうです。今後はlint系のコミットは冗長なので省略するかもしれません。

コミット

2018/06/11のコミットです。

pep8 · pymc-devs/pymc4@bb4de21 · GitHub

以下ファイルが修正・追加されていますが、基本的にpep8のコードスタイルの修正なので、説明は割愛します。

.travis.yml
pymc4/distributions/base.py
pymc4/inference_sampling_sample.py
requirements-dev.txt
setup.cfg
tests/test_model.py
tests/test_nothing.py

python3.5削除（3.6以上対応の方針）

.travis_yml を見ると、pythonのバージョン 3.5 が削除されています。

以下の議論を見ると、基本的に3.6対応していく流れみたいです。

see: Supported Python versions? · Issue #2 · pymc-devs/pymc4 · GitHub

3.6の新機能として、主に以下の f-strings literal と variable annotation があるみたいなので、これらを使用するなら3.5は捨てる必要があります（正確には annotation は3.5では単に無視されるみたいですが、 f-strings が使えないということです）

PEP 498, formatted string literals.
PEP 526, syntax for variable annotations.

Debian 9 (stable) では3.5が使われているようですが、結論としては miniconda使えばいいじゃん ということになったみたいです。

参考資料

2018-09-28

テストサンプルの生成 [e334115, d07338e, 93bc07b] - pymc4のソースコード読んでみた

pymc4 python

f:id:yukinagae:20171122095115p:plain

概要

Model クラスのサンプル生成のメソッドを読んでみます。
- test_point

コミット

2018/06/09から2018/06/11の間のコミットです。

Model - pymc4/model/base.py

test_point

まずはこの test_point メソッドがでどのように呼ばれるか調べます。pymc4のプロジェクト内では以下の tests/test_model.py のみで使用されています。

@pm.inline
def model(cfg):
    ed.Normal(0., 1., name='normal')

testval_random = model.test_point()
testval_mode = model.test_point(sample=False)

assert testval_mode['normal'] == 0.
assert testval_mode['normal'] != testval_random['normal']

ここでは2種類の呼び方があるみたいです。

test_point(): デフォルトでsample=True
test_point(sample=False)

事前情報として、ed.Normal(0., 1., name='normal') というのは、正規分布で平均が0かつ標準偏差が1という意味です。

test_point(sample=True)

Modelに設定されているRandomVariableのインスタンスの確率分布に沿って、1つの乱数が生成されます。

例) 正規分布で平均が0かつ標準偏差が1の確率分布から乱数を生成する

test_point(sample=False)

Modelに設定されているRandomVariableのインスタンスの確率分布に沿って、中央値が返されます。

例) 正規分布で平均が0かつ標準偏差が1の確率分布の中央値は0なので、常に0が返される

さっと test_point メソッドの実装も見てみます。

以下がそのメソッドの全容です。

def test_point(self, sample=True):
    # 1
    def not_observed(var, *args, **kwargs):
        return kwargs['name'] not in self.observed
    values_collector = interceptors.CollectVariables(filter=not_observed)
    chain = [values_collector]

    # 2
    if not sample:
        def get_mode(state, rv, *args, **kwargs):
            return rv.distribution.mode()
        chain.insert(0, interceptors.Generic(after=get_mode))

    # 3
    with self.graph.as_default(), ed.interception(interceptors.Chain(*chain)):
        self._f(self.cfg)

    # 4
    with self.session.as_default():
        returns = self.session.run(list(values_collector.result.values()))
    return dict(zip(values_collector.result.keys(), returns))

一つひとつバラバラに見ていきます。

この辺りの挙動はあまりよくわかっていないので飛ばします。

    def not_observed(var, *args, **kwargs):
        return kwargs['name'] not in self.observed
    values_collector = interceptors.CollectVariables(filter=not_observed)
    chain = [values_collector]

test_point(sample=False) の場合の処理です。やはり予想通り、中央値を返すようになっています。

    if not sample:
        def get_mode(state, rv, *args, **kwargs):
            return rv.distribution.mode()
        chain.insert(0, interceptors.Generic(after=get_mode))

以前、 edward2のinterception処理の記事で確認したので飛ばします。

    with self.graph.as_default(), ed.interception(interceptors.Chain(*chain)):
        self._f(self.cfg)

tensorflowのsessionを実行して、実際にGraph内の演算処理を実行しています。この辺りも今後もう少し理解できてくるはずです。

    with self.session.as_default():
        returns = self.session.run(list(values_collector.result.values()))
    return dict(zip(values_collector.result.keys(), returns))

The tf.Session.run method is the main mechanism for running a tf.Operation or evaluating a tf.Tensor. You can pass one or more tf.Operation or tf.Tensor objects to tf.Session.run, and TensorFlow will execute the operations that are needed to compute the result.

see: Graphs and Sessions | TensorFlow

このコミットを読み続けるのも疲れたので、次回から次のコミットに進みます。

TODO

[ ] #1 の not_observed の処理の内容確認

参考資料

2018-09-27

tensorflowのグラフ構造 [e334115, d07338e, 93bc07b] - pymc4のソースコード読んでみた

pymc4 python

f:id:yukinagae:20180927082722g:plain

概要

まずは Model クラスの初期化処理系のメソッドを読んでみます。
- _init_variables: 今回はここの self.graph.as_default() の処理を読みます

コミット

2018/06/09から2018/06/11の間のコミットです。

Model - pymc4/model/base.py

_init_variables

今回は _init_variables 内の self.graph.as_default() の処理を見ていきます。

info_collector = interceptors.CollectVariablesInfo() # info_collector() で呼べるcallable
with self.graph.as_default(), ed.interception(info_collector):
    self._f(self.cfg)

self._variables = info_collector.result

self.graph は Model クラス内で以下のように定義されているので、実際には self.session.graph == tf.session.graphと同等です。

@property
def graph(self):
    return self.session.graph

また、今回の読む対象に絞って [余計な情報を省略] + [変数を読み替える]、を行うと以下のコードと同等になります。

with tf.Session.graph().as_default():
    ed.Normal(0., 1., name='normal') <= ここはed.XXX()

上記の XX 部分は tensorflow_probability の distribution クラスをedward2の RandomVariable クラスでwrapしたものです。

see: probability/generated_random_variables.py at master · tensorflow/probability · GitHub

結局これらを理解するには tensorflow のGraphの動作を理解する方が早そうです。

f:id:yukinagae:20180927082722g:plain

see: Graphs and Sessions | TensorFlow

実際に as_default() の使用例を見てみると、以下のように with スコープで実行された tf.XXX の処理の結果がそのスコープのGraphに追加されていくようです。

g_1 = tf.Graph()
with g_1.as_default():
  # Operations created in this scope will be added to `g_1`.
  c = tf.constant("Node in g_1")

  # Sessions created in this scope will run operations from `g_1`.
  sess_1 = tf.Session()

g_2 = tf.Graph()
with g_2.as_default():
  # Operations created in this scope will be added to `g_2`.
  d = tf.constant("Node in g_2")

# Alternatively, you can pass a graph when constructing a <a href="../api_docs/python/tf/Session"><code>tf.Session</code></a>:
# `sess_2` will run operations from `g_2`.
sess_2 = tf.Session(graph=g_2)

assert c.graph is g_1
assert sess_1.graph is g_1

assert d.graph is g_2
assert sess_2.graph is g_2

つまり以下のコードを簡略的に理解すると、 g1(graph) -> c(constant) という紐付きをGraphとして構造化することができるということです。

with g_1.as_default():
  # Operations created in this scope will be added to `g_1`.
  c = tf.constant("Node in g_1")

同様に以下のコードも、tf.Session.graph() -> ed.Normal() と紐付けることができます。

with tf.Session.graph().as_default():
    ed.Normal(0., 1., name='normal') <= ここはed.XXX()