Chapter 04

確率・統計入門

mebiusbox
mebiusbox
2021.02.24に更新

📌 はじめに

コンピュータグラフィックスでは,様々なところで確率と統計が使用されています.
例えば,レイトレーシングによるレンダリングでは,モンテカルロ法を使って光輸送のシミュレーションを行い,計算を高速に行えます.モンテカルロ法は乱数を使って数値積分を行う方法です.また,物理ベースレンダリングでよく使用されるマイクロファセット反射モデルでは,マイクロファセットの分布に確率が使用されています.他にもシャドウマップの影生成手法には確率を使ってソフトシャドウをシミュレーションするバリアンスシャドウマッピングがあります.このように確率と統計はコンピュータグラフィックスにおいて重要な学問ではあるのですが,扱っている範囲が広く,個人的には難しい印象でしたし,今でもそう感じています.確率・統計は用語が多く,しかも似たような名前も少なくないので,混乱しがちです.そのため,用語は出来るだけ絞っています.目標はモンテカルロ法によって数値積分できるようになるところです.

📌 数式処理ソフトウェア

数式の計算や,シミュレーションのために Maxima と Scilab というソフトウェアを利用しています.本文中に各ソフトウェアのコマンドを記載しているところがあります.実行するにはインストールする必要があります.下記からダウンロードできます.

📌 確率

まずは高校で学ぶ確率についてです.これは数学者ラプラスが定義した確率論に基づいています.
確率の例として必ず出てくるサイコロを使いましょう.よくある 6 面で 1 から 6 の目があるものです.24 面とか 32 面とかではありませんよ.あとサイコロは 1 個です.

試行,標本,事象

今,サイコロを 1 回振りました.出た目は「1」でした.もう一度振ってみると「2」でした.細工をしない限り1回目と2回目の出る目は偶然によって決まります.このように同じ条件で何度でもサイコロを振ることを「試行」といいます.試行によって起こりえる結果はサイコロの目である「1, 2, 3, 4, 5, 6」のことで,1つひとつを「標本」といいます.そして,すべての標本(ここでは 1, 2, 3, 4, 5, 6 の6つ)を「全事象」といいます.さらに全事象から取り出した標本の集まりを「事象」といいます.事象は1つの標本でもいいし,全部(すなわち全事象)でも構いません.

「全事象」は「標本空間」ともいいます.また,標本が1つだけの事象を「根元事象」といいます.

集合,要素

互いに識別できる「もの」の集まりを,「集合」といいます.たとえば,1から5までの自然数 1, 2, 3, 4, 5 の集まりは集合です.この場合,\{1,2,3,4,5\},あるいは A=\{1,2,3,4,5\} と表します.1つひとつを「要素」といいます.要素をまったくもたない集合のことを「空集合」といい,\phi で表されます.ある集合 A の要素の個数は n(A) と表します.

共通部分と和集合

集合 B の要素がすべて集合 A の要素であるとき,BA の「部分集合」であるといいます.そして B \subset A, または A \supset B と表します.

一般に,集合 U を全体の集合とし,この中で部分集合を考えるとき,U を「全体集合」といいます.部分集合 AB の要素を合わせた集合を「和集合」といい,A \cup B と表します.また,部分集合 AB のどちらにも入っている共通の要素を「共通要素」といい,A \cap B と表します.全体集合 U のうち,部分集合 A ではない要素を A の「補集合」といい,\overline{A} と表します.

全体集合 U,部分集合 A,補集合 \overline{A} は次のような関係になります.

A \cup \overline{A} = U, \quad A \cap \overline{A} = \phi

確率の定義

数学的確率

標本空間 U が,n 個の根元事象で構成されていて,どの2つの根元事象も”重複して起こらず”,また,どの根元事象の起こることも,”同様に確からしい”とします.ある事象 Aa 個の根元事象から成り立つとき,つまり,A の起こる場合が a 通りであるとき,

p = \frac{a}{n},

のことを,事象 A の起こる確率といい,P(A) と表します.これを数学的確率と呼びます.

確率 P(A)=0 のとき,事象 A は絶対に起こりません.P(A)=1 のときは,必ず起こります.サイコロの例で1の目が出る確率は,事象 A=\{1\} で,全事象 U=\{1,2,3,4,5,6\} なので,確率は次のようになります.

P(A) = \frac{1}{6}.

統計的確率

