sigma

統計に関すること

Edwardをさわってみた

こちらはTensorFlow Advent Calendar 2017 - Qiita14日目の記事です。

 

今年何度かEdwardという単語を目にして、気になっていたのでさわってみることにしました。Edwardについては、以下のスライドで紹介されています。

確率的プログラミングライブラリEdward

公式ページはこちらです。

Edward – Home

tensorflowについても初心者なこともあって、公式のTutorialも難しいと感じました。

そこで以下のgithubのページのコード例を参考に二項分布のパラメータ推定を書いてみようというのがこの記事です。

edward/examples at master · blei-lab/edward · GitHub  

APIページを見ながら、Data, Model, Inference, Criticismという流れに沿ってコードを書いていきます。

 

まずはDataです。

Dataはnumpyのarrayやtensorflowのtensorで与えます。

今回はコイン投げ4回分の結果が与えられたものとします。表が1で裏が0ということにしましょう。

ここではデータを適当に手で与えます。

後のメソッドに渡すときに、データの形を合わせないとうまくいかなかったのでreshapeしてあります。

 

# DATA

n = 4
x_train = np.array([0, 1, 1, 1]).reshape((n,1))

 

次にModelを記述します。

ここでは確率変数を定義します。

確率変数xがサンプル数n=4で確率pの二項分布にしたがい、確率変数pはベータ分布Beta(1,1)に従うとしています。

pの方は事前分布ということでしょうか。

# MODEL

alpha = 1.0
beta = 1.0
p = Beta(concentration1=tf.Variable([alpha]),concentration0=tf.Variable([beta]))
x = Bernoulli(probs=p, sample_shape=(n))

  

次にInferenceです。Inferenceには大きく分けてVariational Inference, Monte Calro, exactの三種類あり、exactは共役分布を用いて解析的に事後分布を求めるものです。

今回の二項分布、ベータ分布のモデルでは共役分布を用いた推定及びMonte Carloを用いた推定が

edward/examples at master · blei-lab/edward · GitHub

のbeta_bernoulli.py, beta_bernoulli_conjugate.pyにあるので、Variational Inferenceを試してみました。

ここではパラメータの事後分布を表す確率変数を定義します。 今回はパラメータの確率変数pに対して事後分布の確率変数qpを定義するのですが、まずはMAP推定をするために一点のみ確率1でとる分布を設定します。

Inferenceのためのメソッドには、引数としてパラメータの確率変数と事後分布の確率変数、データの確率変数と実際のデータをdict形式で渡します。

データの確率変数xに対して実際のデータx_trainを与えます。

# INFERENCE by MAP

qp = PointMass(params=tf.Variable([0.1]))
inference = ed.MAP({p: qp}, data={x: x_train})
inference.run()

ここでqpにパラメータを設定しているのですが、これがどういうことなのかまだ理解できていません。 MAP推定なので事後分布の最大値を求めるわけですが、その過程で必要になるということでしょうか。

最後にCriticismです。 結果を見ます。

# CRITICISM

sess = ed.get_session()
print(sess.run(qp.mean()))
print(sess.run(qp.params))
print(qp.eval())

値はどれで見ても0.8124になります。 事前分布がBeta(1.0, 1.0)であればMAP推定量は今回だと0.75になるはずなので、何か違うのかもしれません。

x_post = ed.copy(x, {p: qp})
print(ed.evaluate('binary_accuracy', data={x_post: x_train}))

こちらは0.75になります。 事後分布で予測すると正解率が0.75になるということのようです。

x_trainを別のデータに置き換えると正解率は変わりますが、その都度乱数を発生しているということではなさそうです。   Variational Inferenceの別のメソッドとしてKallback Leibler divergenceを用いた変分推定も試してみます。

# INFERENCE by KLqp
alpha = 1.0
beta = 1.0
qp = Beta(concentration1=tf.Variable([alpha]),concentration0=tf.Variable([beta]))

inference = ed.KLqp({p: qp}, data={x: x_train})
inference.run()

こちらも結果を見ます。

# CRITICISM

sess = ed.get_session()
print(sess.run(qp.mean()))
print(sess.run(qp.concentration0), sess.run(qp.concentration1))
print(qp.eval())

出力は

[ 0.64722526]

[ 2.47322822] [ 4.53755856]

[ 0.91249388]

となります。

qp.eval()は乱数を発生しているようです。

単純に事後分布を計算するとBeta(2.0, 5.0)になるはずですが、それと比べると推定できていると言えるでしょうか。

x_post = ed.copy(x, {p: qp})
print(ed.evaluate('binary_accuracy', data={x_post: x_train}))

こちらは0.25となりました。

ちゃんとコードを動かすだけでもなかなか時間がかかりましたが、ひとまず動いたのでこれからも色々試していきたいと思います。

今回のコードは以下にあります。

gist53d988d72b7ce1ce964dba87f1d066c1