【备战美赛】蒙特卡洛算法

蒙特卡洛算法是一种基于随机数的算法，思想非常直观但是与其他算法相比精确度较低。其实我们在中学已经初步了解过蒙特卡洛算法的应用，例如通过在坐标轴中随机取点来估算圆周率pi的值。蒙特卡洛算法的其他应用本质上也和随机取点估值问题类似。

数值方法求定积分

对于定积分 $\int _a ^bf(x)dx$ ，一般我们会找原函数然后通过牛顿莱布尼兹公式求解。但是遇到一些复杂的表达式可能就很难求出原函数。例如计算定积分： $I=\int _{0.8}^{3.3}\frac{1}{1+\sin{x}+ln^2x}dx$ ，很难通过一般的方法求解。这时候就要用到蒙特卡洛算法计算出一个近似的值。

因为： $\int _a^bf(x)dx=\displaystyle\lim_{n\rightarrow \infty}\frac{b-a}{n}\displaystyle\sum_{i=1}^nf(x_i)$ ，我们可以在积分区间[a,b]上均匀地随机生成一列xi，之后求出积分的近似值。比如计算 $I=\int _{0.8}^{3.3}\frac{1}{1+\sin{x}+ln^2x}dx$ ，就可以先在[0.8,3.3]上随机生成n个点，当n足够大时，便可以用近似值代替估算值。

Matlab代码如下：

clc
clear all;
close all;
syms x
f=1/(1+sin(x)+(log(x))^2)
n=1e3
x=unifrnd(0.8,3.3,[1,n]);
ans=0;
for i=1:ny(i)=1/(1+sin(x(i))+(log(x(i))^2))ans=ans+y(i)
end
ans=2.5*ans/n

最终运行结果是1.12~1.13左右，已经很接近真实值。但是显然需要n很大时才能保证更大的精度。

下雨问题

如果一周每天都有50%的可能下雨，求连续三天下雨的概率？

这里假设每一天是否下雨是相互独立的，记 $Y_i=1_{\{the\quad i^{th} day \quad rains\}}$ ，则{Yi}独立。通过蒙特卡洛随机数的方法来模拟Yi的值：在[0,1]之间均匀分布任取一个数a，若a小于0.5则不下雨，否则下雨。之后再累加下雨天数，记录连续的下雨天数即可。对应的matlab程序如下：

n=1e4;ans=0;
for i=1:nc=0;y=0;p=0.5;for t=1:7r=unifrnd(0,1)if r<pc=c+1elsec=0endif c>=3y=1endendans=ans+y
end
ans=ans/n

n=1e4表示重复实验10000个星期，得出的结果为有3657个星期是有连续三天下雨的。

因此可以估算一周出现有连续三天下雨的概率约为36.57%；同理，一周出现连续三天晴天的概率也约为36.57%。而在代码中改变p的值，也可以得到其他的结果。下雨问题运用蒙特卡洛算法的精髓在于通过取随机数的方法来模拟某天是否下雨，无论是代码实现还是思路的理解都非常直观。

其他应用

经典的模型还有布冯投针问题等，本质上还是选取随机数进行近似地估计。

原理及误差分析

例如在下雨问题中，对于模拟参数Y，我们得到了一列结果Y1,Y2,Y3,……,Yn，由于它们独立同分布，由强大数定律知： $n\rightarrow \infty$ 时， $\frac{Y_1+\cdots+Y_n}{n}\rightarrow EY$ ，因此可用Y1,Y2,Y3,……,Yn的平均值来估计Y的真实值。同时由中心极限定理知，令 $S_n=Y_1+\cdots+Y_n$ ，那么观测平均值Sn/n与真实的均值μ=EY的差值为： $\frac{S_n}{n}-\mu=\frac{\sigma}{\sqrt{n}}$\frac{S_n-n\mu}{\sigma\sqrt{n}} $$ ,并且有 $\frac{S_n-n\mu}{\sigma\sqrt{n}} \rightarrow N(0,1)$ 。所以期望的观测平均值趋于0的变化速度为 $o(\frac{1}{\sqrt{n}})$ ,要使得精度增加两位小数点就需要一万次之多，误差相对来说比较大，精度不够高。但是有一些低精度就能够解决的问题使用蒙特卡洛算法是非常有效的。