マルチカノニカル法で個数推定

はじめに

去年の11,12月にマルコフ連鎖モンテカルロ法(MCMC)の集中講義を受けてきました。 講義の中で紹介されたマルチカノニカル法を応用してレアイベントをサンプリングする手法が特に面白かったので、復習も兼ねて実装してみました。

マルチカノニカル法

まずレアイベントサンプリング用のマルチカノニカル法を説明します。

状態xの確率分布が\pi(x)が与えられているときに、統計量\xi (x)の分布が知りたいとします。 単純にxをサンプリングしてもいいのですが、統計量の周辺分布P(\xi)に極めて確率の小さい区間があるとその部分をうまくサンプリングすることができません。

この問題を回避するために適当な重みQ(\xi )を用意して、\pi (x)の代わりに重みQ(\xi (x)) \pi(x) に従ってxをサンプリングすることを考えます。 これは普通のメトロポリスヘイスティングス法で実現できます。 すなわち、状態xから状態x' に遷移する確率を

 \mathrm{min}(1, \frac{Q(\xi (x')) \pi(x') p(x \to x')}{Q(\xi (x)) \pi(x)
p(x' \to x)})

とすればいいわけです。

重みQ(\xi )として、求めたい周辺分布P(\xi)に近い分布\tilde{P}(\xi)の逆数を用いるのがマルチカノニカル法です。

重み \tilde{P}(\xi (x))^{-1} \pi (x)に従って xをサンプリングしたのときの統計量 \xi の周辺分布 P'(\xi )

 P'(\xi ) \propto P(\xi)  \tilde{P}(\xi )^{-1}

となって、この分布はほぼ平坦であることがわかります。 平坦な周辺分布 P'(\xi )からサンプリングすることで、普通にサンプリングすると確率0と推定されてしまうようなレアイベントにもMCMCで訪れられるのがポイントです。

Wang-Landau 法

残る問題は、求めたい周辺分布P(\xi)に近い分布\tilde{P}(\xi)をどのように得るかです。

周辺分布P(\xi)を大まかに求める方法としてWang-Landau 法というものがあります。 詳細はwiki擬似コードを読むのが一番分かりやすいと思います(https://en.wikipedia.org/wiki/Wang_and_Landau_algorithm)。

大まかに言うと、 \xiヒストグラムが平坦になるように出現回数の低いところの重みを逐次大きくしていき、ある程度平坦になったらより細かく重みを調整していく感じです。

演習問題

具体例として、次のような数え上げ問題を考えます。

L×Lの正方格子において、角と向かい側の角とをつなぐ自己交差のない経路の個数を求めよ。

文章だと分かりづらいですが、要は↓の動画で扱われている問題です。 www.youtube.com

一辺の長さがLのとき、辺は全部で2L(L+1)本あるので、各辺について使うか否かを考えると辺の配置は全部で2^{2L(L+1)}通りです。探索する配置の数が指数関数で増大するので、愚直に全通りを試していてはすぐに計算が終わらなくなってしまいます。

ただ、この問題には効率的なアルゴリズムがあり、L=26までの厳密解が求まっているようです(http://oeis.org/A007764?language=japanese)。 なので個数推定をする意義はあまりないのですが、演習問題を解いているということでそこら辺は気にせずいきます。

この問題をマルチカノニカル法で扱うために、正方格子の各辺について使うか否か決めたものC (以下これを「状態C 」と呼ぶ)に対してエネルギーを次式で定めます。

 \displaystyle {
E(C) = \sum_{i : \mbox{左上と右下以外の頂点}} |\mbox{(iの次数)} - 2| + \sum_{i : \mbox{左上と右下の頂点}} |\mbox{(iの次数)} - 1|
+ (\mbox{2個以上の頂点からなる連結要素の数)}
}

ここで、E(C) = 0であることと状態Cが左上と右下の頂点をつなぐ自己交差のない経路であることが同値になっています。 先のマルチカノニカル法を用いてエネルギーの分布P(E)を求めて、P(E=0)に状態の総数2^{2L(L+1)}をかければ今求めたい経路の数を推定することができます。

具体的な実装はこちら。 github.com

結果

以上の方法で自己交差のない経路の数を100回推定し、その平均と不偏分散の平方根を真値と比較したものが下図です。 有効数字2桁くらいの精度が得ることができました。

ちなみにL=6が手元のMacBook Air で4時間くらいかかっています。 L=6の状態の数が 2^{2 \times 6 \times 7} \sim 3 \times 10^{50}で、その中で自己交差のない経路になっているものはおよそ 6.4 \times 10^{42}分の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