モンテカルロ積分

Pythonによるモンテカルロ法入門（2014/6/20）のつづき。3章のモンテカルロ積分について実験した。

モンテカルロ積分を使うと統計や機械学習で頻繁に出てくる期待値を求める積分が乱数生成で簡単に計算できる。期待値を求める積分とは下の形をした積分。

$\displaystyle E[f] = \int p(x) f(x) dx$

ここで、f(x) は任意の関数、p(x) は確率分布を表す。たとえば、 $f(x)=x$ のときは確率変数Xの平均、 $f(x)=(x-\bar{x})^2$ のときは確率変数Xの分散だった。実際、f(x) は上の二つに限らずどんな関数でもよい。

モンテカルロ積分のアルゴリズムは以下のとおり。確率分布から生成したサンプルの平均値で積分を近似するのがポイント。

確率分布 p(x) からサンプル X = [x_1, x_2, ..., x_N] を生成
E[f] の近似として $\displaystyle \frac{1}{N} \sum_{n=1}^N f(x_n)$ を計算

PRML（パターン認識と機械学習）の19ページにも下のような記述がある。

ある関数 f(x) の確率分布 p(x) の下での平均値を f(x) の期待値（expectation）と呼び、E[f] と書く。離散分布に対しては
$\displaystyle E[f] = \sum_{x} p(x) f(x)$
連続分布の場合は
$\displaystyle E[f] = \int p(x) f(x) dx$
どちらの場合も、確率分布や確率密度から得られた有限個のN点を用いて、期待値はこれらの点での有限和で近似できる。
$\displaystyle E[f] \simeq \frac{1}{N} \sum_{n=1}^N f(x_n)$
サンプリング法について議論する際にこの結果を頻繁に使うことになる。

具体例で本当に成り立つか検証してみよう。Rによるモンテカルロ法入門では練習問題3.1で正規・コーシーベイズ推定量の積分を求めろという難しい例からいきなり始まるけれど自分で検算できるもっと簡単な積分から始めよう。

簡単な積分計算例

関数 f(x) の区間 [a, b] の以下の定積分をモンテカルロ法で求めてみよう。

$\displaystyle I = \int_a^b f(x) dx$

この式は先の期待値を求める式と違って確率分布 p(x) が入っていない。入っていないなら入れてしまおう！区間 [a, b] の一様分布を p(x) とすると I は下のように変形できる。

$\displaystyle I = (b - a) \int_a^b p(x) dx \int_a^b f(x) dx = (b - a) \int_a^b p(x) f(x) dx \simeq (b - a) \frac{1}{N} \sum_{n=1}^N f(x_n)$
上の式変形は間違い！正しい式変形は↓（コメント参照）

$\displaystyle I = (b - a) \int_a^b f(x) \frac{1}{b - a} dx = (b - a) \int_a^b f(x) p(x) dx \simeq (b - a) \frac{1}{N} \sum_{n=1}^N f(x_n)$

ここで、一様分布の [a, b] での~~定積分~~確率密度が 1/(b-a) であることを利用している。モンテカルロ積分のサンプル生成は与えられた区間の一様分布から生成することになる。手元の公式集の定積分の例をいくつか確かめてみよう。scipy.integrateで計算した場合と答えが一致することも確認しておいた。

(1) $\displaystyle I = \int_{-1}^1 x^2 dx$

f:id:aidiary:20140728223506p:plain

# 任意の関数 f(x) の [a, b] での積分をモンテカルロ法で求める
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import uniform
import scipy.integrate

a, b = -1, 1
f = lambda x: x ** 2

# 被積分関数をプロット
x = np.linspace(-1.5, 1.5, 1000)
y = f(x)
plt.plot(x, y)

# 積分範囲を色付け
ix = arange(a, b, 0.001)
iy = f(ix)
verts = [(a, 0)] + list(zip(ix, iy)) + [(b,0)]
poly = Polygon(verts, facecolor='0.8', edgecolor='k')
gca().add_patch(poly)

# scipy.integrateで積分を計算
I = scipy.integrate.quad(f, a, b)[0]
print "scipy.integrate:", I

# モンテカルロ積分
N = 100000
x = uniform(loc=a, scale=b-a).rvs(size=N)
I = (b - a) * np.mean(f(x))
print u"モンテカルロ積分:", I

実行すると

scipy.integrate: 0.666666666667
モンテカルロ積分: 0.669414303329

(2) $\displaystyle I = \int_0^1 x e^{x^2} dx$

f:id:aidiary:20140728224724p:plain

a, b = 0, 1
f = lambda x: x * np.exp(x**2)
# 以下同

scipy.integrate: 0.85914091423
モンテカルロ積分: 0.855417623185

(3) $\displaystyle I = \int_0^1 x^2 e^{2x} dx$

f:id:aidiary:20140728223513p:plain

a, b = 0, 1
f = lambda x: x**2 * np.exp(2*x)
# 以下同

scipy.integrate: 1.59726402473
モンテカルロ積分: 1.60222532927

(4) $\displaystyle I = \int_1^e \frac{\log x}{x} dx$

f:id:aidiary:20140728223508p:plain

a, b = 1, np.exp(1)
f = lambda x: np.log(x) / x
# 以下同

scipy.integrate: 0.5
モンテカルロ積分: 0.499921834232

(5) $\displaystyle I = \int_1^2 x^3 \log x dx$

f:id:aidiary:20140728223515p:plain

a, b = 1, 2
f = lambda x: x**3 * np.log(x)
# 以下同

scipy.integrate: 1.83508872224
モンテカルロ積分: 1.8328253236

(6) $\displaystyle I = \int_e^{e^2} (\log x)^2 dx$

f:id:aidiary:20140728223510p:plain

a, b = np.exp(1), np.exp(2)
f = lambda x: np.log(x) ** 2
# 以下同

scipy.integrate: 12.0598303694
モンテカルロ積分: 12.0453421225

scipy.integrateに比べると少し近似精度が落ちるが積分が計算できていることがわかる。サンプル数Nを増やすと近似精度が増すこともわかる。

数値積分のアルゴリズムには、台形公式、シンプソン法、ガウス求積法などがある。scipy.integrateの実装を調べてみたが、QUADPACKに実装されているガウス求積法を使っているようだ。

実際のところ今回実験したような関数ではscipy.integrateの方が精度が高く、モンテカルロ積分を使う利点があまり感じられなかった。モンテカルロ法の利点は何なのだろう？従来のアルゴリズムの詳細を理解していないので完全に受け売りだが、モンテカルロ積分は統計や機械学習で頻繁に扱う多重積分においてより効率的とのこと。多重積分はあとで出てくるのかな？

次回はRによるモンテカルロ法入門の練習問題3.1の正規・コーシーベイズ推定量の積分を計算してみたい。複雑なのでビビるがけっこう簡単に求まる。