2017年12月26日 | ブログ記事

RustでMCMC(Metropolis-Hasting)

David

皆さんこんばんわ。Davidです。
今(12/26 01:01)この記事を書き始めたワケですが、実は12/24の投稿だったんですね〜。
いやホント申し訳ない。なぜかクリスマス・イブなら僕はassignされてないだろうとか言う妙な思い込みがあったんですね。申し訳ない。

はい、という訳で、一応書いていくんですけども、今回はRustでMCMCの一手法である、Metropolis-Hasting(M-H)法を書きました。

皆さんMCMCってご存知ですか?
もちろん、MiッCuMi〜Cuにし〜てやんよ〜のMCMCではありません。
Marcov Chain Monte CarloのMCMCです。

さて、すでに講義で数値計算を履修した人なら知ってることですが、Monte Carlo(法)というのは、何か計算したいものを、乱数を使ってゴリゴリ近似していく計算方法です。オーソドックスなのだと、円周率の計算をr1U(0,1),r2U(0,1)r_1 \sim U(0, 1), r_2 \sim U(0, 1)を満たすようなサンプルから計算するやつとかですかね!

MCMCもMonte Carlo法の一種です。ではこれは一体、何を計算するためのモノなのでしょうか?
それはある確率分布PPに従うサンプルを計算するための技法です。

一つ例を上げましょう。確率分布に関しては大丈夫ですね?これはナイーブには、あるドメインDがあって、そのドメインの要素それぞれに対して実数(確率)を返す関数です。確率なので、全ての要素に対応する確率の総和は1.0となります。型を書くと、P::DRP::D \to Rとなります.

今、我々にはある確率分布PPが与えられており、それを使い任意の要素xDx \in Dに対する確率または正規化されていない確率は計算できるとしましょう。問題は、&P&に従うようなサンプルxの多重集合を計算しなさい、というものになります。

これってどうやって計算するんでしょうか?

逆関数を求めて、確率からサンプルを求めればいいじゃないか!というのはあります。しかし厳密には、全単射じゃないと逆関数が定義できないので、そこからだけでも一部の分布にしかこの考え方による手法が適用できないことわかります。
また、実際に逆関数が存在することと、それを構成する難しさは別で、U(0,1)U(0, 1)に従うサンプルを求めたい分布に従うサンプルに写す関数を、確率分布の定理を使い積分をかますことで構成することができます。が、これは一般的な変な形をした複雑な分布に対しては、解析に求めることは難しすぎて、できません。

そこで、我々にはもっと一般的なサンプリング手法が必要になるのですが、一変数の分布に対しては効率的にサンプリングできるアルゴリズム(棄却サンプリング等)や、サンプルを得る訳ではなくサンプルを使って計算するアルゴリズム(重点サンプリング(サンプリングってついてるのにちょっと違うんだなこれが))はあるのですが、サンプルが多次元の場合にそれらの手法を適用しようとすると、色々問題が出てきます。まぁほぼ次元の呪いが悪いんですけどね。

なので、一般的な高次元のサンプルをサンプリングするアルゴリズムが必要になる訳ですが、その中で一つ大きな分野がMCMCになります。

MCMCはその名前の通り、Monte Carlo法なので乱数を使います。そしてまた、その名前の通り、Marcov Chain、マルコフ連鎖が裏の理論に大きく、重要な仕事をしています(なんかこの書き方英語っぽくてクサイですね)。

さて、今サンプリングした対象の分布をPPとすると、MCMCでは、このPPと共に、独立に選ぶことができる別の分布Qを使います。MCMCにおけるサンプリング過程の大雑把なキモを説明すると、各サンプルをPPと遷移確率分布QQで計算される分布のマルコフ連鎖の系列とすることで、効率的にサンプリングを行う、というモノです。

なぜマルコフ連鎖が絡んでくるのでしょう?こいつは出来の悪い、2gramベースでの言語生成モデルにしか出てこない、HMMを説明するための概念上あったほうが嬉しい、くらいの道具では無かったのでしょうか(まぁ言語モデルはNNがさいっっきょうなのでHMMとかくらいじゃね...iHMMとかiTHMMとかのディリクレ過程系を持ってこないと脳筋さいつよNNにはちょっと勝てないからね)?