ある試行の全事象を \Omega と表します.この試行の1つの事象を A とします.
全事象の個数を n(\Omega),事象 A の個数を n(A) とすると,事象 A が起こる確率 P(A) は次のようになります.

P(A) = \frac{n(A)}{n(\Omega)}

この確率は統計的確率といいます.数学的確率は同様に確からしいという状況のもとに考えるのに対して,統計的確率では実際に試行した結果をもとに考えます.たとえば,サイコロに細工をして1の目が他の目よりも出やすくなっているものは同様に確からしいとは言えません.この場合,1の目が出る確率は \frac{1}{6} ではないので,実際に試行した結果をもとに確率を算出します.

ある試行を n 回観測します.このとき事象 A が起きた回数を r としたとき,

\frac{r}{n},

これを「相対頻度」といいます.相対頻度が,n を限りなく大きくしていったときに,一定の値 \alpha に限りなく近づくとすれば,事象 A の確率 P(A)

P(A) = \alpha

と定められます.大数の法則(Law of Large Numbers,またはベルヌーイの定理)によって,このことが証明されています.これはモンテカルロ法の根源となる定理です.

排反

事象 A と事象 B が同時に起こらないとき,AB は「排反である」といいます.たとえばサイコロの場合,事象 A が偶数 \{2,4,6\} で,事象 B が奇数 \{1,3,5\} だと排反であるといえます.しかし,事象 A が奇数 \{1,3,5\} で,事象 B が3以下 \{1,2,3\} の場合は「排反ではない」ということになります.

加法定理

事象 A と事象 B が排反であるならば,1回の試行で事象 A と事象 B が起こる確率は以下のようになります.

P(A \cup B) = P(A) + P(B)

条件付き確率

事象 B が事象 A が起きたことで確率が変わる場合,事象 B のことを「条件付き確率」といいます.この条件付き確率は P_A(B) または P(B|A) と表します.

例えば,10本のくじがあり,3本の当たりくじがあったとします.先に太郎が1本引き,次に次郎が1本引くとします.太郎が当たりくじを引く確率は10本中の3本なので,\frac{3}{10}です.次郎は,もし太郎が当たりくじを引いていたら残り9本中の2本となるので,当たりくじを引く確率は \frac{2}{9} です.太郎が当たりくじを引いていなかったら9本中の3本となるので,当たりくじを引く確率は \frac{3}{9} となります.つまり,次郎の当たりくじを引く確率は条件付き確率ということになります.

乗法定理

条件付き確率のときの当たりくじにおいて,太郎と次郎が二人とも当たりくじを引く確率は次のようになります.

P(A \cap B) = P(A) P_A(B)

📌 統計

統計には,一部のデータから全体を推測する「推測統計」と,全体のデータから役に立つ情報を得る「記述統計」があります.モンテカルロ法は確率によって真値を推測するため,推測統計に分類されます.

📌 データ,階級,度数,度数分布

確率でいう試行を統計では観測(実験や調査)といいます.観測で得られた観測値を「データ」といいます.データをいくつかの「階級」に分けて考えたとき,各階級ごとのデータの数を「度数」といいます.この度数を縦軸に,階級を横軸にして表示したグラフを「ヒストグラム」といいます.ヒストグラムによってデータがどのように分布されているのかが分かりやすくなります.

📌 確率変数

試行の結果によって値が変わるものを「確率変数」といい,X, Y などと大文字を用いて表されます.ある試行の結果 \omega によって X の値が決まる場合,X(\omega) と書くことができます.確率変数には数値が入ります.たとえばコインの場合,標本空間は表と裏ですが,表が 1,裏が 0 のように具体的な数値が入ります.

📌 確率分布

たとえば,サイコロを振ったときの出る目を確率変数とすると,確率変数に入る値は \{1,2,3,4,5,6\} のどれかになります.それぞれ出る頻度は確率で表されます.ある確率変数 X に対して,値とその確率をまとめたものを「確率分布」といい,それを表にしたものを「確率分布表」といいます.例えば,サイコロの場合,

Xの値 1 2 3 4 5 6
確率 1/6 1/6 1/6 1/6 1/6 1/6

サイコロに細工をして 1 の目だけ出るようにした場合は次のようになります.

Xの値 1 2 3 4 5 6
確率 1 0 0 0 0 0

確率分布の特徴は「期待値(平均)」と「分散」でわかります.

