重点サンプリング

モンテカルロ積分

確率密度関数 p(θ) で表される確率分布に従う連続確率変数 θ を考える. θ についてのなんらかの関数 f(θ) の期待値計算は次式のように表される:

(1)E[f(θ)]=f(θ)p(θ)dθ

この期待値計算における積分を p(θ) の独立同分布なサンプル {θk}k=1N を用いた総和:

(2)1Nk=1Nf(θk)

によって近似する方法をモンテカルロ積分という.大数の法則により,N で式 (2) の総和は式 (1) の期待値に確率収束する.

重点サンプリング

上記のモンテカルロ積分をコンピュータで計算するには確率分布 p(θ) に従う乱数を生成する必要があるが,これは困難な場合が多い.

重点サンプリングはモンテカルロ積分における確率分布 p(θ) からの乱数生成を別の代理分布 q(θ) に肩代わりさせてモンテカルロ積分を評価する方法である. q(θ) を使うとモンテカルロ積分について次の恒等式が得られる:

(3)f(θ)p(θ)dθ=(f(θ)p(θ)q(θ))q(θ)dθ

つまり,『確率変数 f(θ) の分布 p(θ) に関する期待値』の代わりに『確率変数 f(θ)p(θ)q(θ) の分布 q(θ) に関する期待値』を評価しようというわけである.

(3) の右辺は q(θ) からの独立同分布なサンプル {θk}k=1N を用いた総和:

1Nk=1Nf(θk)p(θk)q(θk)

で近似できる.

代理分布 q(θ) として一様分布や標準正規分布といった乱数生成が容易な分布を用いれば分布 p(θ) の乱数生成を作らずにモンテカルロ積分を直接評価可能になる.

例:ラプラス分布の平均を計算する

ラプラス分布 p(x)=12exp(|x|) の平均の計算,つまり次の積分:

xp(x)dx

を重点サンプリングによるモンテカルロ積分によって計算することを考える.代理分布として標準正規分布:

q(x)=12πexp(x22)

を用いる.標準正規分布からの乱数はPythonであればNumpyの関数 np.random.randn を利用すれば容易に生成できる. Pythonに限らず,MATLABやJulia等の多くの言語でこのような randn 関数はサポートされている.

問題のモンテカルロ積分の計算をPythonで書くと以下のようになる:

import numpy as np


def p(x):
    return np.exp(-abs(x)) / 2


def q(x):
    return np.exp(-x ** 2 / 2) / np.sqrt(2 * np.pi)


N = 1000
x = np.random.randn(N)

print(np.sum(x * p(x) / q(x)) / N)

参考文献

  • 杉山将,『統計的機械学習 : 生成モデルに基づくパターン認識』,オーム社,2009.