マルチカノニカル法で個数推定
はじめに
去年の11,12月にマルコフ連鎖モンテカルロ法(MCMC)の集中講義を受けてきました。 講義の中で紹介されたマルチカノニカル法を応用してレアイベントをサンプリングする手法が特に面白かったので、復習も兼ねて実装してみました。
マルチカノニカル法
まずレアイベントサンプリング用のマルチカノニカル法を説明します。
状態の確率分布がが与えられているときに、統計量の分布が知りたいとします。 単純にをサンプリングしてもいいのですが、統計量の周辺分布に極めて確率の小さい区間があるとその部分をうまくサンプリングすることができません。
この問題を回避するために適当な重みを用意して、の代わりに重み に従ってをサンプリングすることを考えます。 これは普通のメトロポリス・ヘイスティングス法で実現できます。 すなわち、状態から状態 に遷移する確率を
とすればいいわけです。
重みとして、求めたい周辺分布に近い分布の逆数を用いるのがマルチカノニカル法です。
重みに従ってをサンプリングしたのときの統計量 の周辺分布は
となって、この分布はほぼ平坦であることがわかります。 平坦な周辺分布からサンプリングすることで、普通にサンプリングすると確率0と推定されてしまうようなレアイベントにもMCMCで訪れられるのがポイントです。
Wang-Landau 法
残る問題は、求めたい周辺分布に近い分布をどのように得るかです。
周辺分布を大まかに求める方法としてWang-Landau 法というものがあります。 詳細はwikiの擬似コードを読むのが一番分かりやすいと思います(https://en.wikipedia.org/wiki/Wang_and_Landau_algorithm)。
大まかに言うと、のヒストグラムが平坦になるように出現回数の低いところの重みを逐次大きくしていき、ある程度平坦になったらより細かく重みを調整していく感じです。
演習問題
具体例として、次のような数え上げ問題を考えます。
L×Lの正方格子において、角と向かい側の角とをつなぐ自己交差のない経路の個数を求めよ。
文章だと分かりづらいですが、要は↓の動画で扱われている問題です。 www.youtube.com
一辺の長さがのとき、辺は全部で本あるので、各辺について使うか否かを考えると辺の配置は全部で通りです。探索する配置の数が指数関数で増大するので、愚直に全通りを試していてはすぐに計算が終わらなくなってしまいます。
ただ、この問題には効率的なアルゴリズムがあり、L=26までの厳密解が求まっているようです(http://oeis.org/A007764?language=japanese)。 なので個数推定をする意義はあまりないのですが、演習問題を解いているということでそこら辺は気にせずいきます。
この問題をマルチカノニカル法で扱うために、正方格子の各辺について使うか否か決めたもの (以下これを「状態 」と呼ぶ)に対してエネルギーを次式で定めます。
ここで、であることと状態が左上と右下の頂点をつなぐ自己交差のない経路であることが同値になっています。 先のマルチカノニカル法を用いてエネルギーの分布を求めて、に状態の総数をかければ今求めたい経路の数を推定することができます。
具体的な実装はこちら。 github.com
結果
以上の方法で自己交差のない経路の数を100回推定し、その平均と不偏分散の平方根を真値と比較したものが下図です。 有効数字2桁くらいの精度が得ることができました。
ちなみにL=6が手元のMacBook Air で4時間くらいかかっています。 L=6の状態の数がで、その中で自己交差のない経路になっているものはおよそ分の1しかないことを考えると、これほどレアな事象をそれなりの精度でサンプリングできるのはすごいと思えます。
L | 推定値 | 誤差 | 真値 |
---|---|---|---|
3 | 185.072 | 2.6228 | 184 |
4 | 8241.03 | 275.432 | 8512 |
5 | 1.2242e+06 | 75237.9 | 1262816 |
6 | 5.58901e+08 | 6.03206e+07 | 575780564 |
参考文献
http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0125062 https://sites.google.com/site/iwanamidatascience/vol2/number-place http://www.springer.com/us/book/9783642031625