📌 期待値(平均値)

特徴を表す代表値に「期待値」というのがあります.ある確率変数 X の期待値は E[X] と書きます.

確率変数の取る値を x_1,x_2,\ldots,x_n とし,x_1 が出る確率を p_1 だとします.表にすると次のようになります.

Xの値 x_1 x_2 \ldots x_n
確率 p_1 p_2 \ldots p_n

期待値を数式で表すと

E[X] = \sum_{i=1}^{n}x_i\cdot p_i.

この期待値は平均とも考えられます.

確率変数 X の値に対応した確率の合計は 1 になります.

\sum_{i=1}^n p_i = 1.

📌 分散,偏差,標準偏差

確率変数 X と期待値 E[X] との差を「偏差」といいます.

|x_1 - E[X]|, |x_2 - E[X]|, ..., |x_n - E[X]|.

絶対値を含むと取り扱いがしづらいので2乗します.2乗したものの総和を「分散」といいます.確率変数 X の分散は V[X] と書きます.期待値を \mu = E[X] とおくと

V[X] = \sum_{i=1}^n (x_i - \mu)^2 \cdot p_i.

期待値を使って表すと

V[X] = E[(X-E[X])^2].

次のように変形することができます.

V[X] = E[X^2-E[X]^2].

分散は2乗しているため,平方根を取ることで確率変数 X の値との単位を一致させます.これを「標準偏差」といい,S[X] で表します.

S[X] = \sqrt{V[X]}

📌 確率密度

これまでの確率は離散的な標本を対象にしていました.今度は連続的な標本を対象にします.たとえば,高校のあるクラスの身長を考えてみます.160cm や 160.5cm とか 175.4cm とか様々です.あるクラスの人数が30人ならば,最大で30標本あることになりますが,標本としては任意の値が取れることになります.この場合,例えば 160cm~165cm の身長の人がどれくらいの確率なのかというのがわかるとよさそうです.このように連続的な値に対する確率を表すときは「確率密度関数(Probability Density Function: pdf)」を使います.ある区間[a,b]の確率は,

P(a \leq X \leq b) = \int_{a}^{b}f(x)dx.

ここで f(x) が確率密度関数です.

確率密度関数はすべての領域で積分したときに1になっていなければなりません.

\int_{-\infty}^{\infty}f(x)dx = 1.

ある値 x までの確率密度を知りたい場合は -\infty から x までの領域で定積分します.これを「累積分布関数(Cumulative Distribution Function: cdf)」といいます.

cdf(x) = \int_{-\infty}^{x}f(x)dx = P(X \leq x).

累積分布関数を微分すると,確率密度関数になります.

f(x) = \frac{d}{dx}cdf(x).

連続的な確率変数 X の期待値(E[X]=\mu),分散(V[X]=\sigma^2)は次のように積分で求めます.

\begin{aligned} \mu &= \int_{a}^{b}x f(x)dx, \\[2ex] \sigma^2 &= \int_{a}^{b}(x-\mu)^2f(x)dx. \end{aligned}

標準偏差(S[X]=\sigma)は次のようになります.

\sigma = \sqrt{\sigma^2}.

📌 様々な確率分布

ベルヌーイ分布

Bernouilli distribution.

コイン投げを考えます.標本空間は {表,裏} でそれぞれ {1,0} ,表の確率を p,裏の確率を 1-p とします.このように起こり得る結果が2通りで確率が等しくない試行を表す確率分布が「ベルヌーイ分布」です.このときの試行を「ベルヌーイ試行」ともいいます.

ベルヌーイ分布の期待値は p で,分散は p(1-p) です.

二項分布

Binomial distribution.

ベルヌーイ試行を n 回行い,それが成功する確率を p としたときの確率分布を「二項分布」といいます.Bi(n,p) と表すことがあります.ベルヌーイ分布は Bi(1,p) のことです.

二項分布の期待値は np で,分散は np(1-p) です.

一様分布

Uniform distribution.

確率変数 X のとる値がすべて同じ確率のとき,この確率分布を「一様分布」といいます.連続的な確率変数の場合,確率密度関数が f(x)=1 であることを意味しています.この場合,期待値と分散は次のようになります.

