Monte Carlo Method

數學, 機率論

介紹 Monte Carlo Method 的採樣法 Inverse Method 與 Rejection Sampling，與 Monte Carlo Method 在期望值與積分估計的應用。

2024/09/13

介紹

Monte Carlo Method 是一種隨機採樣方法，能用於期望值估計與積分估計，核心概念是隨機漫步 (random walk)，這邊介紹幾種 Monte Carlo 方法

隨機採樣

隨機採樣 (Random Sampling) 假設機率分部已知，並透過採樣去分析機率分布的特徵，例如期望值，能應用於機率密度複雜、變數間不獨立或甚至問題不存在 closed form 的情況

介紹幾種常用的方式

Inverse method: 一種透過計算分布反函數以產生採樣的方式，常用於透過均勻分布 $\unif (0, 1)$ 去生成其他分布
Rejection sampling: 一種透過設定拒絕範圍的採樣方式，適用於反函數不存在或分布過於複雜的情況
Importance sampling: 一種調整 rejection sampling 的方案，透過調整權重以提高採樣效率的方案

Inverse Method

我們僅有的隨機生成工具一般只有均勻分布 $\unif (0, 1)$ ，但能透過轉換使得採樣服從指定分布，給定 cdf (cumulative distribution function) $P (x)$ ，我們會用上 cdf 的性質 $P : \bb R \to [0, 1]$ 的反函數 $P^{-1} : [0, 1] \to \bb R$ 並以下定義

\begin{align*} P^{-1} (x) = \inf \{ t : P (t) \geq x \} \end{align*}

使得

\begin{align*} P^{-1} (U) \sim P \end{align*}

其中 $U$ 為均勻 $(0, 1)$ 分布

其優點為計算簡單，缺點則是需計算 $P^{-1}$ ，而並非所有分佈的 $P^{-1}$ 都容易計算，例如常態分佈

演算法

目標：給定 cdf $P$ ，生成服從 $P$ 的採樣
輸入：cdf $P$ 、採樣數量 $n$
輸出：服從 $P$ 的採樣 $\contia{x}{n}$

計算 $P^{-1} (x) = \inf \{ t : P (t) \geq x \}$
生成 $\contia{u}{n} \sim \unif (0, 1)$
計算 $x_i = P^{-1} (u_i)$ ，其中 $i = \conti{n}$

Rejection Sampling

給定 pdf $p (x)$ 且此分佈難以透過 inverse method 採樣，我們可以尋找另一個較容易採樣的分佈 $q (x)$ ，且存在 $c > 0$ 使得對任意 $x$ ， $c q (x) \geq p (x)$ 。在對 $q$ 採樣一個 $x$ 後有機率

\begin{align*} \frac{p(x)}{c q(x)} \end{align*}

接受此採樣，反之則拒絕此次採樣

其優點為能適應複雜的 pdf $p (x)$ ；但缺點是須找到一個已知採樣 $q (x)$ ，且滿足 $c q (x) \geq p (x)$ ，若 $p (x)$ 在 $c p(x)$ 占比不高則會有大量的點被拒絕，導致採樣效率低下

演算法

目標：給定 pdf $P$ ，生成服從 $P$ 的採樣
輸入：pdf $P$ 、採樣數量 $n$
輸出：服從 $P$ 的採樣 $\contia{x}{n}$

找到已知採樣 $q (x)$ 且存在 $c > 0$ 使得任意 $x$ ， $c q (x) \geq p (x)$
生成 $x \sim p$ 與 $u \sim \unif (0, 1)$ ，若 $u \leq p (x) / [c q(x)]$ 則接受此採樣；反之拒絕此採樣
重複步驟 2 直到接受採樣至 $n$ 次

絕對值指數分布生成常態分佈

已知常態分佈的 pdf 為

\begin{align*} p (x) = \frac{1}{\sqrt{2 \pi}} \exp \left( - \frac{x^2}{2} \right) \end{align*}

其 cdf 不存在 closed form，因此難以使用 inverse method 採樣；但能透過 rejection sampling 生成，給定絕對值指數分佈 pdf 與 cdf

\begin{align*} q (x) & = \frac{1}{2} e^{-|x|} \\ Q (x) & = \begin{cases} 1 - \frac{1}{2} e^{-x} & \text{if } x \geq 0 \\ \frac{1}{2} e^{x} & \text{otherwise} \end{cases} \end{align*}

透過 inverse method 轉換

\begin{align*} Q^{-1} (x) & = \begin{cases} \ln (2x) & \text{if } 0 < x \leq \frac{1}{2} \\ - \ln (2 - 2x) & \text{if } \frac{1}{2} < x < 1 \end{cases} \end{align*}

以此獲得滿足 $q$ 分佈的採樣，接著尋找 $c > 0$ 使得 $c q(x) \geq p (x)$ ，等價於 $c q (x) / p (x) \geq 1$

\begin{align*} 1 \leq \frac{c q (x)}{p (x)} = \frac{\frac{c}{2} \exp (- |x|)}{\frac{1}{\sqrt{2 \pi}} \exp(-\frac{x^2}{2})} = c \sqrt{\frac{\pi}{2}} \exp \left( \frac{(|x| - 1)^2 - 1}{2} \right) \end{align*}

即

\begin{align*} c \geq \sqrt{\frac{2}{\pi}} \exp \left( - \frac{(|x| - 1)^2 - 1}{2} \right) \end{align*}

因上式須對任意 $x$ 均成立，故

\begin{align*} c \geq \max_{x \in \bb R} \sqrt{\frac{2}{\pi}} \exp \left( - \frac{(|x| - 1)^2 - 1}{2} \right) = \sqrt{\frac{2e}{\pi}} \approx 1.3 \end{align*}

接著生成 $x \sim p$ 與 $u \sim \unif (0, 1)$ ，若 $u \leq p (x) / [c q(x)]$ 則接受此採樣；反之拒絕此採樣，並重複此步驟直到接受採樣至 $n$ 次

期望值估計

給定隨機變數 $X$ 與其 pdf $p (x)$ ，再給定函數 $f (x)$ ，試計算期望值 $E_{p} [f (X)]$ ，可以透過採樣 $\contia{x}{n} \sim p$ 與計算樣本平均

\begin{align*} \hat E_{p} [f (X)] = \frac{1}{n} \sum_{i = 1}^{n} f (x_i) \end{align*}

在大數法則 (law of large number) 下，此方式是收斂的

\begin{align*} \hat E_{p} [f (X)] \to E_{p} [f (X)] \quad \text{as} \quad n \to \infty \end{align*}

Monte Carlo Integration

給定積分函數 $h (x)$

\begin{align*} \int_{\mathcal X} h (x) dx \end{align*}

若 $h (x)$ 能被分解成一個 pdf $p (x)$ 與另一函數 $f (x)$ 使得 $h (x) = p (x) f (x)$ ，則

\begin{align*} \int_{\mathcal X} h (x) dx = \int_{\mathcal X} f (x) p (x) dx = E_p [f (x)] \end{align*}

將此問題轉換成期望值問題，即可以用期望值估計處理

參考資料

LI, Hang. Machine Learning Methods. Springer Nature, 2023.
Lange, Kenneth, J. Chambers, and W. Eddy. Numerical analysis for statisticians. Vol. 1. New York: springer, 2010.