蒙特卡洛积分

该文是对《Stochastic Numerical Methods: An Introduction for Students and Scientists》的总结。

Hit and miss法

给定一个区域\(\Omega\)上的非负函数\(f\)，如果 $0 \leq f \leq c $ ，则可以考虑一个随机变量 \((x, y)\) ，它均等可能地取值于 \(\Omega \times [0,c]\) ，则\(y \leq f(x)\)的概率是\(I/(m(\Omega) \times c)\)，其中\(I = \int_{\Omega} f\)。通过取点并计算\(y \leq f(x)\)的比例，可以得到对积分的估计。
或者认为这里有一个伯努利随机变量，取1(\(y \leq f(x)\))的概率为\(p=I/(m(\Omega) \times c)\)，该随机变量均值是\(p\)。通过实验估计均值，我们可以得到对积分\(I\)的区间概率估计。
该方法的误差以\(M^{-1/2}\)的速度下降，其中\(M\)是随机变量的实验总数。

Uniform sampling法

该方法是重要性采样的特例，只要随机变量在区域内均匀取值即可。

重要性采样

给定一个区域\(\Omega\)上的函数\(g\)，我们想要计算\(I = \int_{\Omega} g\)。可以将积分写为\(I = \int_{\Omega} Gf\)，其中\(f\)是某个随机变量的pdf，并且\(Gf=g\)。如果找到合适的pdf函数\(f\)，我们可以将积分视为期望\(E[G(\xi)]\)，\(\xi\)服从pdf\(f\)。这样我们可以用样本均值估计积分，且误差按\(M^{-1/2}\)的速度下降。
当\(f\)取为区域\(\Omega\)上的常量函数时（由于pdf的要求，这个常量大于0），即为uniform sampling方法。考虑到Riemann积分的定义，可以认为这是一个对取点更随意、误差下降更慢的按定义计算方法。
重要性采样并未指定\(f\)的选择方法，可以证明，使误差最小的\(f\)和\(|g|\)有相同的形状，即只要将\(|g|\)归一化，就是\(f\)的最优选择，这时\(G=c \times sign(g)\)。实际应用中，完全和\(|g|\)形状一致是很难的，需要在形状近似和易于生成之间权衡。
重要性采样也可以用于离散数列的求和，无论是单变量或多变量的\(\sum\)。
重要性采样需要根据指定pdf生成值，如何生成可以参见https://www.cnblogs.com/arcueid-take-me-go/p/18703874。

带重复的拒绝法

对于高维分布，通过拒绝法生成随机变量时，接受的概率\(\epsilon\)可能会非常小，以致于需要很多次尝试才能获得一个值。如果只关注积分的问题，带重复的拒绝法仍能获得\(M^{-1/2}\)的误差下降速度，而它有效率得多。
带重复的拒绝法步骤如下：

1. 根据拒绝法的要求选取函数\(g、h\)。
2. 选取某个初始值\(x_0\)，或者从某个分布中生成一个初值。
3. 已有\(x_i\)，根据拒绝法生成一个随机变量的值\(y\)
4. 根据拒绝法，如果接受生成的值，令\(x_{i+1}=y\)；如果拒绝，令\(x_{i+1}=x_i\)
5. 重复3-5步

该方法获得的点列不是相互独立的，并且不一定遵循希望的分布。但可以证明，\(x_i\)遵循的分布\(f_i\)和期望的分布\(f\)之间的距离，以指数形式下降；且拒绝法的接受概率越大，下降越快。因此，经过预热（舍弃前部分项）后，可以认为近似服从分布\(f\)。尽管点列之间有相关性，但仍然可以用于计算积分，且误差下降的速度不变；但标准差被放大为\(\sqrt{ \frac{2-\epsilon}{\epsilon}}\)倍，因此仍然要求拒绝法的接受概率不能太小，否则会影响效率。

MCMC法（Markov Chain-Monte Carlo Method）

这里的核心定理如下：

定理： 假设\(p(u, w)\)是一个马尔可夫链的概率转移函数，\(\mu\)是该马尔科夫链的不变测度，其pdf为\(\rho\)。如果该马尔可夫链是遍历的(ergodic)，则对任意有界连续函数\(\varphi\)，

\[\frac {1}{N} \sum_{n=1}^{N} \varphi (u^n) \rightarrow E^\mu \varphi \]

对\(\mu\)几乎处处成立。特别地，如果存在概率\(p\)以及\(\epsilon > 0\)使得，对任意\(u\)和Borel集\(A\)，\(p(u, A) \geq \epsilon p(A)\)，则对任意\(u\)，我们有

\[d_{TV} ( p^n (u, \centerdot ), \mu ) \leq 2(1-\epsilon )^n \]

并且，此时存在\(K=K(\varphi )>0\)，使得

\[\frac {1}{N} \sum_{n=1}^{N} \varphi (u^n) = E^\mu \varphi + K \xi_N N^{-1/2} \]

其中，\(xi_N\)弱收敛到\(N(0,1)\)。

这个定理的强大之处在于，如果我们想要在重要性采样中计算积分，我们只需要选择好pdf\(f\)，然后构造一个使该pdf是不变测度的马尔科夫链，然后随便选个值，根据马尔科夫链得到一系列值，将这些值代入函数\(G\)，求个平均就好了！你知道的，遍历性常常会被满足。
定理的其他部分告诉我们，如果该马尔可夫系统从任意状态出发，转移到某个其它区域（集合）的概率都能被某个与初始状态无关的概率支撑，那么状态转移给出的概率分布会逐渐靠近不变测度，并且误差按指数下降。定理的剩下部分告诉了我们计算积分时，误差下降的速度。你可以想象在这种情况下，我们在观测某个马尔可夫链，我们不知道初值是什么，但是我们在状态空间某个区域看到它的概率就是不变测度\(\mu\)！你有一种感觉，只有这个粒子的轨迹服从不变测度\(\mu\)的分布时，这才能实现。而既然它的轨迹就是我们感兴趣的分布，拿它计算能得到积分也就不足为奇了——反正代入计算平均值的随机变量取值和其它方法获得的随机变量取值是一样的，最多顺序不同而已。

我们剩下的问题就是，给定测度\(\mu\)，如何构造马尔可夫链使其为不变测度？