\begin{aligned} \mu &= \int_{a}^{b}x f(x)dx = \int_{a}^{b}xdx = \frac{1}{2}, \\[2ex] \sigma^2 &= \int_{a}^{b}\left(x-\frac{1}{2}\right)^2 f(x)dx = \int_{a}^{b}\left(x-\frac{1}{2}\right)^2 dx = \frac{1}{12}. \\[2ex] \end{aligned}

正規分布

Normal distribution.

ガウス分布(Gaussian distribution)ともいわれています.下図のような形をしています.

tuto-probability-normal-distribution

この曲線をベル曲線といいます.この図は maxima でプロットしたもので,次のようなコマンドです.

plot2d(%e^(-x^2),[x,-3,3],[ylabel,"y"],[title,"Normal distribution"]);

ベル曲線を数式で表すと

e^{-x^2}.

正規分布の確率密度関数は

f(x) = \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(x-\mu)^2}{2\sigma^2}}.

\frac{1}{\sqrt{2\pi}\sigma} は以下を満たすための正規化定数です.

\int_{-\infty}^{\infty} f(x)dx = 1

これは

\int_{-\infty}^{\infty} e^{-\frac{(x-\mu)^2}{2\sigma^2}} = \sqrt{2\pi}\sigma,

から来ています.
正規分布の期待値は

E[X] = \int_{-\infty}^{\infty} x \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(x-\mu)^2}{2\sigma^2}}dx = \mu.

分散は

V[X] = \int_{-\infty}^{\infty} (x-\mu)^2 \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(x-\mu)^2}{2\sigma^2}}dx = \sigma^2.

つまり,期待値と分散は指定した \mu\sigma^2 そのものとなります.正規分布を,期待値(平均)\mu,分散 \sigma^2 から N(\mu,\sigma^2) と表すことがあります.

その他の分布

指数分布,ポアソン分布,ガンマ分布,ベータ分布,\chi^2(カイ二乗)分布,コーシー分布,対数分布,パレート分布,ワイブル分布,t分布,etc.

📌 モンテカルロ法

モンテカルロ法は乱数を使って数値計算を行う方法です.コンピュータグラフィックではレイトレーシングでよく使われています.その特徴に,高次元積分でも適用できる,計算が早いといったメリットがありますが,真値への収束が遅いといった欠点があります.

円周率を求めてみる

よくある例として円周率 \pi を求めてみます.
半径が 1 の円を考えます.この円の面積は S = r^2\pi です.また,円を囲む正方形は一辺の長さが 2r となります.この正方形内に無作為に N 個の点を描いたとき,点が円の内部に入っている割合は,円の面積と正方形の面積の比になると推定されます.大数の法則によれば,N を十分に大きくすると,実際の比と考えることができます.

\frac{r^2\pi}{(2r)^2} = \frac{\pi}{4}

試しにやってみます.以下は Scilab のコマンドです.

n = 10;
x = rand(1,n)*2-1;
y = rand(1,n)*2-1;
r = x.*x + y.*y;
ans=length(find(r<=1))/n*4
disp(ans);

n を変えて計算した結果は次の通りです.

n answer
10 4
1000 3.128
1000000 3.14118

n の回数を増やせば \pi に近づいていることが分かります.また,乱数のため実行する度に結果が変わります.

モンテカルロ法による数値積分の基礎

例えば1次元の積分で考えてみます.

I = \int_{a}^{b}f(x)dx.

この値を求めるために,区間 [a,b] 上での一様乱数 v を発生させます.このような乱数は区間[0,1]の一様乱数 u から v = a+u(b-a) で求まります.次に

X = (b-a)f(v),

を計算してみます.確率変数 X の期待値 E[X]v の確率密度関数が

p(v) = \begin{cases} \frac{1}{b-a} && (v\in[a,b]) \\ 0 && (v\notin[a,b]) \end{cases},

なので,

E[X] = \int_{a}^{b}X p(v)dv = \int_{a}^{b}(b-a)f(v)\frac{1}{b-a}dv = \int_{a}^{b}f(x)dx = I,

となって,求める積分の値に等しくなります.

乱数 v を独立に N 個発生させて,それぞれに対して X_i = (b-a)f(v_i) を計算して,それらの N 個の X_i の平均をとると次のようになります.

\hat{X}_N = \frac{1}{N}\sum_{i=1}^N X_i = \frac{b-a}{N}\sum_{i=1}^N f(v_i).

その期待値と分散は次のようになります.

