まえがき

本記事は「日本語プログラミング言語「プロデル」 Advent Calendar 2021」18日目の記事です．
17日目は @jpl_produire さん，19日目は @sanei_media さんです．

プロデルは日本語プログラミング言語の中でも特に自然な日本語に近い表現ができる言語です．
本記事ではプロデルで乱数サンプリングアルゴリズムを実装し，コードからアルゴリズムを伝えたいという試みです．

※コンピューター上で（擬似）乱数を得るためのアルゴリズム（線形合同法やメルセンヌ・ツイスタなど）の内容については扱わず，組み込み機能として一様乱数が得られるものとします．
プロデルでは引数なしの乱数を求める手順で一様乱数が得られます¹．

なお，Qiita上ではプロデルのシンタックスハイライトが効かないためスクリーンショットを貼り付けています．ソードのテキストなどは GitHub に掲示しています．

サンプリングアルゴリズム

サンプリングとは，ある確率分布関数 $f(x)$ が与えられたときに，それに従う確率変数 $x$ を得る手法のことです．
特に，コンピューターで得やすい $[0,1]$ の一様分布（またはそこから高速に計算できるボックスミュラー法などが知られている正規分布）から，求めたい任意の確率分布関数を与えたときのサンプルを得たいというのがよくあるシチュエーションです．モンテカルロシミュレーションやベイズ推定など，数値的な解析によく用いられます．
よく用いられる有名なものとして

逆関数法
(フォン・ノイマンの)棄却法
マルコフ連鎖モンテカルロ法

などが挙げられます．
本記事ではこれらについてそれぞれ解説していきます．

逆関数法

確率分布関数 $f(x)$ と，その積分である累積分布関数 $F(x)=\int f(x) dx$，およびその逆関数 $F^{-1}(y)$ が計算できるときに使える手法です．

例として

