モンテカルロ積分おかわり

次回から社内 PRML 読書会が 11章「サンプリング法」に突入予定。
昔(と言っても1年前)、自分が初読の時にあれこれ試した記録も参照してもらってて、光成さんからいろいろツッコミをもらう。

PRML 11章の重点サンプリングと SIR を試す: http://d.hatena.ne.jp/n_shuyo/20100506/sampling

中で、ガウス分布に対してモンテカルロ積分を試しているところで、1000000 サンプル使ってもたいした精度が出てない点について、「ちなみに PRML には「十分な精度で期待値を推定するには、実際には 10 ないし 20 個の独立なサンプルで十分であろう」と書いてある。そうかなあ。」と書いた部分についてちょい議論。

普通に数値積分するなら、区分求積法の眷属(シンプソン公式とかガウス・ルジャンドル積分公式とか)で高精度で計算できるわけで(ここらへんは光成さんの受け売り)、サンプリング法のメリットはどの辺にあるのか、というあたり。
で、つらつら考えると、やっぱり台の空間が非常に高次元になる場合だろう。それからベイジアンの文脈なら、積分したい事後分布はもっと尖った(分散の小さい)ものになっていることが期待できる。
上の記事で試してみた Gamma(3,1) はどのどちらにも当てはまってなさ過ぎるなあ、とちょっと反省。

というわけで、題材を変えてもう一度試してみた。
多次元でお手軽な分布ということなら、みんなのアイドル、ディリクレ分布をおいて他にないだろう。
で、すこしは事後分布っぽいあたりということで、 p(z_1,...,z_5) = Dir(11,12,13,14,15) なあんて代物を考えてみて、E[ln z_1] を求めてみることにする。

まずは、ディリクレ分布に対しては E[ln z_k] = ψ(z_k)-ψ(Σz_i) であることがわかっているので、R を使って計算すると、

> digamma(11)-digamma(65)
[1] -1.814923

が正解。

続いてサンプリング法で積分してみよう。
R でディリクレ分布からのサンプリングを行うには MCMCpack の rdirichlet 関数を使うのがお手頃。*1

> library(MCMCpack)
> (samples = rdirichlet(20, c(11,12,13,14,15)))
            [,1]      [,2]      [,3]      [,4]      [,5]
 [1,] 0.06962782 0.1964396 0.2100311 0.2828679 0.2410336
 [2,] 0.11101838 0.1857581 0.2561643 0.1971526 0.2499067
 [3,] 0.18754173 0.2378527 0.1843784 0.1550071 0.2352201
 [4,] 0.22508286 0.1967883 0.2203731 0.1952157 0.1625401
 [5,] 0.09536442 0.2609473 0.1176665 0.2677637 0.2582580
 [6,] 0.13672580 0.2050953 0.1625757 0.3161957 0.1794076
 [7,] 0.17322767 0.1842071 0.2635450 0.1338937 0.2451265
 [8,] 0.29949577 0.1657692 0.1580508 0.1595159 0.2171684
 [9,] 0.13084545 0.1696426 0.1417823 0.2664633 0.2912665
[10,] 0.12282590 0.1726242 0.3050201 0.2197265 0.1798033
[11,] 0.21913593 0.1877809 0.1814171 0.2420822 0.1695838
[12,] 0.14696155 0.1795109 0.1807812 0.1641738 0.3285725
[13,] 0.13350812 0.1736853 0.1638994 0.2741285 0.2547786
[14,] 0.21677257 0.1897497 0.1858936 0.2061213 0.2014628
[15,] 0.17997192 0.1637597 0.2500945 0.2043236 0.2018504
[16,] 0.16981303 0.1809795 0.2256946 0.2361699 0.1873430
[17,] 0.17809172 0.1597740 0.1883114 0.2332019 0.2406209
[18,] 0.19542136 0.1540445 0.2130940 0.2102640 0.2271763
[19,] 0.16000754 0.1808641 0.2284115 0.1535779 0.2771390
[20,] 0.19726227 0.2107288 0.1529770 0.1955782 0.2434536

さてこのサンプルの先頭 n 個を使って E[ln z_1] を求めるとこうなる。

> cumsum(log(samples[,1])) / 1:20
 [1] -2.664591 -2.431325 -2.178802 -2.006923 -2.075548
 [6] -2.061253 -2.017238 -1.915790 -1.928896 -1.945705
[11] -1.906828 -1.907725 -1.915868 -1.888228 -1.876677
[16] -1.870200 -1.861686 -1.848959 -1.848094 -1.836851

20個目くらいで、正解 -1.814923 にだいぶ近くなってきている！
喜んで、じゃあサンプル 50 個くらいだとかなり精度が高くなるのでは、とやってみると……

> mean(log(rdirichlet(50, c(11,12,13,14,15))[,1]))
[1] -1.858163
> mean(log(rdirichlet(50, c(11,12,13,14,15))[,1]))
[1] -1.778846
> mean(log(rdirichlet(50, c(11,12,13,14,15))[,1]))
[1] -1.816499
> mean(log(rdirichlet(50, c(11,12,13,14,15))[,1]))
[1] -1.772064
> mean(log(rdirichlet(50, c(11,12,13,14,15))[,1]))
[1] -1.844697
> mean(log(rdirichlet(50, c(11,12,13,14,15))[,1]))
[1] -1.809723
> mean(log(rdirichlet(50, c(11,12,13,14,15))[,1]))
[1] -1.884225
> mean(log(rdirichlet(50, c(11,12,13,14,15))[,1]))
[1] -1.713337

うーん、もっと尖ってて、もっと高次元でないとダメかあ。
Dir(101,...,110) とかだとどうだろう？

> alpha <- 101:110
> digamma(alpha[1]) - digamma(sum(alpha))
[1] -2.35066          # ←★正解★

> mean(log(rdirichlet(10, alpha)[,1]))
[1] -2.434258
> mean(log(rdirichlet(20, alpha)[,1]))
[1] -2.344468
> mean(log(rdirichlet(50, alpha)[,1]))
[1] -2.362748
> mean(log(rdirichlet(50, alpha)[,1]))
[1] -2.354506
> mean(log(rdirichlet(50, alpha)[,1]))
[1] -2.346873

だいぶ良くなった。
今度は分布をなまらせて Dir(1,...,10) でやってみる。

> alpha <- 1:10
> digamma(alpha[1]) - digamma(sum(alpha))
[1] -4.57543          # ←★正解★

> mean(log(rdirichlet(10, alpha)[,1]))
[1] -5.198124
> mean(log(rdirichlet(20, alpha)[,1]))
[1] -4.630841
> mean(log(rdirichlet(50, alpha)[,1]))
[1] -4.488265
> mean(log(rdirichlet(50, alpha)[,1]))
[1] -4.512672
> mean(log(rdirichlet(50, alpha)[,1]))
[1] -4.72877
> mean(log(rdirichlet(50, alpha)[,1]))
[1] -4.327961

やっぱり分布が尖ってないとダメ。

というわけで、PRML の「十分な精度で期待値を推定するには、実際には 10 ないし 20 個の独立なサンプルで十分であろう」(下巻 p238)という記述は、分布がそこそこ尖っている(台のほとんどで≒0)という場合の話で、高次元だと区分求積法のお仲間よりも計算量で(圧倒的に)優位になっていく、と読めば良さそう。

実を言うと、精度を上げたいと思っても、サンプルをそこそこ増やしたくらいではあまり効果はない。
というのも、推定量の分散が次の形をしていることから。