\begin{aligned} E[\hat{X}_N] &= \frac{b-a}{N}\sum_{i=1}^N f(v_i) \\[2ex] &= (b-a)\cdot E[f(v_i)] \\[2ex] &= (b-a)\cdot \int_{a}^{b} f(v)p(v)dv \\[2ex] &= (b-a)\cdot \int_{a}^{b} f(v)\frac{1}{b-a}dv \\[2ex] &= \int_{a}^{b} f(v)dv = I. \\[2ex] V[\hat{X}_N] &= V\left[\frac{1}{N}\sum_{i=1}^{N}X_i\right] \\[2ex] &= \frac{1}{N^2}V\left[\sum_{i=1}^{N}X_i\right] \\[2ex] &= \frac{1}{N^2}\sum_{i=1}^{N}V[X_i] \\[2ex] &= \frac{1}{N}\sigma^2. \end{aligned}

まとめると

E[\hat{X}_N] = I, \quad V[\hat{X}_N] = \frac{1}{N}\sigma^2.

標準偏差 \sigma[\hat{X}_N] は次のようになります.

\sigma[\hat{X}_N] = \frac{1}{\sqrt{N}}\sigma[f(X)]

実際にモンテカルロ法で積分してみる

次の積分を求めてみます.

I = \int_{0}^{2}x^2dx

モンテカルロ法による一様乱数 u での推定値は次のようになります.

\hat{I} = \frac{2}{N} \sum_{i=1}^{N} f(X_i)

確率変数 X は次のようになります.

X = 2u

これでシミュレーションしてみます.以下は scilab のコードです.

N=[10, 1000, 1000000];
for i=1:length(N)
    x=2*rand(1,N(i));
    ans=2*mean(x.*x);
    disp(ans);
end

結果は次のとおりです.乱数を使っているため,実行するたびに答えが変わるので注意してください.

n answer
10 3.0309054
1000 2.7869074
1000000 2.6679645

今度は maxima を使って積分を計算してみます.

integrate(x^2,x,0,2);

\frac{8}{3} と表示されました.float(8/3) と入力すると 2.666666666666667 でしたのでほぼ真値になっています.

分散を減らす方法

モンテカルロ法の弱点は収束が遅いことです.そのため,サンプリング数を大きくしないと真値に近づきません.収束を早くする,つまり分散値を0に近づけることが重要になります.

層化サンプリング

Stratified sampling.

積分区間を小区間に分けて,小区間内で無作為に1つのサンプルを抽出します.

円周率を層化サンプリングで求めてみます.ここでは正方形を n\times n 個に分割します.

以下は scilab のコマンドです.

n=10;
x=rand(n,n);
y=rand(n,n);
r=zeros(n,n);
for i=1:n
    for j=1:n
        X=2*((i+x(i,j))/n)-1;
        Y=2*((j+y(i,j))/n)-1;
        r(i,j) = X*X+Y*Y;
    end
end
ans=4*length(find(z<=1))/(n*n);
disp(ans);

実行結果は次の通りです.

n answer
10 2.76
100 3.1264
1000 3.141336

通常のモンテカルロ法よりも少ないサンプルで真値に近づいていることがわかります.
実装が比較的簡単ではありますが,次元が増えるとサンプル数も指数的に増加していきますので,高次元では使用しにくい問題があります.

重点的サンプリング

Importance sampling.

被積分関数 f(x) の値の大きさに比例するような確率分布の乱数を使うことで分散を小さくすることができます.

一様分布ではない確率密度関数 p(x) の確率変数を X とすると

I = \int_{a}^{b}\frac{f(x)}{p(x)}dx.

この推定値 \hat{X}_N は次のようになります.

\hat{X}_N = \frac{1}{N}\sum_{i=1}^{N}\frac{f(X_i)}{p(X_i)}.

次に期待値 E[\hat{X}_N] は次のようになります.

\begin{aligned} E[\hat{X}_N] &= E\left[\frac{1}{N}\sum_{i=1}^{N}\frac{f(X_i)}{p(X_i)}\right], \\[2ex] &= \frac{1}{N}\sum_{i=1}^{N}E\left[\frac{f(X_i)}{p(X_i)}\right], \\[2ex] &= \frac{1}{N}\sum_{i=1}^{N}\int_{a}^{b}\frac{f(x)}{p(x)}p(x)dx, \\[2ex] &= \frac{1}{N}\sum_{i=1}^{N}\int_{a}^{b}f(x)dx, \\[2ex] &= \int_{a}^{b}f(x)dx = I \end{aligned}

