Monte Carlo Method
介紹
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.