f(x)=\left\{
\begin{array}{ll}
2x & (0 \leq x \leq 1)\\
0 & (\text{otherwize})
\end{array}
\right. , \quad
F(x)=\left\{
\begin{array}{ll}
0 & (x < 0)\\
x^2 & (0 \leq x \leq 1)\\
1 & (x>1)
\end{array}
\right. , \quad
F^{-1}(y)=\sqrt{y}

という三角状の確率分布について計算を行います．

例示用確率分布

実装は以下のようになります．

本題の乱数を生成するアルゴリズムは5～7行目の目的の乱数を求める手順だけです．
すなわち，一様乱数を累積分布関数の逆関数に入れるだけで目的の乱数を得ることが出来ます．
確率密度関数と逆累積分布関数を差し替えれば任意の確率分布を計算することが出来ます．

これは確率変数の変換公式
$$
X\sim f_X(x),\; Y=h(X)\sim f_Y(y)
$$
のとき
$$
f_Y(y)=f_X(h^{-1}(y))\frac{dh^{-1}(y)}{dy}
$$

において $f_X=\mathrm{Unif}(0,1),\; h=F^{-1}$ を適用した結果を利用したものです．

逆関数法の確率変数の変換のイメージ（[2]を元に作成）

ただし，累積分布関数の逆関数を簡単に求めることができるという状況は限られてきます．
重要な例として正規分布は累積分布関数が初等関数で書くことが出来ないため逆関数の値を求めるには数値的に求める必要があり，numpyなどのライブラリを用いても後述する他の方法よりも計算が遅くなってしまいます．
さらに複数の分布が重なった多峰性の分布であったり定義域を制限した切断正規分布などバリエーションが加わるとより計算が困難になります．

棄却法

同様に上に挙げた $f(x)$ でのサンプリングの棄却法での例を示します．

2つの一様乱数を使って矩形領域から目的の分布を抽出する方法です．
この方法では逆関数法と異なり累積分布関数や逆関数を用いず，確率分布関数の順方向の計算だけで実行することが出来ます．

イメージ的には適当な大きめの領域（図青＋赤）のうち確率分布関数の下側の部分（図青）だけを取り出せば求める確率分布になるというものです．

棄却法のイメージ

一方，値域の範囲を求める必要があるためベイズ推定など確率分布関数の規格化計算が困難な場合には用いづらいです．
まだ，定義域からの一様分布を用いていることからわかるように，正規分布のような定義域が実数全体に渡る確率分布ではそのまま適用することは出来ません．
定義域を一定の有限範囲で近似するといった方法が考えられますが，近似精度を上げるため範囲を広げるほど「Y≦（Xの確率密度関数）」の条件をパスする確率が低くなり（全体に対し青の面積が小さくなり），確率的な無限ループがなかなか終わらなくなります．

なお，一般には求める確率分布関数を被覆するような分布を使えばよいため，実数全体に台を持つ分布であってもサンプリングは可能です．
しかし，ループを抜ける確率を保ちつつ容易にサンプリングできる分布を見つけるには目的の分布の形状をよく知っておく必要があり，設計が大変です．
さらに，高次元の変数をサンプリングしようとすると，次元の呪いによりさらにループを抜ける確率が低くなります．

↓冗長ですが， $f(x)$ を正規分布で被覆した場合の実装です．

正規分布での被覆

提案分布乱数，提案分布確率密度関数，提案分布倍率を適切に設定することで一般の場合の棄却法の実装となります．

ネイピア数＝１の自然対数乗
【正規分布乱数器】は，ボックスミュラー乱数生成器を作ったもの
【提案分布倍率】は，２×（（２×π）の平方根）

１０００回，『
	目的の乱数を出力して改行する
』ことを繰り返す

目的の乱数を求める手順
	『
		Xは，提案分布乱数
		Yは，０から（Xの提案分布確率密度関数×提案分布倍率）までの一様乱数
		もしY≦（Xの確率密度関数）ならXを返す
	』ことを繰り返す
終わり

［数値］の，確率密度関数を求める手順
	もし数値が0より小さいまたは数値が１より大きいなら
		０を返す
	そうでなければ
		２×数値を返す
	もし終わり
終わり

提案分布乱数を求める手順
	正規分布乱数器の乱数＋１を返す
終わり

［数値］の，提案分布確率密度関数を求める手順
	(－(数値－1)＾2)の自然対数乗÷((2×π)の平方根)を返す
終わり

［下限］から，［上限］までの，一様乱数を求める手順
	（乱数×（上限－下限）＋下限）を返す
終わり

ボックスミュラー乱数生成器とは
	＋X：浮動小数
	＋Y：浮動小数
	フラグ：真偽値＝×
	乱数を求める手順
		もしフラグならば
			フラグ＝×
			Yを返す
		そうでなければ
			L＝－２とネイピア数で（０から１までの一様乱数）の対数の積の平方根
			W＝２×円周率×（０から１までの一様乱数）
			X=L×cos(W)
			Y=L×sin(W)
			フラグ＝○
			Xを返す
		もし終わり
	終わり
終わり

ｘの自然対数を求める手順
	ネイピア数でｘの対数を返す
終わり

被覆イメージ

青の面積が全体の $\frac{1}{2\sqrt{2\pi}}$ しかないため計算に時間がかかる

マルコフ連鎖モンテカルロ法

上記の手法で現れた欠点を解決したのがマルコフ連鎖モンテカルロ法（MCMC）です．

複雑な逆累積分布関数計算をする必要がない（確率分布関数の計算だけでよい）
確率分布関数の規格化定数を求めなくても良い（確率分布関数の定数倍の関数の計算だけでよい）
高次元でもループを抜ける確率（採択率）を落としづらい

といった利点があります．

$f(x)$ の分布を得たいけれど規格化定数が分からないので尤度関数

\mathcal{L}(x)=\left\{
\begin{array}{ll}
x & (0 \leq x \leq 1)\\
0 & (\text{otherwize})
\end{array}
\right.

だけが得られているとしましょう．

※簡単な関数では馬鹿げた例に思えますが，ベイズモデルなどでは

p(x|\Theta) = \frac{1}{\int p(x|\theta_1, \ldots ,\theta_M)\cdots p(x|\theta_M)\; d\theta_1 \cdots d\theta_M}p(x|\theta_1, \theta_2,\ldots ,\theta_M)p(x|\theta_2,\ldots ,\theta_M)\cdots p(x|\theta_M)

というような高次元の積分を実行しなければ規格化定数が求められない場合も多く，尤度関数だけから計算できることは非常に重要な性質です．

以下はMCMCの中で最も基本的な Metropolis-Hastings 法の実装です．

ただし，MCMCも万能ではなく欠点があります．次のサンプルが前のサンプルの±提案幅の範囲でしか動けないことからも明らかなように，サンプル列は相関を持っています．
提案幅が小さすぎれば大きく動くことが出来ず自己相関は大きくなる一方，大きすぎれば０から１までの一様乱数＜次のサンプルの尤度関数÷前のサンプルの尤度関数の条件をパスする確率が低くなり，逆にまた自己相関が大きくなってしまいます．
提案幅を使わないよう改良されたアルゴリズムもありますが，何れにせよサンプル列に自己相関が発生するので相関が切れるようサンプリングに間隔を開けること，また初期値からの相関が切れ定常分布に収束しているかの検証なども必要です．

サンプルに相関がある様子．提案幅が大きすぎると棄却率が高く値が変化しない一方小さすぎると分布全体を動くことが出来ない．

あとがき

プログラム及び図は GitHub で公開しています．

参考文献

[1] シ「公開コピー誌サンプリング法のノート」暗黒通信団 (2014) 2021-12-18 閲覧, https://ankokudan.org/d/dl/pdf/pdf-sampling.pdf
[2] 樋口さぶろお「逆関数法・棄却法による連続型乱数の生成」 (2019) 2021-12-18 閲覧, https://www.data.math.ryukoku.ac.jp/course/compscib_2019/w10.pdf

ver. 1.7.1090 以降 ↩