PRML 11章の重点サンプリングと SIR を試す

PRML 機械学習 統計 sampling R

PRML 11章の予習中。
p(z) = Gamma(3, 1) について、 E[ln z] を求めよう。


まずは「正解」。
PRML Appendix B を見ると、Gam(tau | a,b) に対して、E[ln tau] = ψ(a) - ln b とある(ψはディガンマ関数)。
R なら簡単に計算できる。

> digamma(3) - log(1)
[1] 0.9227843

サンプリングによる期待値の推論(モンテカルロ積分)

p(z) から独立にサンプリングした \bf{z}^{(l)} \; (l=1,...,L) を用いて、 \mathbb{E}[f] = \int f(\bf{z})p(\bf{z})d(\bf{z}) \approx \frac{1}{L}\sum_{l=1}^L f(\bf{z}^{(l)}) と求めるのがモンテカルロ積分(であってるよね？)。
R でのガンマ分布からのサンプリングは rgamma 関数を用いれば出来るので、p(z) = Gamma(3, 1) について E[ln z] を求めるには次のように書けばいい。

E_ln_z <- function(L) sum(log(rgamma(L, 3, 1))) / L;

もちろんサンプルの個数 L を大きくすればするほど精度が良くなる。

> E_ln_z(10);
[1] 1.080437
> E_ln_z(20);
[1] 0.8831613
> E_ln_z(30);
[1] 0.8005687
> E_ln_z(50);
[1] 0.9022031
> E_ln_z(100);
[1] 0.8526046
> E_ln_z(200);
[1] 0.9655233
> E_ln_z(1000);
[1] 0.9480774
> E_ln_z(10000);
[1] 0.927171
> E_ln_z(100000);
[1] 0.9216669
> E_ln_z(1000000);
[1] 0.9228172

ちなみに PRML には「十分な精度で期待値を推定するには、実際には 10 ないし 20 個の独立なサンプルで十分であろう」と書いてある。そうかなあ。

重点サンプリング(importance sampleing)

rgamma() 関数が使えない、つまり「容易に p(z) からサンプリングできない」状況を考える。
そう言う場合でも、適当な提案分布を仮定して、PRML 11.1.4 の重点サンプリングで求められるはず。


p(z) からサンプリングするのが難しいとき、提案分布 q(z) と、q(z) から独立にサンプリングした \bf{z}^{(l)} \; (l=1,...,L) を用いて、 \mathbb{E}[f] \approx \sum_{l=1}^L w_l f(\bf{z}^{(l)}) と求めるのが重点サンプリング。
ここで w_l=\frac{p(\bf{z}^{(l)})/q(\bf{z}^{(l)})}{\sum_m p(\bf{z}^{(m)})/q(\bf{z}^{(m)})} であり、p(z), q(z) は正規化できてなくてもいいというのが魅力。


提案分布 q(z) は p(z) に近くて、簡単にサンプリングできるものを選ぶ。
最初コーシー分布にしようかと思ったけど(一様分布を tan() で変換すればサンプリングできる)、ガンマ分布の定義域は z>0 なので、指数分布を使うとしよう。
どちらも PRML 11.1.1 に載っている。

q(z) = \lambda \exp(-\lambda z)

q(z) からサンプリングするには、(0,1) から一様乱数 u をとり、-ln(1-u)/λ を求めればよい。
指数分布のパラメータであるλには、とりあえず 0.3 を選んでおこう。


さっくり実装。

importance <- function(L, lambda=0.3) {
	z_l <- -log(1 - runif(L)) / lambda;
	w_l <- dgamma(z_l, 3, 1) / (lambda * exp(-lambda * z_l));
	sum(w_l / sum(w_l) * log(z_l));
}

計算。

> importance(10);
[1] 0.5703455
> importance(20);
[1] 1.132268
> importance(100);
[1] 0.9373564
> importance(200);
[1] 0.8507546
> importance(1000);
[1] 0.932799
> importance(10000);
[1] 0.9362712
> importance(100000);
[1] 0.91805
> importance(1000000);
[1] 0.9231349
> importance(10000000);
[1] 0.9227575

