所谓蒙特卡罗方法（Monte Carlo method），也称为统计模拟方法，指的是一系列随机模拟某个分布，然后近似计算某些量的方法。蒙特卡罗方法在金融，计算物理，机器学习等领域有着广泛的应用。蒙特卡罗方法的命名来自于大数学家冯诺依曼(John von Neumann)，其将该算法命名为一家摩纳哥的世界著名赌场的名字，给该算法增加了一层神秘色彩。
1. 逆概率变换法（inverse probability transform）。
蒙特卡罗方法的核心问题其实是给定一个分，如何按照该分布模拟生成随机数(向量)，我们先来看一个最简单情况。已知一个取值非负的函数（实际问题中一般是连续函数）$p:mathbb{R}longrightarrowmathbb{R}^{+}$，如果其满足$int_{-infty}^{+infty}p(x)dx=1$，那么如何模拟生成随机数$X$使得$Xsim p$?
一种最简单的办法是直接利用分布函数（Cumulative Distribution Function, 简称CDF）F的反函数，这种方法称为逆概率变换法（inverse probability transform）。
首先，计算机可以利用线性同余发生器（Linear congruential generator）生成伪随机数, 来模拟均匀分布$U([0,1])$。下面我们看下面的一个简单结论：
命题1.1. 如果$F$是某个概率分布的分布函数，其反函数$F^{-1}$存在，如果随机变量$Usim U([0,1])$，则$F^{-1}(U)sim F$。
以上结论从下面等式立马可以被观察到，当$Usim U([0,1])$的时候，对于任意$cinmathbb{R}$：
egin{equation}P(F^{-1}(U)leq c)=P(Uleq F(c))=F(c),end{equation}
所以我们只要知道了某个分布函数的反函数，则模拟该分布生成随机数就很容易了，所幸的是常见的分布族求解反函数都还比较容易。
2. Box-Muller变换生成模拟正态分布
我们现在想模拟二维标准正态分布，其概率密度是$p(x,y)=frac{1}{2pi}e^{-frac{x^{2}+y^{2}}{2}}$，对任意的$(x,y)inmathbb{R}^{2}$。除了上面的逆概率变换法外还有一种更加方便的算法。
我们用极坐标表示$mathbb{R}^{2}$中的点：$x=rcos heta$,$y=rsin heta$, $rin[0,infty)$, $ hetain[0,2pi)$, 概率计算的时候的积分项（测度）在极坐标下有：
egin{equation}egin{split}&frac{1}{2pi}e^{-frac{x^{2}+y^{2}}{2}} ext{d}xcdot ext{d}y ewline =&frac{1}{2pi}e^{-frac{r^{2}}{2}}r ext{d}rcdot ext{d} heta ewline =&e^{-frac{r^{2}}{2}}r ext{d}rcdot ext{d}(frac{ heta}{2pi}) ewline =& ext{d}(e^{-frac{r^{2}}{2}})cdot ext{d}(frac{ heta}{2pi})end{split},end{equation}
所以我们做一个变量替换:egin{equation}egin{split}u=&e^{-frac{r^{2}}{2}} ewline v=&frac{ heta}{2pi},end{split}end{equation}
或者等价于: egin{equation}egin{split}x=&sqrt{-2log(u)}cos(2pi v) ewline y=&sqrt{-2log(u)}sin(2pi v)end{split}end{equation}
的时候，我们必然有：egin{equation} ext{d}ucdot ext{d}v=frac{1}{2pi}e^{-frac{x^{2}+y^{2}}{2}} ext{d}xcdot ext{d}y.end{equation}
由上面的式子我们立即知道，如果定义一个变换：
$$BM: [0,1] imes [0,1]longrightarrow mathbb{R}^{2},$$
egin{equation}BM(u,v) riangleq (sqrt{-2log(u)}cos(2pi v), sqrt{-2log(u)}sin(2pi v)),end{equation}
则对于任意独立的$U_{1}sim U([0,1])$， $U_{2}sim U([0,1])$，$(V_{1}, V_{2}) riangleq BM(U_{1}, U_{2})sim mathcal{N}(0, ext{I}).$
3.拒绝-接受算法：
3.1. 拒绝-接受算法的推导
拒绝-接受算法的大体思路是先模拟某个比较容易模拟的提议分布(proposal probability)生成随机数（向量），然后用一定的概率来二次把关，决定是保留还是不要该数（向量），以达到保留该数（向量）的条件下该数（向量）的分布为预期要模拟的分布。
具体来说是这样：
假设给定了某个待模拟的$K$维分布的概率密度函数$p$, 现在选择一个容易模拟的分布，其概率密度$q$，生成一个随机数$X$，然后生成一个$Ysim U([0,1])$, 两过程完全独立，我们想设定一个接受$X$之概率$r(X)$，也就是当$Yleq r(X)$的时候才接受$X$，否则拒绝。
这个时候我们定义事件：egin{equation}A=lbrace Yleq r(X) brace,end{equation}
该事件就是"采样被接受"。我们希望$X$在$A$上的条件概率密度为$p$。也就是希望：
egin{equation}P(Xleq cmid A)=int_{-infty}^{c}p(x) ext{d}x,end{equation}
对于任意的$cinmathbb{R}^{K}$，其中$ ext{''}leq ext{''}$表示每个分量都不大于。
我们容易计算：
egin{equation}egin{split}P(Xleq cmid A)=&P(lbrace Xleq c,Yleq r(X) brace)/P(lbrace Yleq r(X) brace) ewline =&int_{lbrace (x,y)mid xleq c,yleq r(x) brace}Mq(x) ext{d}x ext{d}y ewline =&int_{-infty}^{c}M q(x)(int_{0}^{r(x)} ext{d}y) ext{d}x ewline =&int_{-infty}^{c}Mq(x)r(x) ext{d}x,end{split}end{equation}
其中$M=(int_{mathbb{R}^{K}}q(x)r(x) ext{d}x)^{-1}$, 所以我们有：
$$p(x)=M q(x)r(x),$$
所以我们立即得到：egin{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)如果$uleq frac{p(x)}{Mq(x)}$，令$x_{0}=z$，否则回到步骤1)继续，直到得到一个被接受的向量。
3.2. 拒绝-接受采样的缺陷：
拒绝接受采样的缺陷很明显，那就是当$M$过大的时候，算法的效率极其低下。因为每次试验生成的随机数或者向量被接受的概率仅仅为$frac{1}{M}$，所以我们必须平均生成$M$个数(向量)才能最终成功采样其中一个。当我们采样的维度$K$非常大的时候(例如在统计物理中的模拟问题)，往往首先$q$难以寻找，齐次就算找到$M$相应也会很大，显然拒绝-接受采样不合适，所以一般它只适用于低维分布。
4.重要性采样
现在假设某$K$维随机变量$X$有概率密度$p$，我们现在考虑用一种蒙特卡罗方法计算期望，也就是右边的积分：
egin{equation}E(f(X))=int_{mathbb{R}^{K}}f(x)p(x) ext{d}x.end{equation}
现在如果选取一个提议分布，密度函数为$q$，则对于任意的$ ilde Xsim q,$
egin{equation}egin{split}E(f(X))=&int_{mathbb{R}^{K}}f(x)p(x) ext{d}x ewline =&int_{mathbb{R}^{K}}frac{f(x)p(x)}{q(x)}q(x) ext{d}x ewline =&E(frac{f( ilde X)p( ilde X)}{q( ilde X)}) ewline =&E(f( ilde X)w( ilde X)),end{split}end{equation}
在这里我们定义函数$w riangleq 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$。
现在我们来考察一下方差egin{equation}delta_{q}^{2}=Var(f( ilde X)w( ilde X))end{equation}，因为他直接影响着$Ilongrightarrow E(f(X))$的收敛速率。由中心极限定理，当$N$足够大的时候$frac{I-E(f)}{delta_{q}/sqrt{N}}$近似服从标准正态分布$mathcal{N}(0, ext{I})$，由此我们可以近似给出$E(f(X))$的一个近似$1-alpha$置信区间：
egin{equation}(I-frac{delta_{q}}{sqrt{N}}z_{alpha/2},I+frac{delta_{q}}{sqrt{N}}z_{alpha/2}),end{equation}
而$Var(f( ilde X)w( ilde X))=E(f^{2}( ilde X)w^{2}( ilde X))-E(f( ilde X))^{2}$，于是我们只需要估计一下$E(f^{2}( ilde X)w^{2}( ilde X)).$由Cauchy-Schwarz不等式：
egin{equation}E(f^{2}( ilde X)w^{2}( ilde X))geq (E(f( ilde X)w( ilde X)))^{2}=E(f(X))^{2},end{equation}
等号成立当且仅当：$vert fwvertequiv ext{const}$，也就是：egin{equation}q(x)=q^{ast}(x) riangleqfrac{vert f(x)vert p(x)}{int_{mathbb{R}^{K}}vert f(u)vert p(u) ext{d}u}end{equation}
的时候$E(f^{2}( ilde X)w^{2}( ilde 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) ext{d}u}$。但是往往这是很难办到的，尤其是分布的维数$K$很大的时候，该算法收敛可能还是会出现比较慢的问题。
下面我们看一个例子，看看提议分布的选取如何影响收敛速度，该例子取自于[1]的例25.6。
例4.1 如上图，我们令概率密度$p(x)=(2pi)^{-frac{1}{2}}e^{-frac{x^{2}}{2}}$，表示标准正态分布的概率密度函数，$f(x)=100$如果$xgeq 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)=(2pi)^{-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