確率論 6

[1][2]のおかげで[3]が読めるようになってきた．

[1]結城浩, 数学ガール(乱択アルゴリズム)
[2]平岡和幸, 堀玄プログラミングのための確率統計
[3]伊庭幸人, 種村正美, 計算統計 2 マルコフ連鎖モンテカルロ法とその周辺 (統計科学のフロンティア 12)

モンテカルロ法(例題)

モンテカルロ法とは，目的の数値や関数の近似解を得るための，乱数を用いたアルゴリズムの総称．

例題．「半径1の円の面積をπを使わないで求めよ」

モンテカルロ法による解法
-1)[-1,1]×[-1,1]の二次元領域に一様分布に従う点列x0, x1,...xNを発生させる．
-2)円の面積 S = (円内に入った点数) / (総点数N) × (領域全体の面積) とする．

static double calcPI(int ptNum = 10000000)
{
    //[-1,1]×[-1,1]の範囲に一様分布に従う点をptNum個発生させる
    int numOfInsidePt = 0;
    for( int i=0; i<ptNum; ++i){
        double x = 2 * rand()/(double)RAND_MAX -1.0;
        double y = 2 * rand()/(double)RAND_MAX -1.0;
        if( x*x + y*y <= 1.0 ) ++numOfInsidePt;
    }

    //円の面積(1*1*π) = (半径1の円に入ったサンプル点数) / 全てのサンプル点数 * 全体の長方形数 
    return numOfInsidePt / (double) ptNum * (2.0*2.0);
}

実行結果
calcPI(10^2) = 3.27999999
calcPI(10^4) = 3.14639999
calcPI(10^6) = 3.14144800
calcPI(10^8) = 3.14101904
少数第4桁以上の精度は出なかった．

モンテカルロ積分

モンテカルロ積分法

-- f(x)を，d次元空間R^dに定義されたスカラー関数とする．(f(x)∈R , x∈R^d)
-- 領域ΩをR^dの部分空間とする．
一様分布に従うN個のサンプル点 x₁,...,x_Nに対して以下が成り立つ．
$\displaystyle I = \int_\Omega f({\b x})d{\b x} \simeq \frac{V}{N}\sum_{i=1}^{N} f({\b x}_i)$
ただし，Vは領域Ωの体積である．
$\displaystyle V = \int_\Omega d{\b x}$

つまり，領域Ω(d次元空間)の積分をサンプル点の和で近似できる．
サンプル点が多いほど近似制度が上がるが，計算コストもかかる．
　

この近似が成り立つ証明．

まず，Uniform Law of Large Numbersより，十分大きなNに対して，
$\displaystyle \frac{f(x_1) +... + f(x_N)}{N} \simeq E(f(X))=\int_\Omega f({\b x}) p({\b x})d{\b x}$
ただし，x1,...,x2は，領域Ω内の一様分布に従うサンプル点であり，p(x)は領域Ω内の一様分布の確率密度関数．
つまり，
$\displaystyle p({\b x}) = \frac{1}{ \int_\Omega d{\b x} } = \frac{1}{V}$
以上より，　　 $\displaystyle \int_\Omega f({\b x})d{\b x} \simeq \frac{V}{N}\sum_{i=1}^{N} f({\b x}_i)$ 　　となる．

先の円の面積の解法は，このモンテカルロ積分法を以下のセットアップで使ったものである．
領域Ω 　: [-1,1]×[-1,1]
関数f(x,y): f(x,y) = 1 if (x,y) が円の内部
　　　　　　　　 0 if (x,y) が円の外部

(本題と関係ないが，久しぶりにやったら苦戦したので，この積分問題を解いてみる．)
$\displaystyle S = \int_{-1}^1 \int_{-1}^1 f(x,y) dydx$
$\displaystyle = \int_{-1}^1 \left( \int_{-\sqrt{1-x^2}}^{+\sqrt{1-x^2}}\;\;\; 1 dy \right) dx$ 　　　　...(円の内部でf(x,y)=1, 外部で0なので)
$\displaystyle = \int_{-1}^1 2 \sqrt{1-x^2} dx$ 　　
$\displaystyle = \int_{-\pi/2}^{\pi/2} 2 \cos^2 \theta d\theta$ 　　　　...x=sinθと置くと, dx = cosθdθより
$\displaystyle = \int_{-\pi/2}^{\pi/2} 1 + \cos2\theta d\theta$ 　　...倍角の公式
$\displaystyle = \pi$

モンテカルロ法の一般形

前述のモンテカルロ積分では，一様分布に従うサンプル点群で，積分を近似していた．
より一般的に，任意の確率分布Q(X)に従うサンプル点群を用いる手法を以下に示す．

解きたい問題

■ A 離散的な確率変数を考えている場合.
確率変数Xが，正規化定数の不明な確率分布P(X=x)に従うとし，以下の期待値を求めたい．
$\displaystyle E(f(X)) = \sum_{{\bf x}\in\Omega}f({\bf x})P({\bf x})$