SIR (sampling importance resampling)

重点サンプリングをやったから、ついでに SIR もやっておく。


重点サンプリングに出てきてた重みを使ってリサンプリングするのが SIR。
だが実のところ、積分近似値を求めるだけならそのまま重点サンプリングすればいい(PRML にもそう書いてある)。が、あえて SIR で積分してみよう。
提案分布は同じ λ=0.3 の指数分布。

sir <- function(L, lambda=0.3) {
	z_l <- -log(1 - runif(L)) / lambda;  # sampling
	w_l <- dgamma(z_l, 3, 1) / (lambda * exp(-lambda * z_l));  # weight
	z_l <- sample(z_l, L, replace=T, prob=w_l);  # resampling
	sum(log(z_l)) / L;
}

結果。
明らかに重点サンプリングより悪い。処理量も多くなるので、モーメントを求めるのに SIR はいいところなしだな。

> sir(10)
[1] 0.4273788
> sir(100)
[1] 0.939323
> sir(1000)
[1] 0.899571
> sir(10000)
[1] 0.9253118
> sir(100000)
[1] 0.9269862
> sir(1000000)
[1] 0.9237318
> sir(10000000)
[1] 0.9230488

他の本だと、 SIR のリサンプリングの個数はサンプリングの個数の 1/20 とかにするように書かれている。
実際、L 個から L 個サンプリングしたら同じ値が多数抽出されてしまうので、直感的にも良くなさそうなのだが、PRML はどうしてそうなっていないのだろう。

評価

ざっくり雰囲気がわかったところで、まじめに精度を評価してみよう。
提案分布をいろいろ変えて、積分値を 1000 回求め、平均、分散、L1 & L2 誤差を求めた。
提案分布にはガンマ分布 G(a,1) (a=1〜10)、指数分布 (λ=0.1, 0.2, ..., 1.0)、一様分布を用意してみた。
ガンマ分布 G(3,1) を選んだ場合は、重点サンプリングは \frac{1}{L}\sum_{l=1}^L f(\bf{z}^{(l)}) と一致する。
一様分布は定めた区間における数値積分と同等の精度になると思っていいのかな？

f <- log;
p <- function(z) dgamma(z, 3, 1);
exact_value <- digamma(3) - log(1);

# 重点サンプリング。提案分布も与えられるように。
importance <- function(f, p, L, dproposal, rproposal) {
	z_l <- rproposal(L);
	w_l <- p(z_l) / dproposal(z_l);
	sum(w_l * f(z_l)) / sum(w_l);
}
report <- function(name, results) {
	cat(sprintf("%s, mean=%f, sd=%f, L1=%f, L2=%f\n",
		name,
		mean(results), sd(results),
		sum(abs(results - exact_value)),
		sum((results - exact_value) ^ 2)
	));
}

L <- 10000;
N <- 1000;

# proposal = gamma distribution
for(a in seq(1,10)) {
	results <- sapply(1:N, function(i){
		importance(f, p, L, 
			function(z) { dgamma(z, a, 1) },
			function(L) { rgamma(L, a, 1) });
	})
	report(sprintf("gamma(%d,1)", a), results);
}

# proposal = exponential distribution
for(lambda in seq(0.1,1,0.1)) {
	results <- sapply(1:N, function(i){
		importance(f, p, L, 
			function(z) { lambda * exp(-lambda * z) },
			function(L) { -log(1 - runif(L)) / lambda });
	})
	report(sprintf("exp_dist(%.1f)", lambda), results);
}

# proposal = unified distribution
for(a in c(1,2,3,5,10,20,30,50,100)) {
	results <- sapply(1:N, function(i){
		importance(f, p, L, 
			function(z) { 1 / a },
			function(L) { runif(L, 0, a) });
	})
	report(sprintf("unif(0,%d)", a), results);
}

