前几天在COS主站上看到一篇谈Monte Carlo求定积分的方法的文章。记得去年我学matlab时，求定积分便用到了这两种方法，而在上概率论课时，老师也略微涉及到相关的内容。当时的感觉时，定积分这种一般用确定型算法计算的东西，竟然还可以用随机的方法来实现，确是一种非常巧妙的思想，(Monte Carlo是摩纳哥的一个赌城，二战时Von Neumanm一伙人研究原子弹时弄出来的方法。）于是被震撼，也没有再深究其中更详细的内容。
但是昨天一看到这篇文章，便想起寒假时候看的那本《计算方法》中关于Monte Carlo平均值方法的论述，而最近在上计算方法的课，正纠结与一堆无聊的误差计算中，于是便再次翻开此书看看，竟发现了一些一直没太注意的很重要的章节。
具体如下：
- 平均值方法是根据弱大数定律证明得来的。构造如此
<bblatex>I(f)\approx 1/N\sum_{i=1}^{N}f(X_i)=I_N(f)</bblatex>,
均方误差
<bblatex>E|e_N|^2=E(I_N(f)-I(f))^2=1/N*Var(f(X))</bblatex>，
而由Schwartz不等式知
<bblatex>E|e_N|\leq sqrt(E|e_N|^2)=sqrt(Var(f)/N)</bblatex>,
从中便可知 MC 法具有半阶收敛性，且收敛阶不依赖问题的维数，这样在计算高维积分的时候，比较优越。但是收敛速度还是很慢的，是一个比较坏的收敛速度。因此在有其他高效率、高精度的确定型算法时，就不要采用MC法，那是在没有更好的算法时更好的抉择，比如超高维积分。
- 由均方误差可知，要提高 MC 法的精确度，那么久必须得增大N，减小VAR(f)，也就是说，要增加样本量，同时改善抽样方法的选择，这时便有很多方法选择，书中提到了四种方法，我就谈谈我的理解和试验。
- 分层抽样。核心就是将原来的[0,1]区间 M 等分后在各个小区间上分别再以均匀分布抽样 <bblatex>N/M</bblatex> 个点，然后以
<bblatex>!I_N=1/N \sum_{k=1}^{M}\sum_{i=1}^{N_k}f(X_i^{(k)}</bblatex>。
理论证明知这样方差减少，具体证明式子较多就不再写，否则就是一本教材的翻版了。
具体实验结果如下
<br />
M<-100<br />
N<-1000<br />
A=matrix(as.numeric(gl(M,N/M))-1,ncol=N/M,byrow=T)<br />
B=replicate(M,runif(N/M))<br />
x=1/M*t(B)+A/M<br />
y=exp(-x^2/2)/sqrt(2*pi)<br />
I=mean(y)<br />
可知只取了1000个样本点的情况下，精确度已经十分好了，而一般的均值模拟达到此精确度则需要<bblatex>10^7</bblatex>个点。
- 对偶变量法。
这是一种特殊的针对于定义域具有一定对称性的函数性质的技巧。
基于这样的命题:如果f(x)是单调的，则<bblatex>Cov(f(X),f(1-X))\leq0</bblatex>,这里<bblatex>X\sim \mathcal{u}[0,1]</bblatex>。则有
<bblatex>!I_N=1/(2*N)*sum_{i=1}^{N}(f(X_i)+f(1-X_i))</bblatex>,
则
<bblatex>EI_N=I(f)</bblatex>,
经过证明知 <bblatex>Var(I_N)\leq \frac{1}{2*N}Var(f)</bblatex>，方差显然减少。
但是这种方法只是对于定义域有一定对称性的函数，如果对一种函数先验了解不多，而且对称性不知，则运用会出现失误，所以具有一定局限性。模拟<bblatex>\int_{0}^{1}exp(-x^2/2)/sqrt(2*pi)dx</bblatex>的代码如下：
<br />
a=runif(10^6)<br />
y=exp(-a^2/2)/sqrt(2*pi)+exp(-(1-a)^2/2)/sqrt(2*pi)<br />
b=sum(y)/(2*10^6)<br />
- 重要性抽样：核心就是构造一个新的积分
<bblatex>!I(f)=\int_{0}^{1}f(x)dx=\int_{0}^{1}\frac{f(x)}{p(x)}*p(x)dx</bblatex>。
其中p(x)为[0,1]的概率密度，<bblatex>\int_{0}^{1}p(x)dx=1</bblatex>,<bblatex>p(x)\geq 0</bblatex>。这样得到类似于最初最简单的那种均值法的新的形式积分近似<bblatex>I(f)\approx 1/N *sum_{i=1}^{1}f(Y_i)/p(Y_i) </bblatex>，<bblatex>Y_i</bblatex>为遵循p(y)的独立同分布随机变量(i.i.d)
而之所以这样构造，可以根据这幅图来理解，![]()
</p>
当<bblatex>\int_{0}^{a}+\int_{b}^{1}f(x)dx</bblatex>占 I(f) 很少比例(即图像的边缘部分)，尽可能的不让随机点去计算这些值，这样的“拟合”积分较好（面积）,否则很多的随机点在[a,b]外，为计算小比例的东西做贡献，这样必定精度较差，所以我们要期望尽可能的值出现在峰值处。我们称其为“重要性”抽样。理论分析，
<bblatex>!Var_X(f)=\int_{0}^{1}(f-i(f))^{2}dx=\int_{0}^{1}f^{2}dx-I^{2}(f)</bblatex>,
<bblatex>!Var_Y(f/p)=\int_{0}^{1}(f/p)^{2}pdy-I^{2}(f)=\int_{0}^{1}f^2/pdy-I^{2}(f)</bblatex>.
只要选取适当的p(y)就可以使方差减少。
当然不可能是var(f/p）=0，这在图的模拟便可知，p(y)不可能与f(x)相合的同时，还有<bblatex>\int_{0}^{1}p(x)dx=1</bblatex>。但是p难以选择，需要对问题先验有了一定的了解程度才能适当选取。比如对于正态分布求<bblatex>\int_{0}^{1}exp(-x^2/2)/sqrt(2*pi)dx</bblatex>，我不太好选p，就没有做相关的模拟了。
- 控制变量法。
核心也是构造，数学式子<bblatex>I_N(f)=i/N\sum_{i=1}^{N}(f-g)(X_i)+I(g)</bblatex>.
需要I(g)已知，但是这是同重要性抽样一样，都不知道，也需要对问题先验有了一定的了解程度，才能好的选择，使I(g)已知,并且有<bblatex>Var(f-g)\leq Var(f)</bblatex>。相关模拟我也没有再做
总之，正如书中所言，减小误差都只是具有指导思想的意义，具体问题，往往需要我们对被积函数有先验的了解，才能够减小误差，提高精度。确实如此，我的例子都是比较特殊或者简单的，如果对于复杂的问题，而且对该函数比较怪异少见，那么在以上的方法中寻找什么p(x)，都是一项非常艰巨或者不可能完成的任务。因此，我下一结论，奇技淫巧之所以巧无外乎我们有了个千里眼和顺风耳，瞎碰瞎撞弄出一个巧妙的解法那是很难的，大部分都是前人解法突然激发灵感得来（此时此刻竟想起了贝叶斯来...)。所以对于专门解决高维积分强有力算法——Metropolis算法我是敬佩有加，这一算法2000年惊被评为20世纪十大算法之一。后面产生的什么模拟退火，拟Monte Carlo都与他有莫大的渊源。而伟大的东西终究充满了高深而深刻的东西，我目前楞死没有看太明白，已经被符号给堆死了。但是我想过段时间在查阅相关资料应该可以弄明白吧。
[attachment=200929,45]