分散 V[\hat{X}_N] は次のようになります.

V[\hat{X}_N] = \frac{1}{\sqrt{N}}\sigma\left[\frac{f(X)}{p(X)}\right].

次の積分を重点的サンプリングで求めてみます.

\int_{0}^{2}x^2dx

区間[0,2]の確率密度関数 p(x) を考えてみます.p(x)x に比例するとします.数式で表すと

P(0 \leq x \leq 2) = \int_{0}^{2}Cp(x)dx.

C は定数で傾きになります.
確率密度関数は積分すると 1 になることがわかっていますので

\int_{0}^{2}p(x)dx = 1,

を満たさなければなりません.これを満たす定数 C を求めます.p=Cx とすると

\int_{0}^{2}Cxdx = \left[\frac{1}{2}Cx^2\right]_{0}^{2} = 2C = 1,

よって

C = \frac{1}{2},

だから

p(x) = \frac{x}{2}.

被積分関数は x^2 で,確率密度関数は \frac{x}{2},区間[0,2] を図にしてみると次のようになります.

tuto-probability-pdf

確率密度が一様分布(緑)より,求めた確率密度関数(青)の方が被積分関数(赤)に形状が似ています.このように被積分関数の形状に近い確率密度関数を決めることが重点的サンプリングの重要なことです.

確率密度関数 p(x) の累積分布関数 cdf(x)

cdf(x) = \int_{0}^{x}\frac{x}{2}dx = \frac{1}{4}x^2 \quad (0 \leq x \leq 2)

となります.例えば cdf(1.0) = \frac{1}{4} となり,これは x が 1.0 のとき,確率が 25% になるということを表しています.今,一様乱数 u を使って確率密度が \frac{1}{4}x^2 となる関数 g(u) を求めたいわけです.これは g(0.25) = 1.0 となるはずです.

g(cdf(u)) = u

なので

g(u) = cdf^{-1}(u)

cdf(x) = \frac{1}{4}x^2 なので,

\begin{aligned} u &= \frac{1}{4}x^2 \\[2ex] 4u &= x^2 \\[2ex] x &= \sqrt{4u} \end{aligned}

このようにある分布の確率密度関数を逆関数を使って求める方法を逆関数法といいます.

求めた確率密度関数を使ってシミュレーションしてみます.以下は scilab のコードです.

function ret = pdf(x)
    ret = 0.5*x;
endfunction

N=[10, 1000, 1000000];
for i=1:length(N)
    x=sqrt(4*rand(1,N(i)));
    y=(x.*x)./pdf(x);
    ans=mean(y);
    disp(ans);
end

結果は次の通りです.

n answer
10 3.0600588
1000 2.6592209
1000000 2.6657773

重点的サンプリングにおいて,理想的な確率密度関数は p(x) = \frac{f(x)}{I} です.すでに I8/3f(x) = x^2 と分かっているので

p(x) = \frac{3}{8}x^2.

逆関数法を使って,先ほどの確率密度に従う確率密度関数を計算します.まず,累積分布関数を求めると

\begin{aligned} cdf(x) &= \int_{0}^{x} p(x) \\[2ex] &= \int_{0}^{x}\frac{3}{8}x^2 \\[2ex] &= \left[\frac{1}{8}x^3\right]_{0}^{x} \\[2ex] &= \frac{1}{8}x^3. \end{aligned}

次に逆関数を求めます.

cdf^{-1}(x) = 8x^{\frac{1}{3}}.

分散が0なので,サンプル数は1つで十分です.実際にシミューレーションしてみます.以下は scilab のコードです.

function ret = pdf(x)
    ret = 3*x*x/8;
endfunction

N=[1];
for i=1:length(N)
    x=(8*rand(1,N(i)))^(1/3);
    y=(x.*x)./pdf(x);
    ans=mean(y);
    disp(ans);
end

結果は次のとおりです.

n answer
1 2.6666667

📌 乱数の生成

重点的サンプリングするために,ある確率分布に従った乱数を作成する必要があります.既に,逆関数法を使ったやり方を説明しました.しかし,必ずしも逆関数法が使えるわけではありません.条件に累積分布関数が単調増加(または単調減少)である必要があります.

正規分布乱数

