蒙特卡罗模拟

背景

随机模拟也可以叫做蒙特卡罗模拟(Monte Carlo Simulation)。这个方法的发展始于20世纪40年代，和原子 弹制造的曼哈顿计划密切相关，当时的几个大牛，包括乌拉姆、冯.诺依曼、费米、费曼、Nicholas Metropolis， 在美国洛斯阿拉莫斯国家实验室研究裂变物质的中子连锁反应的时候，开始使用统计模拟的方法 并在最早的计算机上进行编程实现。

均匀分布

随机模拟中有一个重要的问题就是给定一个概率分布p(x)，我们如何在计算机中生成它的样本。一般而言均匀分 布 Uniform(0,1)的样本是相对容易生成的。 通过线性同余发生器可以生成伪随机数，我们用确定性算法生成 [0,1]之间的伪随机数序列后，这些序列的各种统计指标和均匀分布 Uniform(0,1) 的理论计算结果非常接近。 这样的伪随机序列就有比较好的统计性质，可以被当成真实的随机数使用。例如可以如下方式生成：

$$x_{n+1} =(ax_{n}+c) \ mod \ m$$

正太分布

利用获得的均匀分布样本，Box-Muller变换可以生成满足正太分布的样本。
Box-Muller 变换
如果随机变量 $U1, U2$ 独立且 $U1,U2 \sim Uniform(0,1)$,

$$Z_{0}=\sqrt{-2\ln U_{1}}\cos(2\pi U_{2})$$ $$Z_{1}=\sqrt{-2\ln U_{1}}\sin(2\pi U_{2})$$ 则 $Z_{0}, Z_{1}$ 独立且 $Z_{0}, Z_{1} \sim N(0,1)$ .

Matlab代码实现如下：

for i =1:10000
    U1 = rand;
    U2 = rand;

    Z0 (i) = sqrt((-2)*log(U1))*cos(2 * pi * U2);
    Z1 (i) = sqrt((-2)*log(U1))*sin(2 * pi * U2);
end

subplot(2,1,1); 
histfit(Z0)   
xlabel('Z0'); 
ylabel('出现的次数'); 

subplot(2,1,2); 
histfit(Z1)   
xlabel('Z1'); 
ylabel('出现的次数');   

实验结果：

Monte Carlo数值积分

Monte Carlo 抽样计算随机变量的期望值是接下来内容的重点：$z$ 表示随机变量，服从概率分布 $p(z)$，计算 $p(z)$ 的期望， 只需要我们不停地从 $p(z)$ 中抽样 $z_{i}$，然后对这些 $f(z_{i})$ 取平均值即可近似 $f(z)$ 的期望。

$$E_{N}(f)=\frac{1}{N}\sum_{i=1}^{N}f(z_{i})$$

求解定积分

要求解 $f(x)$ 在区间 $[a,b]$ 上的定积分：$\int_{a}^{b} f(x)$，如果$f(x)$ 形式比较复杂不好求解，就可以从过数值模拟的方法得到近似解，如下所示。

$$\int_{a}^{b} \frac{f(x)}{p(x)} p(x) \ dx$$
其中 $p(x)$ 为区间 $[a,b]$ 上的任一概率分布（密度函数），即满足 $\int_{a}^{b} p(x) \ dx = 1$. 上式等同于求解随机变量 $\frac{f(x)}{p(x)}$ 的期望。根据大数定理，可以根据 $p(x)$ 抽取 $N$ 个样本，当样本数足够大时，可以用样本均值来近似期望值。
$$\int_{a}^{b} \frac{f(x)}{p(x)} p(x) \ dx \approx \frac{1}{N} \sum_{i=1}^{N}\frac{f(x_{i})}{p(x_{i})}$$

例子

如果我们要求解定积分 $\int_{1}^{2}x^{3} \ dx$，采用的概率分布为区间 $[1,2]$ 上的均匀分布，概率密度函数为：

$$p(x) = \frac{1}{3-1}= \frac{1}{2}$$
matlab 代码

a = 1;
b = 3;

sum = 0;
N =10000;

for i=1:N
    s = a + (b-a) .* rand;
    p = 1 / (b-a);
    sum = sum + s*s / p;
end
sum/N

实验结果为：

ans =

    8.6403

非常接近于精确解 $\frac{26}{3}=8.6666...$