リサンプリング（復元抽出）で積分値を評価する

R 統計科学

の応用かつ、「これで間違ってねぇよな？」という自分への備忘録。インポータンスサンプリングや測度変換とは微妙に分母の形式が違ってくるので間違いやすいんだ。

数式の整理

以下の積分I
I = \int_1^{L} w(x) f(x) dx
を計算することを考える。積分範囲は以下の関数例にあうように適当に設定してしまった。

これはモンテカルロ法により以下のように近似される。UNIFORMは一様分布のこと。
I = (L-1)\int_1^{L} \frac{1}{L-1} w(x) f(x) dx \sim \frac{L-1}{N} \sum_{i=1}^{N} w(x_i)f(x_i), \,\,\,\,\, x_i \sim UNIFORM(1, L)

ここで、積分Iの式を、以下のように式変形して考える。
I = \frac{L-1}{N} sumw \sum_{i=1}^{N} \frac{w(x_i)}{sumw} f(x_i), \,\,\,\,\, sumw=\sum_{i=1}^{N} w(x_i)

ここで定義した
p(x_i)=\frac{w(x_i)}{sumw}
は各サンプルiに対する重み（確率）とみなすことができるので、この値に比例するようにリサンプリングすると
I = \frac{L-1}{N} sumw \frac{1}{N}\sum_{i=1}^{N} f(x_i)=\frac{L-1}{N} weight \sum_{i=1}^{N} f(x_i), \,\,\,\,\, x_i \sim p(x_i), \,\,\, weight=\frac{sumw}{N}
として計算することができる*1。


できる・・・はずだと信じてやってみる。

やってみる

#適当な関数w, f
w <- function(x){log(x)}
f <- function(x){exp(-x^2)}
#Lを4にしてどんなもんになるか描画



L <- 4
x <- seq(1, L, length.out = 100)
plot(x, w(x)*f(x))
#ふつーのモンテカルロ法（１万乱数）
N <- 10^4
x <- runif(N, min=1, max=L)
(L-1)*sum(w(x)*f(x))/N
#際サンプリングして計算
index <- sample(N, replace=TRUE, prob=w(x))
weight <- sum(w(x))/N
(L-1)*sum(weight*f(x[index]))/N
#実際に積分させてみた答え
stats::integrate(function(x){w(x)*f(x)}, 1, L)

実行結果の箇所だけ示すと

> (L-1)*sum(w(x)*f(x))/N
[1] 0.03627271
> (L-1)*sum(weight*f(x[index]))/N
[1] 0.03558152
> stats::integrate(function(x){w(x)*f(x)}, 1, L)
0.03588273 with absolute error < 4.1e-12

というわけで、ちゃんと一致してそうだ。

補記

ここで定義した
p(x_i)=\frac{w(x_i)}{sumw}
は各サンプルiに対する重み（確率）とみなすことができるので、この値に比例するようにリサンプリングすると

・・・の流れは
I = \int_1^{L} w(x) f(x) dx = (\int_1^{L} w(x)dx) \int_1^{L} \frac{w(x)}{\int_1^{L} w(x)dx} f(x) dx
=(\int_1^{L} w(x)dx) \frac{1}{N}\sum_i^N f(x_i), \,\,\,\, x_i \sim \frac{w(x)}{\int_1^{L} w(x)dx}
\sim \left( \frac{L-1}{N} \sum_i^N w(z_i)\right) \frac{1}{N}\sum_i^N f(x_i), \,\,\,\, x_i \sim \frac{w(x)}{\left( \frac{L-1}{N} \sum_i^N w(z_i)\right)}, z_i \sim UNIFORM(1, L)
\sim \frac{L-1}{N} sumw \frac{1}{N}\sum_i^N f(x_i), \,\,\,\, x_i \sim \frac{w(x)}{\frac{L-1}{N} sumw}, z_i \sim UNIFORM(1, L), sumw=\sum_i^N w(z_i)
\sim \frac{L-1}{N} weight \sum_i^N f(x_i), \,\,\,\, x_i \sim \frac{w(x)}{\frac{L-1}{N} sumw}, z_i \sim UNIFORM(1, L), sumw =\sum_i^N w(z_i), weight = \frac{sumw}{N}
とした方が見通しがよさげ

*1:ここで、"二回Nで割る"操作が入っているのを一回だと勘違いしてひどい目にあった。そのためのメモ