Isa@Diary

ソフトウェア開発やってます。プログラミングとか、US生活とかについて書きます。

Hamiltonian Monte-Carloを実装してみる

前回の記事でHMCの概要は理解したので作ってみることにした。 パラメタ推定のコードをGistに置いた: 1d_Gaussian_mixture_hmc · GitHub

Normal distribution

まずは普通の正規分布からサンプルすることを目標にする。各サンプル回ごとのleap frogによる移動を描画するとこんな感じ

f:id:isa_rentacs:20200802225104p:plain
サンプルごとの軌跡
水色が開始地点で紫が終了地点、番号はサンプルのindex。採択率が高いので基本的に次ステップは前ステップの終了位置から運動量だけランダムサンプルされる(=垂直移動する)

pが正ならxが正方向に移動するので負方向に回っているように見え、上下端に行くほど運動量が大きいのでxの変化も大きい。

また、途中の点はleap frog中の点なので最初の点と最後の点以外はpとxの時間が半ステップ分ずれていて、そのため最初/最後の点は少し軌跡からずれているように見える。

Normal mixture distribution

次に2つの正規分布の和からサンプルしてみる。オレンジが真の分布で、青がサンプルされたもの。

f:id:isa_rentacs:20200802225526p:plain
混合正規分布からのサンプル
先ほどと同じ軌跡を描くとこんな感じ

f:id:isa_rentacs:20200802230828p:plain

y軸に位置エネルギーを取ったものを動画にした。

www.youtube.com

(動画にするのにはcelluloidを使った。matplotlibのwrapperで使いやすい。)

Parameter inference

最後に与えられたデータとモデルからパラメタ推定するケースもやってみた。

 p(x|\mu_1, \mu_2, \pi) = \pi\cdot\mathcal{N}(x_i|\mu_1, 1) + (1-\pi)\cdot\mathcal{N}(x_i|\mu_2, 1) から生成したデータ点から \mu_1, \mu_2, \piを推定するというもの。

\piの事前分布としてBeta(1,1)を使うとleap frog中に容易に[0,1]からはみ出してしまうのでBeta(2,2)にした。この辺りきちんと処理するの難しそう。

まとめ

HMCの実装をした。位置エネルギーの各パラメタごとの偏微分が必要なのでここを自動微分に任せるのかーという気分になった(PyMC3はTheano, PyMC4はTF Probabilityを使っている) 。GPU使うとどのくらい速くなるのかも気になる。

NUTSに関しては調べてある程度理解したので実装まではしない予定。これである程度PyMC/STAN使った推論が裏で何やってるか理解できた気がする。