これを説明するには、不変分布と詳細釣り合い条件、そしてエルゴード性という三つの要素について話す必要があるのですが、長いので、以下に概要を示します。

  1. 不変分布: あるマルコフ連鎖に対して定義される。マルコフ連鎖のあるステップでサンプリングされたサンプルxxが分布PPに従うなら、次のステップでサンプリングされるxx'も分布PPに従うような、マルコフ連鎖が持つ。これ持っているマルコフ連鎖は不変であるとか、定常であるとか言う。
  2. 詳細釣り合い条件: これを満たすようなPPQQがあれば、(1)を持つマルコフ連鎖が作れる.
  3. 任意の分布から始めたマルコフ連鎖が不変分布を持つことができる.

(1)はそれはそうって感じですね、PPからサンプリングしたいわけですから、系列全体をサンプルの多重集合とするなら、系列のステップごとにPPに従うようなサンプルがサンプリングされて欲しいですよね。
(2)は結構面白くて、これを満たせれば、そのマルコフ連鎖は不変になります。これは後々の手法が理論的に不変分布を持つための証明の随所で出てくるんですよ〜。
(3)もそれはそう。だって最初のサンプル点をどこに打てばわからないから、大体のマルコフ連鎖の系列の最初のステップの分布はPPじゃなくておそらく裾の広い一様な分布ですよね。そこからご所望のPPをマルコフ連鎖に持たせないといけない訳ですから。もっと感覚的には、サンプリング点の空間の一部にさえ近づきさえすれば、そこからいい感じに効率的にサンプリングできるので、最初の開始点がクソでも任意の点に動くことができれば、どうにかなるってことですかね。

はい。ここまで話しましたが、じゃあQってどうやって作るのさ?という話になりますね。そんな都合のいい分布そうそう無いんじゃ無イカ(注意: 筆者はSONY SIEに就職なのでイカはやってません。モンハンをやれ。プレステ最高FOOOOOO↑。なお同期達と人事は懇親会の紹介で皆イカを愛遊していることで盛り上がってたので特にそういう制約はありません)?という疑問は出てくるでしょうが、MCMCのすごいのは、任意のQQ'(解析的にサンプリングできて、かつ適当な遷移確率となりうる分布、例えば現在のサンプル点を平均とした標準偏差hogehogeの多次元正規分布)から、(2)を満たすQQを作るアルゴリズムが存在する、ということです。

今回私が実装したのは、M-H法です。

大まかには(眠くなってきた)QQ'から、PPのサンプリング点、の候補点xをサンプリングします。その後、PPQQ'でxをPのサンプリング点として受理する確率を計算し、それに従って受理します。受理された場合は現在の点を更新し、次の点に映りますが、されなかったら現在の点はそのままでxをQQ'からサンプリングし、受理確率を計算し...をやり直します。

さて、流石に寝るのですが、結果を乗せて起きましょう。分散が500の10次元等方正規分布から70000点サンプリングしたサンプルの1次元目の要素のヒストグラム(赤)です。今回は等方なので、多次元ガウス分布に関する一般的な性質から、対応する周辺分布(ここではサンプルの1次元目に相当する変数が従う分布)は分散500に従います。それを確かめるために、すでに知られている計算することができる1変数ガウス分布のグラフ(青)を重ねて書きました。しっかりと合ってますね。hist

リポジトリを張っておくので、コードはこっから見て下さい。結構なRustっぽくなさなのでバンバン指摘くれると嬉しいです!

ちなみに今回Rustで書いたのは特に意味はなく、Rustの勉強のためです。ヒストグラムとか全部自分で書いたけど、結構書きやすくていいですねRust。Haskellみたいな気持ちで、より早い(というかC/C++と同じくらい?)のでガンガン機械学習に使えて良いですね。なおRustに関する記述はここの部分とタイトルだけです(は?)。

明日は誰がAdv書くんでしたっけ?
traP人が多すぎてよく覚えていないですね...修士の学位もらって終了したら、僕もOBか...なんか寂しくなりますね!

オンオンオン!!!!!!!!!

ではみなさん、良いお年を!Merry X'smas!過好新年!

References

PRML上
計算統計Ⅱ

この記事を書いた人
David

この記事をシェア

このエントリーをはてなブックマークに追加

関連する記事

2017年12月8日
Jupyter Notebookと機械学習
hatasa-y
2017年11月14日
IBIS2017参加報告
Keijan
2017年11月4日
文章をよしなに分散表現しよう
David
2018年11月9日
機械学習の初歩についてのお話[アドベントカレンダー2018 16日目]
drun
2018年9月15日
論文の紹介【Universal Transformers】
hukuda222
2018年7月7日
ますにかえる ハッカソン参加記
yusuke

活動の紹介

カテゴリ

タグ