2014-08-20

モンテカルロ積分 (2)

機械学習 PRML

モンテカルロ積分（2014/7/28）のつづき。今回は、Rによるモンテカルロ法入門の例3.3と練習問題3.1のモンテカルロ積分の例を検証した。

例3.3

例3.3は簡単な積分の問題。

f:id:aidiary:20140821205346p:plain

これも前回と同じように確率分布 p(x) が明示的にないため一様分布を追加する。積分区間が [0, 1] なので b - a = 1 となり無視できる。

実行結果は

scipy.integrate: 0.96520093605
モンテカルロ積分: 0.966150133367

f:id:aidiary:20140820212910j:plain

練習問題3.1 正規・コーシー-ベイズ推定量の積分

練習問題3.1はもっと複雑な積分だ。どうやら正規・コーシー-ベイズ推定量と呼ばれるものらしい。何に使うのかはよくわからないけど・・・

f:id:aidiary:20140820212913j:plain

これがなぜ正規・コーシー-ベイズ推定量と呼ばれるのか謎だったのだが、上の式を下の式のように変形すると正規分布とコーシー分布が出てくることがわかった。 $\pi$ が突然現れるけど分子と分母でクリアされるので追加しても問題ない。

f:id:aidiary:20140820212923p:plain

分子と分母に出てくる積分をモンテカルロ積分で計算してみる。分子と分母はほとんど同じ形なので分子だけ説明。この積分は二通りの方法でモンテカルロ積分が適用できることがわかる。

(1) コーシー分布をpとみなしてコーシー分布からサンプリングする場合

f:id:aidiary:20140820212920p:plain

(2) 正規分布をpとみなして正規分布からサンプリングする場合

f:id:aidiary:20140820212918p:plain

初めて一様分布以外からサンプリングするモンテカルロ積分の例が出てきた。両方のパターンで積分を計算し、結果を比較してみよう。ついでにscipy.integrateの結果も出してみた。

scipy.integrate: 3.43506155523
モンテカルロ積分(1): 3.4456557745
モンテカルロ積分(2): 3.43354612482

どちらの場合でもscipy.integrateとほぼ同じ結果が得られることがわかる。何回かやるたびに結果が変わるのだけれど正規分布の方がscipy.integrateの結果に近いことが多いようだ。コーシー分布は正規分布より裾野が広いため平均から外れたとんでもないサンプルが出てしまうことがあるからかな？特にサンプリング数Nが小さいときに違いが出てくる。どちらの分布からサンプリングした方がよいかの検証は次回やってみよう。

2014-07-28

モンテカルロ積分

機械学習 PRML

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の正規・コーシーベイズ推定量の積分を計算してみたい。複雑なのでビビるがけっこう簡単に求まる。

2014-07-17

受理・棄却法 (4)

機械学習 PRML

Pythonによるモンテカルロ法入門（2014/6/20）の受理・棄却法（2014/7/12）4回目。今回は練習問題2.18を解いてみた。下のような〜分布という名前がついていない完全に任意な目標分布 $f(x)$ にしたがう乱数を受理・棄却法で生成してみよう。

$\displaystyle f(x) \propto \exp (-\frac{x^2}{2}) \{ \sin (6x)^2 + 3 \cos (x)^2 \sin (4x)^2 + 1 \}$

提案分布 $g(x)$ には標準正規分布を使う。

$\displaystyle g(x) = \frac{1}{\sqrt{2 \pi}} \exp (-\frac{x^2}{2})$

目標分布の式が示すように左辺は右辺に「比例する」だけで、積分が1になるように正規化されていない。乱数を生成するだけなら正規化は必須でないが、今回は最初から正規化しておいた。正規化するには右辺を全区間で積分した値で $f(x)$ を割ればよい。Pythonで定積分を計算するには、scipy.integrate.quadが使える。

さっそく実装してみよう。

# 任意の確率分布にしたがう乱数を受理・棄却法で生成する

import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize
import scipy.integrate
from scipy.stats import norm

np.random.seed()

# 正規化前の目標分布f
f = lambda x: np.exp(-x**2 / 2) * (np.sin(6*x)**2 + 3 * np.cos(x)**2 * np.sin(4*x)**2 + 1)

# 正規化定数
I = scipy.integrate.quad(f, -Inf, Inf)[0]

# 正規化済み（積分が1）の目標分布
f = lambda x: np.exp(-x**2 / 2) * (np.sin(6*x)**2 + 3 * np.cos(x)**2 * np.sin(4*x)**2 + 1) / I

# 提案分布g
gv = norm(loc=0.0, scale=1.0)
g = gv.pdf

# 分布の上限を指定する定数Mを設定
xopt = scipy.optimize.fmin(lambda x: - f(x) / g(x), 0.0, disp=False)[0]
# そこでの値をMとする
M = f(xopt) / g(xopt)
print "M =", M

# 受理・棄却法
Nsim = 100000

# 提案分布gからの乱数Yを生成
Y = gv.rvs(size=Nsim)

# 一様乱数UをNsim個生成
U = uniform.rvs(size=Nsim)

# Yから受理の条件を満たすサンプルXを残して残りを棄却
X = Y[U <= f(Y) / (M * g(Y))]
print u"サンプル数: %d => %d" % (len(Y), len(X))
print u"実際の受理率  : %f" % (len(X) / float(len(Y)))
print u"理論的な受理率: %f" % (1.0 / M)

# 目標分布を描画
x = np.linspace(-10, 10, 1000)
y = f(x)
plt.plot(x, y, 'r-', lw=2)

# 提案分布を描画
y = M * g(x)
plt.plot(x, y, 'g-', lw=2)

# 受理した乱数の分布を描画
plt.hist(X, bins=50, normed=True)

plt.xlim((-5, 5))
plt.show()

実行すると

M = 1.85606973547
サンプル数: 100000 => 53818
実際の受理率  : 0.538180
理論的な受理率: 0.538773

f:id:aidiary:20140717213835j:plain

赤色が目標分布 $f(x)$ 。緑色が提案分布 $g(x)$ の正規分布。青色が受理されたサンプルのヒストグラム。ヒストグラムは赤色の理論値と一致していることがわかる。

少し気になってたのだけれど、どんな分布でも理論的な受理率が $1/M$ で計算できてしまうのが不思議に感じてきた。もう少し理論的な背景も調べた方がよさそうだな。

今回で受理・棄却法はおしまい。次回から第3章のモンテカルロ積分にいく予定。乱数生成で難しい積分が効率的に計算できてしまうのだ！

人工知能に関する断創録

このブログでは人工知能のさまざまな分野について調査したことをまとめています（更新停止: 2019年12月31日）

モンテカルロ積分 (2)

例3.3

練習問題3.1 正規・コーシー-ベイズ推定量の積分

モンテカルロ積分

モンテカルロ積分

簡単な積分計算例

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

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

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

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

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

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

受理・棄却法 (4)