正規分布に従う乱数を作成する場合,その累積分布関数の逆関数を求めるのは簡単ではないので,中心極限定理を使った方法とボックス・ミュラー法を使って求めるやり方があります.

大数の法則

Law of large numbers.

確率のところで少し出てきましたが,中心極限定理の前に大数の法則について,あらためて見てみます.

コイン投げを考えます.標本空間は {表,裏} で {1,0} とし,10回の試行をしました(2つのうちのいずれかの事象が発生する試行をベルヌーイ試行といいました).ここで,確率変数 X_i は,i 回目のコイン投げの結果です.表の出た回数(頻度)は,確率変数の和となります.

r = X_1 + X_2 + \ldots + X_{10}

表の出た回数の割合 \hat{p} = \frac{r}{10} は観測された成功率で,\hat{p} = 0, 0.1, 0.2 \cdots となります.一般に,試行回数を n とするとき,\frac{r}{n} を相対頻度といいます.r は確率変数で n=10p=0.5 の二項分布 Bi(10,0.5) に従います.期待値と分散は二項分布の計算式を使うとそれぞれ

\begin{aligned} E[r] &= np=5 \\[2ex] V[r] &= np(1-p)=2.5, \end{aligned}

となります.相対頻度の期待値と分散は次のとおりです.

\begin{aligned} E\left[\frac{r}{n}\right] &= p=0.5 \\[2ex] V\left[\frac{r}{n}\right] &= \frac{p(1-p)}{n}=0.025. \end{aligned}

ここで,p=0.5 は真の成功率です.観測された成功率とその確率を表にまとめてみます.二項分布の確率は f(x) = {}_{n} C _{x} p^x(1-p)^{n-x} です({} _{n} C _{r} は組み合わせで \frac{{} _n P _r}{r!}{} _n P _r は順列で \frac{n!}{(n-r)!} です).今回の場合は f _i(x) = {} _{10} C _x (\frac{1}{2})^{10} となります.

観測された成功率 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
確率 0.001 0.010 0.044 0.117 0.205 0.246 0.205 0.117 0.044 0.010 0.001

下記は scilab で確認したコードです.

function y = permutations(n,r)
    y=factorial(n)/factorial(n-r);
endfunction

function y = combinations(n,r)
    y = permutations(n,r)/factorial(r);
endfunction

n=10;
x=zeros(10);
for i=0:n
    x(i+1) = combinations(10,i)*(1/2)^(10);
end
disp(x);

この表を見てもわかりづらいので,ヒストグラムにしたものが下図です.

tuto-probability-large-numbers

なんだかそれっぽい図ですね.

(ヒストグラムは先程の scilab のコードの最後に bar(0:0.1:1,x,0.2); を追加します)

ここで,コイン投げの回数 n を増やしていくと \frac{r}{n} の期待値 E[\frac{r}{n}] = p,およびその周辺が発生する確率が上がっていき,最終的に真の成功率 p=0.5 に集中し,相対頻度は真の成功率に収束していきます.これが大数の法則です.

Bi(1,0.5) のベルヌーイ試行を繰り返していき,相対頻度の変化をグラフにしたものが下図です.真の成功率 p=0.5 近くに収束していくのがわかります.

tuto-probability-bernoilli

以下は scilab のコードです.

n=10000;
x=bool2s(rand(1,n)<0.5);
y=cumsum(x)./(1:n);
plot(y);

a=gca();
a.data_bounds(:,2)=[0.4,0.6];

中心極限定理

Central limit theorem.

期待値 \mu,分散 \sigma^2 である確率変数 X を考えます.この確率分布は一様分布かもしれないし,二項分布かもしれませんし,なんでも構いません.この試行を何度も繰り返していくと,標本平均(確率変数の平均)の確率分布はおおむね正規分布になるというのが中心極限定理です.

標本平均 S_n とすると

S_n = \frac{X_1 + X_2 + \ldots + X_n}{n}.

n が大きくなるにつれて S_n の分散は近似的に N(\mu,\frac{\sigma^2}{n}) となることがわかっています.また,

Z_n = \frac{S_n-\mu}{\sigma/\sqrt{n}},

と定義し,S_n を標準化すると n が大きくなるにつれて,確率変数 Z_n は標準正規分布 N(0,1) に近づいていきます.整理すると

Z = \frac{S-\mu}{\sigma},

この Z は標準化変数といいます.

二項分布による正規分布の近似

