インポータンス・サンプリング(importance sampling)の有難味について

なんとなーく頭でわかっているつもりになっていたものをちゃんと手を動かしてやっておきたい。
まず、インポータンスサンプリングについて一言で説明しておくと
「ある量の期待値を計算したいという状況において、その確率変数が大きな値を取る個所を重点的*1にサンプリングしてやることで、期待値評価の精度をあげる手法」
ということができる。

問題設定としては以下のような期待値評価をモンテカルロ法で計算したい、そんな状況を考える。
$E[f(x)] = \int_{-\infty}^{\infty} f(x) p(x) dx= \int_{-\infty}^{\infty} X e^{\frac{1}{2}(x-\mu)^2} e^{-x}dx$
$f(x) = X e^{\frac{1}{2}(x-\mu)^2}$
$p(x) = e^{-x}$
意味合いとしては「指数分布に従うxを使って正規分布（ガウシアン）的な関数の期待値を評価したい」ということ。
ここで $X$ , $\mu$ は適当に決める定数で、今回はそれぞれ10000, 6と設定した。

この積分自体は解析的に計算できて
$E[f(x)] = \sqrt{2\pi}X e^{-\mu + 0.5} \sim 102.4402$
となって、以下ではこの答えをきちんとモンテカルロ法で再現できるかどうかというのが問題になる。

まずはパス数を1,000パスとしてR言語を使って評価したい関数を定義する等しておく。

N <- 10^3
X <- 10^4
mu <- 6
#評価したい関数
f <- function(x, X, mu)
{
  X * exp(-0.5*(x - mu)^2) 
}

これを通常の素直なモンテカルロ法を使って計算すると

> x <- rexp(N)
> y <- f(x, X, mu)
> mean(y)
[1] 166.6111
> sd(y) / sqrt(N)
[1] 32.9838

となって標準誤差がかなり大きくなる。これは指数分布で発生させた乱数が評価すべき関数（コードだとfと書いたガウシアン的なもの）の極大値付近の値をなかなか（確率的に）取れずにいることに起因する。図示すると

(左軸が指数分布の確率密度関数値、右軸が評価したい関数f(x)の値）
となっており、期待値を正しく評価するためには赤丸で囲った領域付近を重点的に評価すべきなのだけれども、指数分布で生成した乱数がその辺(6辺り）の値を中々出してくれないのでうまく評価できていないということを意味している。

これを回避するためにインポータンスサンプリングを実施する。期待値計算を
$E[f(x)] = \int_{-\infty}^{\infty} f(x) p(x)dx= \int_{-\infty}^{\infty} f(x)\frac{p(x)}{\Pi(x)} \Pi(x) dx= \int_{-\infty}^{\infty} f(x) weight(x) \Pi(x) dx$
$weight(x)=\frac{p(x)}{\Pi(x)}$
と、変形し $\Pi$ を評価したい $f(x)$ と近い/似た関数形の確率密度関数として設定する。ここでは
$\Pi(x) = \frac{1}{\sqrt{2\pi}} e^{\frac{1}{2}(x-4)^2}$
のように平均だけ評価関数と異なる分散１の正規分布として設定した。

このようにインポータンスサンプリングのテクニックは

評価したい $f(x)$ と近い/似た関数形の確率密度関数 $\Pi(x)$ 探してくる
元の $p(x)$ (ここでは指数分布）ではなく $\Pi(x)$ からxをサンプリングする
式変形に伴いお釣りのウェイト $weight(x)=\frac{p(x)}{\Pi(x)}$ が出てくるのでこれを掛け算するのを忘れないようにして期待値計算をする

というステップで構成される。

このインポータンスサンプリングを使って再計算してみると

> nu <- 4
> x <- rnorm(N) + nu
> weight <- dexp(x) / dnorm(x - nu)
> y <- f(x, X, mu) * weight
> mean(y)
[1] 97.36879
> sd(y) / sqrt(N)
[1] 3.869766

となり、モンテカルロパス数は先程の通常のモンテカルロ法を使ったケースと同じである一方標準誤差が激減しており、より期待値を精緻に評価できることがわかる。

*1:インポータンスの所以