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

WAICとWBICによる混合Gauss分布のモデル選択

Stan Advent Calendar 2017 - Qiitaの12/7の記事です。

Stanによる混合Gauss分布のパラメータ推定のモデル選択をWAIC及びWBICで行なってみます。

 

statmodeling.hatenablog.com

上記の記事で1次元の場合にされていることを今回は2次元でやってみました。

 

WAICやWBICについては上記記事や、渡辺先生による以下のページをご覧ください。

http://watanabe-www.math.dis.titech.ac.jp/users/swatanab/waicwbic.html

http://watanabe-www.math.dis.titech.ac.jp/users/swatanab/waicwbic_j.html

 

今回扱うデータを生成する分布は

- 平均は等しく、x方向とy方向で分散が異なる二つの正規分布からの混合分布

でやりたかったのですが、平均を等しくすると収束がうまくいかなかったので、少し平均をずらした二つの正規分布の混合分布としました。

実際には

- 平均(0.1, 0)で共分散行列が(2, 0 \\ 0, 0.1)の正規分布と平均(0, 0.1)で共分散行列が(0.1, 0\\ 0, 2)である正規分布の混合比率1:1の混合分布

でデータを生成しています。

 

これを

- model1:クラスタ数1、つまり単純な正規分布モデル

- model2:クラスタ数2の混合分布モデル、ただし混合比率は0.5で固定

の2通りで平均、分散を推定し、データを複数回生成した上でWAICとWBICの分布の様子、及びモデル選択の結果を調べます。

 

サンプル数としてはN=30, 100の2通りで実験してみました。

サンプル数30ではmodel2の収束がなかなかうまくいっていないことも多いです。

事前分布を設定することでうまく収束させることもできるかもしれませんが、今後試してみます。

 

WAIC, WBIC及びモデル選択の結果は以下の通りです。

まずサンプル数30で100回試行した結果

- WAICではmodel1が49回

- WBICではmodel1が43回

サンプル数100で50回試行した結果

- WAICではmodel1が12回

- WBICではmodel1が7回

選択されるという結果となりました。

 

それぞれの分布の様子は以下のようになりました。

f:id:unaoya:20171207231733p:plain

f:id:unaoya:20171207231759p:plain

 

今回用いたコードは以下のものです。WAICやWBICを計算する部分は、上記記事のものを参考にさせていただきました。

gist308edbde04b774e49287026a0a969b16

 

 

本記事を書くにあたり、上記記事に加えちょうど一年前の同じ日にStan Advent Calendarに書かれた以下の記事も参考にさせていただきました。

 

statmodeling.hatenablog.com

 

また混合分布のコードを書くにあたっては、stanのmanualも参考にしました。

 

本記事では実験した結果の報告のみですが、別のケースで実験したり理論的な部分も含めて今後は考察していきたいと思っています。