所謂蒙特卡羅方法（Monte Carlo method），也稱為統計模擬方法，指的是一系列隨機模擬某個分布，然后近似計算某些量的方法。蒙特卡羅方法在金融，計算物理，機器學習等領域有着廣泛的應用。蒙特卡羅方法的命名來自於大數學家馮諾依曼(John von Neumann)，其將該算法命名為一家摩納哥的世界著名賭場的名字，給該算法增加了一層神秘色彩。
1. 逆概率變換法（inverse probability transform）。
蒙特卡羅方法的核心問題其實是給定一個分，如何按照該分布模擬生成隨機數(向量)，我們先來看一個最簡單情況。已知一個取值非負的函數（實際問題中一般是連續函數）$p:\mathbb{R}\longrightarrow\mathbb{R}^{+}$，如果其滿足$\int_{-\infty}^{+\infty}p(x)dx=1$，那么如何模擬生成隨機數$X$使得$X\sim p$?
一種最簡單的辦法是直接利用分布函數（Cumulative Distribution Function, 簡稱CDF）F的反函數，這種方法稱為逆概率變換法（inverse probability transform）。
首先，計算機可以利用線性同余發生器（Linear congruential generator）生成偽隨機數, 來模擬均勻分布$U([0,1])$。下面我們看下面的一個簡單結論：
命題1.1. 如果$F$是某個概率分布的分布函數，其反函數$F^{-1}$存在，如果隨機變量$U\sim U([0,1])$，則$F^{-1}(U)\sim F$。
以上結論從下面等式立馬可以被觀察到，當$U\sim U([0,1])$的時候，對於任意$c\in\mathbb{R}$：
\begin{equation}P(F^{-1}(U)\leq c)=P(U\leq F(c))=F(c),\end{equation}
所以我們只要知道了某個分布函數的反函數，則模擬該分布生成隨機數就很容易了，所幸的是常見的分布族求解反函數都還比較容易。
2. Box-Muller變換生成模擬正態分布
我們現在想模擬二維標准正態分布，其概率密度是$p(x,y)=\frac{1}{2\pi}e^{-\frac{x^{2}+y^{2}}{2}}$，對任意的$(x,y)\in\mathbb{R}^{2}$。除了上面的逆概率變換法外還有一種更加方便的算法。
我們用極坐標表示$\mathbb{R}^{2}$中的點：$x=r\cos\theta$,$y=r\sin\theta$, $r\in[0,\infty)$, $\theta\in[0,2\pi)$, 概率計算的時候的積分項（測度）在極坐標下有：
\begin{equation}\begin{split}&\frac{1}{2\pi}e^{-\frac{x^{2}+y^{2}}{2}}\text{d}x\cdot\text{d}y\newline =&\frac{1}{2\pi}e^{-\frac{r^{2}}{2}}r\text{d}r\cdot\text{d}\theta\newline =&e^{-\frac{r^{2}}{2}}r\text{d}r\cdot\text{d}(\frac{\theta}{2\pi})\newline =&\text{d}(e^{-\frac{r^{2}}{2}})\cdot\text{d}(\frac{\theta}{2\pi})\end{split},\end{equation}
所以我們做一個變量替換:\begin{equation}\begin{split}u=&e^{-\frac{r^{2}}{2}}\newline v=&\frac{\theta}{2\pi},\end{split}\end{equation}
或者等價於: \begin{equation}\begin{split}x=&\sqrt{-2\log(u)}\cos(2\pi v)\newline y=&\sqrt{-2\log(u)}\sin(2\pi v)\end{split}\end{equation}
的時候，我們必然有：\begin{equation}\text{d}u\cdot\text{d}v=\frac{1}{2\pi}e^{-\frac{x^{2}+y^{2}}{2}}\text{d}x\cdot\text{d}y.\end{equation}
由上面的式子我們立即知道，如果定義一個變換：
$$BM: [0,1]\times [0,1]\longrightarrow \mathbb{R}^{2},$$
\begin{equation}BM(u,v)\triangleq (\sqrt{-2\log(u)}\cos(2\pi v), \sqrt{-2\log(u)}\sin(2\pi v)),\end{equation}
則對於任意獨立的$U_{1}\sim U([0,1])$， $U_{2}\sim U([0,1])$，$(V_{1}, V_{2})\triangleq BM(U_{1}, U_{2})\sim \mathcal{N}(0,\text{I}).$
3.拒絕-接受算法：
3.1. 拒絕-接受算法的推導
拒絕-接受算法的大體思路是先模擬某個比較容易模擬的提議分布(proposal probability)生成隨機數（向量），然后用一定的概率來二次把關，決定是保留還是不要該數（向量），以達到保留該數（向量）的條件下該數（向量）的分布為預期要模擬的分布。
具體來說是這樣：
假設給定了某個待模擬的$K$維分布的概率密度函數$p$, 現在選擇一個容易模擬的分布，其概率密度$q$，生成一個隨機數$X$，然后生成一個$Y\sim U([0,1])$, 兩過程完全獨立，我們想設定一個接受$X$之概率$r(X)$，也就是當$Y\leq r(X)$的時候才接受$X$，否則拒絕。
這個時候我們定義事件：\begin{equation}A=\lbrace Y\leq r(X)\rbrace,\end{equation}
該事件就是"采樣被接受"。我們希望$X$在$A$上的條件概率密度為$p$。也就是希望：
\begin{equation}P(X\leq c\mid A)=\int_{-\infty}^{c}p(x)\text{d}x,\end{equation}
對於任意的$c\in\mathbb{R}^{K}$，其中$\text{''}\leq\text{''}$表示每個分量都不大於。
我們容易計算：
\begin{equation}\begin{split}P(X\leq c\mid A)=&P(\lbrace X\leq c,Y\leq r(X)\rbrace)/P(\lbrace Y\leq r(X)\rbrace)\newline =&\int_{\lbrace (x,y)\mid x\leq c,y\leq r(x)\rbrace}Mq(x)\text{d}x\text{d}y \newline =&\int_{-\infty}^{c}M q(x)(\int_{0}^{r(x)}\text{d}y)\text{d}x\newline =&\int_{-\infty}^{c}Mq(x)r(x)\text{d}x,\end{split}\end{equation}
其中$M=(\int_{\mathbb{R}^{K}}q(x)r(x)\text{d}x)^{-1}$, 所以我們有：
$$p(x)=M q(x)r(x),$$
所以我們立即得到：\begin{equation}r(x)=\frac{p(x)}{M q(x)},\end{equation}
別忘記了$r(x)\in [0,1],$所以$M$還得足夠大使得：$Mq(x)\geq p(x)$。
綜上所述我們總結一下接受-拒絕采樣得到一個向量的方法：
算法(接受-拒絕采樣)：
1)選擇一個合適的，容易模擬的具有概率密度$q(x)$的提議分布，選擇一個足夠大的$M$使得$Mq(x)\geq p(x)$；
2)模擬提議分布，生成$z$；
3)從均勻分布$U([0,1])$抽樣得到$u$；
4)如果$u\leq \frac{p(x)}{Mq(x)}$，令$x_{0}=z$，否則回到步驟1)繼續，直到得到一個被接受的向量。
3.2. 拒絕-接受采樣的缺陷：
拒絕接受采樣的缺陷很明顯，那就是當$M$過大的時候，算法的效率極其低下。因為每次試驗生成的隨機數或者向量被接受的概率僅僅為$\frac{1}{M}$，所以我們必須平均生成$M$個數(向量)才能最終成功采樣其中一個。當我們采樣的維度$K$非常大的時候(例如在統計物理中的模擬問題)，往往首先$q$難以尋找，齊次就算找到$M$相應也會很大，顯然拒絕-接受采樣不合適，所以一般它只適用於低維分布。
4.重要性采樣
現在假設某$K$維隨機變量$X$有概率密度$p$，我們現在考慮用一種蒙特卡羅方法計算期望，也就是右邊的積分：
\begin{equation}E(f(X))=\int_{\mathbb{R}^{K}}f(x)p(x)\text{d}x.\end{equation}
現在如果選取一個提議分布，密度函數為$q$，則對於任意的$\tilde X\sim q,$
\begin{equation}\begin{split}E(f(X))=&\int_{\mathbb{R}^{K}}f(x)p(x)\text{d}x\newline =&\int_{\mathbb{R}^{K}}\frac{f(x)p(x)}{q(x)}q(x)\text{d}x\newline =&E(\frac{f(\tilde X)p(\tilde X)}{q(\tilde X)})\newline =&E(f(\tilde X)w(\tilde X)),\end{split}\end{equation}
在這里我們定義函數$w\triangleq \frac{p}{q}.$
稱$w(X_{i})$為$X_{i}$的重要性權重(importance weights)。令$I=\sum_{i=1}^{N}f(X_{i})w(X_{i})$，由大數定理，$I$幾乎處處收斂於$E(f(X_{1})w(X_{1}))=E(f(X))$，所以我們可以取足夠大的$N$, 這時就可以近似計算得到：$E(f)\approx I.$。我們總結一下該算法（重要性采樣，Importance Sampling）：
重要性采樣：
1)選擇某個容易模擬的提議分布，具有概率密度$q$，然后模擬該分布生成$x_{1},...,x_{N}$；
2)計算$I=\sum_{i=1}^{N}f(x_{i})p(x_{i})/q(x_{i})$，得到近似估計值$E(f)\approx I$。
現在我們來考察一下方差\begin{equation}\delta_{q}^{2}=Var(f(\tilde X)w(\tilde X))\end{equation}，因為他直接影響着$I\longrightarrow E(f(X))$的收斂速率。由中心極限定理，當$N$足夠大的時候$\frac{I-E(f)}{\delta_{q}/\sqrt{N}}$近似服從標准正態分布$\mathcal{N}(0,\text{I})$，由此我們可以近似給出$E(f(X))$的一個近似$1-\alpha$置信區間：
\begin{equation}(I-\frac{\delta_{q}}{\sqrt{N}}z_{\alpha/2},I+\frac{\delta_{q}}{\sqrt{N}}z_{\alpha/2}),\end{equation}
而$Var(f(\tilde X)w(\tilde X))=E(f^{2}(\tilde X)w^{2}(\tilde X))-E(f(\tilde X))^{2}$，於是我們只需要估計一下$E(f^{2}(\tilde X)w^{2}(\tilde X)).$由Cauchy-Schwarz不等式：
\begin{equation}E(f^{2}(\tilde X)w^{2}(\tilde X))\geq (E(f(\tilde X)w(\tilde X)))^{2}=E(f(X))^{2},\end{equation}
等號成立當且僅當：$\vert fw\vert\equiv \text{const}$，也就是：\begin{equation}q(x)=q^{\ast}(x)\triangleq\frac{\vert f(x)\vert p(x)}{\int_{\mathbb{R}^{K}}\vert f(u)\vert p(u)\text{d}u}\end{equation}
的時候$E(f^{2}(\tilde X)w^{2}(\tilde X))$取極小值$E(f(X))^{2}$，此時$\delta_{q}^{2}$取最小值$Var(f(X))$。所以為了保證我們能用更小的生成向量個數$N$以得到足夠小的置信區間，我們理想狀態是選擇某個提議分布使得其概率密度$q(x)$接近於$\frac{\vert f(x)\vert p(x)}{\int_{\mathbb{R}^{K}}\vert f(u)\vert p(u)\text{d}u}$。但是往往這是很難辦到的，尤其是分布的維數$K$很大的時候，該算法收斂可能還是會出現比較慢的問題。
下面我們看一個例子，看看提議分布的選取如何影響收斂速度，該例子取自於[1]的例25.6。
例4.1 如上圖，我們令概率密度$p(x)=(2\pi)^{-\frac{1}{2}}e^{-\frac{x^{2}}{2}}$，表示標准正態分布的概率密度函數，$f(x)=100$如果$x\geq 3$, 否則$f(x)=0$。這時候如果我們先直接取$q=p$，生成隨機數$X_{i}\sim p$，$i=1,...,N$，$N=100$來近似計算$E(f(X))= 10P(X > 3) =0.0013$，看看會發生什么。首先，如果$N$取得太小，例如是這里的100,那么我們抽樣得到右邊尾部$[3,\infty)$的概率極其低，幾乎抽不到這里的點，算出來的$I$大概率是$0$。我們再計算一下得到$Var(I)=0.039$。這時這種抽樣效率顯然過低下。這時候其實最好是取一個概率分布使得尾部$[3,\infty)$的概率顯著增大，而$Var(I)$顯著減小，例如如果我們取$q(x)=(2\pi)^{-\frac{1}{2}}e^{-\frac{(x-4)^{2}}{2}}$，這時$Var(I)=0.0002$，方差大大降低了，收斂的速率也就越高。一般來說，在高維的時候，如果密度函數$p$，$f$很復雜的時候，尋找$q$就沒那么容易了，高維度的困難似乎是這里面提到的所有蒙特卡羅方法的通病。
5.參考文獻：
[1]Larry Wasserman, All of Statistics: A Concise Course in Statistical Inference,
[2]Kevin P. Murphy, Machine Learning A Probabilistic Perspective,2012 Massachusetts Institute of Technology