結果。

gamma(1,1), mean=0.922300, sd=0.018370, L1=14.571208, L2=0.337342
gamma(2,1), mean=0.922952, sd=0.007344, L1=5.811185, L2=0.053902
gamma(3,1), mean=0.922810, sd=0.006362, L1=5.058234, L2=0.040436
gamma(4,1), mean=0.923525, sd=0.011643, L1=9.274079, L2=0.135974
gamma(5,1), mean=0.923237, sd=0.035869, L1=25.980229, L2=1.285489
gamma(6,1), mean=0.926092, sd=0.098243, L1=61.730362, L2=9.652991
gamma(7,1), mean=0.966826, sd=0.113952, L1=94.781181, L2=14.911841
gamma(8,1), mean=1.014625, sd=0.167312, L1=154.870891, L2=36.399990
gamma(9,1), mean=1.095750, sd=0.178085, L1=222.254030, L2=61.599665
gamma(10,1), mean=1.154239, sd=0.215338, L1=287.018583, L2=99.895272
exp_dist(0.1), mean=0.923024, sd=0.007576, L1=6.016154, L2=0.057402
exp_dist(0.2), mean=0.922997, sd=0.006018, L1=4.826845, L2=0.036220
exp_dist(0.3), mean=0.922509, sd=0.005523, L1=4.387463, L2=0.030545
exp_dist(0.4), mean=0.922490, sd=0.005659, L1=4.488707, L2=0.032081
exp_dist(0.5), mean=0.922978, sd=0.006480, L1=5.103097, L2=0.041986
exp_dist(0.6), mean=0.922867, sd=0.007380, L1=5.957213, L2=0.054418
exp_dist(0.7), mean=0.922470, sd=0.008590, L1=6.832207, L2=0.073806
exp_dist(0.8), mean=0.922205, sd=0.010699, L1=8.574897, L2=0.114686
exp_dist(0.9), mean=0.922436, sd=0.013524, L1=10.652125, L2=0.182838
exp_dist(1.0), mean=0.923340, sd=0.018350, L1=14.579029, L2=0.336689
unif(0,1), mean=-0.403001, sd=0.003507, L1=1325.784979, L2=1757.718097
unif(0,2), mean=0.205763, sd=0.003619, L1=717.021016, L2=514.132220
unif(0,3), mean=0.512206, sd=0.004189, L1=410.578790, L2=168.592475
unif(0,5), mean=0.792722, sd=0.004868, L1=130.062036, L2=16.939811
unif(0,10), mean=0.918626, sd=0.006703, L1=6.319061, L2=0.062176
unif(0,20), mean=0.922781, sd=0.009517, L1=7.523265, L2=0.090476
unif(0,30), mean=0.922955, sd=0.011500, L1=9.166735, L2=0.132140
unif(0,50), mean=0.922383, sd=0.014572, L1=11.700022, L2=0.212299
unif(0,100), mean=0.922691, sd=0.020807, L1=16.530611, L2=0.432507

ガンマ分布の場合は p(z) = gamma(3,1) を提案分布に取ったものが一番精度がいいし、一様分布はこういう有界でない積分の場合に区間の選び方が難しいというあたりは直感に一致した結果が出ているが。
おもしろいのは、指数分布の λ=0.2〜0.4 では、gamma(3,1) を提案分布に取ったものより精度が良くなっていること。

赤が p(z)=Gamma(3,1)、青が ln(z)p(z)、緑が λ=0.3 の指数分布。
PRML には「重点サンプリングによるアプローチが成功するか否かも、サンプリング分布 q(z) が求めたい分布 p(z) にどの程度適合しているかに決定的に依存する」と書かれている。とはいえ、p(z) = q(z) が最も精度がいいとは限らないだろうなと思っていたが、グラフで見る限り、指数分布が精度がよいという理由が見あたらない。
こんな簡単な例であっても、良い提案分布を見つけるのは難しいという教訓、ということでいいんだろうか……