■ B 連続的な確率変数を考えている場合.
確率変数Xが, 正規化定数の不明な確率密度関数p_X(X=x)により定義される確率分布に従うとし，以下の期待値を求めたい．
$\displaystyle E(f(X)) = \int_{\Omega}f({\bf x})p_X({\bf x}) d{\b x}$

※xはd次元空間の位置ベクトル ※関数f(x）は何かしらの統計量を表すものとする．
※確率分布P(x)の正規化定数をZとすると， $\displaystyle P({\bf x}) = \hat{P}({\bf x})/Z$ , と書けて， $\hat{P}({\bf x})$ は計算できるものとする(Zはunknown)．
※確率密度関数p_X(x)の正規化定数をZとすると， $\displaystyle p_X({\bf x}) = \hat{p_X}({\bf x})/Z$ , と書けて， $\hat{p_X}({\bf x})$ は計算できるものとする(Zはunknown)．

解法

■ A 離散的な確率分布の場合．
1) 次の条件を満たす確率分布Q(X)を用意する．
- 確率分布Q(X)に従うサンプル点群を効率的に発生させられる．
- サンプル点xを与えると，その確率Q(X=x)が計算できる．
- 確率分布 P(X) に似ている．（似ていなくても計算はできるが効率が悪くなる）

2) 確率分布Q(X)に従うサンプル点群をN個発生させる． x₁, x₂, ..., x_N

3) 期待値と正規化定数を以下の通り計算する．
$\displaystyle w_i = \frac{ \hat{P}({\b x}_i)}{ Q({\b x}_i)}$
$\displaystyle Z \simeq \frac{1}{N} \sum_{k=1}^N w_i$
$\displaystyle E(f(X)) = \sum_{{\bf x}\in\Omega}f({\bf x})P({\bf x}) \simeq \frac{1}{NZ} \sum_{k=1}^N f({\b x}_i) w_i$

■ B 連続的な確率分布の場合．--- (確率分布が確率密度関数になるだけ)
1) 次の条件を満たす確率分布Q(X)を用意する．
- 確率分布Q(X)に従うサンプル点群を効率的に発生させられる．
- サンプル点xを与えると，その確率密度q(x)が計算できる．
- 確率密度関数が求めたいp_X(x)に似ている．（似ていなくても計算はできるが効率が悪くなる）

2) 確率分布Q(X)からサンプル点群をN個発生させる．　x₁, x₂, ..., x_N

3) 期待値と正規化定数を以下の通り計算する．
$\displaystyle w_i = \frac{ \hat{p_X}({\b x}_i)}{ q({\b x}_i)}$
$\displaystyle Z \simeq \frac{1}{N} \sum_{k=1}^N w_i$
$\displaystyle E(f(X)) = \int_{\Omega}f({\bf x})p_X({\bf x}) d{\b x} \simeq \frac{1}{NZ} \sum_{k=1}^N f({\b x}_i) w_i$

証明) 離散の場合
正規化定数の定義　 $\displaystyle P({\bf x}) = \hat{P}({\bf x})/Z$ 　より　 $\displaystyle w_i = \frac{ P({\b x}_i) Z}{ Q({\b x}_i)}$ ．

$\displaystyle \frac{1}{N}\sum_{i=1}^N w_i = \frac{1}{N}\sum_{i=1}^N \frac{ P({\b x}_i) Z}{ Q({\b x}_i)} \simeq \sum_{{\b x}\in \Omega} \left{ \frac{ P({\b x}) Z}{ Q({\b x}) } Q({\b x}) \right}$ 　 $\displaystyle = Z \sum_{{\b x}\in \Omega} P({\b x}) =Z$

Uniform Law of Large Numbersより， $\frac{ \sum_{i=1}^N f({\b x}_i) }{N} \simeq E(f(X) )$ とできる. この変型では $f({\b x})=P({\b x})Z/Q({\b x})$ 　と置いている．

また，期待値に関して，
$\displaystyle \frac{1}{NZ} \sum_{i=1}^N f({\b x_i}) w_i = \frac{1}{NZ} \sum_{i=1}^N f({\b x_i}) \frac{ P({\b x}_i) Z}{ Q({\b x}_i)}$ $\displaystyle \simeq \sum_{{\b x}\in \Omega} \left{ \frac{ f({\b x}) P({\b x})}{ Q({\b x}) } Q({\b x}) \right} = \sum_{{\b x}\in \Omega} f({\b x}) P({\b x})$

連続の場合は，x∈Ω の和分がΩ内での積分になり，確率分布P(x)/Q(x)が確率密度関数p(x)/q(x)に置き換わるだけ．

モンテカルロ法の正しさは，大数の法則(Uniform Law of Large Numbers)に基づいている．
Nが十分大きければ (x1+x2+...+xN)/N ≒ ΣxP(x) とNでの割り算が確率分布Pの掛け算になっていることに注意．

離散の例．イジングモデル

連続の例．上の面積積分が連続確率分布の良い例．
上の例では，p(x)を2次元における一様分布と定義し (p(X) = 1/4，x∊[-1,1]x[-1,1])，
f(x) = 1 if　点xが円の内部
　　= 0 others
とした時の期待値 E(f(X)) を，モンテカルロ法で求めている．

確率論7へ