中心極限定理によれば,n が大きいときに Bi(n,p) は正規分布に近づくので,それを利用して正規分布乱数を作成することができます.

二項分布の期待値は \mu = np,分散は \sigma^2 = np(1-p) なので,標準化
変数に代入すると

Z = \frac{S-np}{\sqrt{np(1-p)}}.

Bi(1000,0.5) を scilab でシミューレーションしたコードと結果は次のとおりです.

p=0.5;
n=1000;
N=10000;
x=rand(N,n);
y=bool2s(x<p);
S=sum(y,"c");
m=n*p;
sigma=sqrt(n*p*(1-p));
Z=(S-m)/sigma;
histplot(50,Z,style=12);

tuto-probability-normal-random-binomial

正規分布乱数の生成

一様乱数から正規分布乱数を生成します.これも中心極限定理に基づきます.
一様乱数の期待値と分散は \mu=\frac{1}{2}\sigma^2=\frac{1}{12} です.標準化変数は

z = \left(\sum_{i=1}^{n}X_i-\frac{n}{2}\right) \frac{1}{ \sqrt{n/12}}.

期待する正規分布 x の期待値 \mu,分散を \sigma^2 とすると,x の標準化の逆変換 x = \sigma z + \mu に代入して

x = \sigma\sqrt{\frac{12}{n}} \left( \sum_{i=1}^{n}X_i-\frac{n}{2} \right) + \mu.

これが,正規分布 N(\mu,\sigma^2) に従う正規分布乱数です.

N(10,5) の正規分布を scilab でシミュレーションしたコードとヒストグラムは次のようになっています.

sigma=5;
mu=10;
n=1000;
N=1000;
x=rand(N,n);
S=sum(x,"c");
Z=sqrt(sigma)*sqrt(12/n)*(S-n/2)+mu;
histplot(50,Z);

tuto-probability-normal-random-uniform

標準化変数で n=12 とすると式を簡単にすることができます.

x = \sigma \left(\sum_{i=1}^{12}X_i-6\right)+\mu.

これもシミュレーションした結果は次のとおりです.

sigma=5;
mu=10;
N=1000;
x=rand(N,12);
S=sum(x,"c");
Z=sqrt(sigma)*(S-6)+mu;
histplot(50,Z);

tuto-probability-normal-random-uniform-opt

ボックス=ミュラー法

Box–Muller's method.

詳しい証明はここでは行いませんが,以下の式を使うと正規分布乱数を生成できます.Z_1, Z_2 どちらかを使います.

\begin{aligned} Z_1 &= \sqrt{-2\log X}\cos(2\pi Y) \\[2ex] Z_2 &= \sqrt{-2\log X}\sin(2\pi Y). \end{aligned}

分散と期待値を指定する場合は

\begin{aligned} Z_1 &= \sigma\sqrt{-2\log X}\cos(2\pi Y)+\mu \\[2ex] Z_2 &= \sigma\sqrt{-2\log X}\sin(2\pi Y)+\mu. \end{aligned}

N(10,5) の正規分布を scilab でシミュレーションしたコードとヒストグラムは次のようになっています.

sigma=5;
mu=10
N=1000;
x=rand(1,N);
y=rand(1,N);
z=sigma*sqrt(-2*log(x)).*cos(2*%pi*y)+mu;
histplot(50,z);

tuto-probability-normal-random-box-muller

📌 最後に

モンテカルロ法による数値計算に必要と思われる確率と統計に関してまとめてみました.参考にした書籍は入門,超入門とかなので,まだまだ入り口部分の内容かもしれません.私自身は専門ではないので,間違いなど見つけましたらご連絡いただけると幸いです.
確率密度関連では図ももう少し入れらたらわかりやすいので,いずれ入れておこうと思います.

📌 参考文献

  • 郡山彬,和泉澤正隆,「確率と統計超入門」, 2001
  • 東京大学教養学部統計学教室,「統計学入門」,2017
  • リチャード・エルウィス著,宮本寿代訳,「マスペディア1000」,2016
  • 涌井良幸,「『数学』の公式・定理・決まりごとがまとめてわかる事典」,2016
  • 伊里正夫,藤野和建,「数値計算の常識」, 2010
  • Shocker,「モンテカルロ法 (Monte Carlo method)」Link
  • Shocker,「モンテカルロ法の分散低減手法」, Link
  • Ushio,「確率密度関数からモンテカルロ積分まで」,Link