📂機械学習

중요도 샘플링

概要¹

重要度サンプリングはモンテカルロ法の一つで、有限和で積分（期待値）を近似するときに使える一種のサンプリングトリックだ。

ビルドアップ

標準正規分布で$z$値が$4$以上である確率は約$3.167 \times 10^{-5}$だ。

$$\begin{equation} \int_{4}^{\infty} \dfrac{1}{\sqrt{2\pi}} e^{-z^{2}/2}dz \approx 3.167 \times 10^{-5} \end{equation}$$

この積分をモンテカルロ積分で計算するジュリアコードは次のとおりだ。

using Distributions

n = 50000                      # sample size
z₀ = 4.0

for i in 1:5
    z = rand(Normal(), n)      # n-samples drawn from a standard normal distribution 
    est_MC = sum(z .> z₀)/n    # estmiates using Monte Carlo
    println(est_MC)
end

4.0e-5
6.0e-5
0.0
4.0e-5
4.0e-5

5回程度やってみると計算された値が真の値に近くなく、時には$0$として計算されることもある。この不正確さは標準正規分布の裾が非常に薄いために生じる。紙上で積分をするときには裾がどんなに薄くてもその計算が積分に含まれるが、コンピュータでサンプリングする際には有限なサンプルしか扱わないため、発生確率が非常に低い事象は計算に含まれないことがある。もちろん$n$を増やせばより正確にはなるが、これを解決方法とは呼ばない。

$$\int f(x)p(x) dx = \int f(x)\dfrac{p(x)}{q(x)} q(x) dx = \int \left[ f(x)\dfrac{p(x)}{q(x)} \right] q(x) dx$$

次に確率密度関数が$p$である確率変数$X$に対して、$f(X)$の期待値を計算する状況を考えてみよう。これは上の式の左辺と同様に計算できるが、ここで任意の$q(x)$を掛けて割っても依然として同じ値だ。式の右辺は左辺と同じ値を持つが、意味は「確率密度関数が$q$である確率変数$X$に対する$f(X)\dfrac{p(X)}{q(X)}$の期待値」になる。ここで$q$を$p$より裾が厚い確率密度関数として選べば、発生確率が非常に低い事象による不正確さを減らすことができる。

定義

$X$を確率密度関数が$p$である確率変数とする。$q$を次を満たすもう一つの確率密度関数とする。

$$\supp p \subset \supp q \qquad (\iff q\ne 0 \text{ if } p \ne 0)$$

重要度サンプリング^{importance sampling}とは、$f(X)$の期待値$E_{p}[f(X)]$を$E_{q}[f(X)\frac{p(X)}{q(X)}]$により計算して求める方法のことをいう。このとき$\supp$は関数のサポートである。

$$E_{p}[f] = \int f(x)p(x) dx = \int \left[ f(x)\dfrac{p(x)}{q(x)} \right] q(x) dx = E_{q}[fp/q]$$

ここで$q$を提案分布^{proposal distribution}、$w = \frac{p}{q}$を重要度重み^{importance weight}という。

説明

つまり重要度サンプリングとはモンテカルロ積分を次のように計算することである。

$$\dfrac{1}{n}\sum f(x_{i})\dfrac{p(x_{i})}{q(x_{i})}, X_{i} \sim q(x) \quad \left( \text{instead of } \dfrac{1}{n}\sum f(x_{i}), X_{i} \sim p(x) \right)$$

次のように使うと役に立つ。

• ビルドアップの例のように、$p(x)$の裾が非常に薄い場合、サポートがより大きい分布を提案関数とする。
• $p(x)$の分布でサンプリングするのが難しいとき、サンプリングしやすい分布を提案関数として$q(x)$に置く。

$q$が$p$の分母として乗じられるので、$\supp q$は$\supp p$を含まなければならない。

ビルドアップの積分$(1)$をコーシー分布を提案分布として重要度サンプリングで再計算した結果は次の通りだ。

# Julia code

using Distributions

n = 50000                      # sample size
z₀ = 4.0

for i in 1:5
    z = rand(Cauchy(), n)                       # n-samples drawn from a Cauchy distribution 
    w = pdf.(Normal(), z) ./ pdf.(Cauchy(), z)
    est_IS = sum(w[z .> z₀])/n                  # estmiates using Importance Sampling
    println(est_IS)
end

3.0330744703037005e-5
3.1808448538011614e-5
3.049259883894785e-5
3.305423127141489e-5
3.1807695818207125e-5

性質

不偏推定量

$\dfrac{1}{n}\sum\limits f(x_{i})\dfrac{p(x_{i})}{q(x_{i})}$は$E_{p}[f(X)]$の不偏推定量だ。

証明

$$\begin{align*} E_{q}\left[ \frac{1}{n}\sum_{i}^{n}f(X_{i})\dfrac{p(X_{i})}{q(X_{i})} \right] &= \frac{1}{n}\sum_{i}^{n} E_{q}\left[ f(X_{i})\dfrac{p(X_{i})}{q(X_{i})} \right] &\text{by linearity of$$$E$} \\ &= E_{q}\left[ f(X)\dfrac{p(X)}{q(X)} \right] \\ &= E_{p}\left[ f(X) \right] \end{align*}

■

분산

중요도샘플링에서의 분산은

$$\begin{align} \Var \left[ \frac{fp}{q} \right] &= E_{q}\left[ \left( \frac{fp}{q} \right)^{2} \right] - \left( E_{q}\left[ \frac{fp}{q} \right] \right)^{2} \\ &= \int \frac{f(x)^{2}p(x)^{2}}{q(x)} dx - \left( \int f(x)p(x) dx \right)^{2} \nonumber \end{align}$$

이므로, 유한한 분산을 갖기 위해선 다음이 성립해야한다.

$$\int \frac{f(x)^{2}p(x)^{2}}{q(x)} dx \lt \infty$$

これは$p(x)/q(x)$が有界であるということであり、これは再び$q(x)$の裾が$p(x)$の裾より厚くなければならないということと同じだ。分散を最小化する$q$は次の通りだ。

$$q^{\ast} = \argmin\limits_{q} \Var [fp/q] = \frac{ \left| f(x) \right| p(x)}{ \int \left| f(x) \right| p(x) dx}$$

証明

$(2)$の第二項は$q$に依存しているように見えるが、▽eq41◁であるため$q$とは無関係な値だ。

$$\argmin\limits_{q} \Var [fp/q] = \argmin\limits_{q} E_{q}\left[ \left( \frac{fp}{q} \right)^{2} \right]$$

イェンセンの不等式

開区間$I$ で関数$\phi$ がコンベックスし、二回微分可能、確率変数$X$ の期待値$\mu$ が存在し、 $X \subset I$ ならば $$\phi [ E(X) ] \le E [ \phi (X)]$$

しかし二次関数はコンベックスなので、

$$\begin{align*} && \left( E_{q}\left[ \frac{fp}{q} \right] \right)^{2} &\le E_{q}\left[ \left( \frac{fp}{q} \right)^{2} \right] \\ \implies && \left( \int f(x)p(x) dx \right)^{2} &\le \int \frac{f(x)^{2} p(x)^{2}}{q(x)} dx \end{align*}$$

ここで$q(x) = \frac{ \left| f(x) \right| p(x)}{ \int \left| f(x) \right| p(x) dx}$の場合、等号が成立する。

■