前几天在$COS$主站上看到一篇谈Monte Carlo求定积分的方法的文章。记得去年我学$matlab$时，求定积分便用到了这两种方法，而在上概率论课时，老师也略微涉及到相关的内容。当时的感觉时，定积分这种一般用确定型算法计算的东西，竟然还可以用随机的方法来实现，确是一种非常巧妙的思想，($Monte, Carlo$是摩纳哥的一个赌城，二战时$Von, Neumanm$一伙人研究原子弹时弄出来的方法。）于是被震撼，也没有再深究其中更详细的内容。

但是昨天一看到这篇文章，便想起寒假时候看的那本《计算方法》中关于$Monte, Carlo$平均值方法的论述，而最近在上计算方法的课，正纠结与一堆无聊的误差计算中，于是便再次翻开此书看看，竟发现了一些一直没太注意的很重要的章节。

具体如下：

• 平均值方法是根据弱大数定律证明得来的。构造如此

  $$I(f)\approx 1/N\sum_{i=1}^{N}f(X_i)=I_N(f)$$

  
  均方误差

  $$E|e_N|^2=E(I_N(f)-I(f))^2=\frac{1}{N}Var(f(X))$$

  
  而由$Schwartz$不等式知

  $$E|e_N|\leq \sqrt{E|e_N|^2}=\sqrt{Var(f)/N}$$

  ,从中便可知 $MC$ 法具有半阶收敛性，且收敛阶不依赖问题的维数，这样在计算高维积分的时候，比较优越。但是收敛速度还是很慢的，是一个比较坏的收敛速度。因此在有其他高效率、高精度的确定型算法时，就不要采用$MC$法，那是在没有更好的算法时更好的抉择，比如超高维积分。
• 由均方误差可知，要提高 $MC$ 法的精确度，那么久必须得增大$N$，减小$VAR(f)$，也就是说，要增加样本量，同时改善抽样方法的选择，这时便有很多方法选择，书中提到了四种方法，我就谈谈我的理解和试验。
  1. 分层抽样。核心就是将原来的 $[0,1]$ 区间 $M$ 等分后在各个小区间上分别再以均匀分布抽样 $N/M$ 个点，然后以
     

     $$I_N=1/N \sum_{k=1}^{M}\sum_{i=1}^{N_k}f(X_i^{(k)})$$

     
     理论证明知这样方差减少，具体证明式子较多就不再写，否则就是一本教材的翻版了。
     模拟$\int_{0}^{1}exp(-x^2/2)/\sqrt{2pi}dx$,具体实验结果如下 
     ^?View Code RUBY

     1 2 3 4 5 6 7

     M<-100 N<-1000 A=matrix(as.numeric(gl(M,N/M))-1,ncol=N/M,byrow=T) B=replicate(M,runif(N/M)) x=1/M*t(B)+A/M y=exp(-x^2/2)/sqrt(2*pi) I=mean(y)

     结果为0.3413287，可知只取了1000个样本点的情况下，精确度已经十分好了，而一般的均值模拟
     达到此精确度则需要$10^7$个点。

  2. 对偶变量法。
     这是一种特殊的针对于定义域具有一定对称性的函数性质的技巧。基于这样的命题:如果$f(x)$是单调的，则$Cov(f(X),f(1-X))\leq0$,这里$X\sim U[0,1]$。则有
     

     $$I_N=1/(2N)\sum_{i=1}^{N}(f(X_i)+f(1-X_i))$$

     ,
     则$EI_N=I(f)$,
     经过证明知 $Var(I_N)\leq \frac{1}{2N}Var(f)$，方差显然减少。 

     但是这种方法只是对于定义域有一定对称性的函数，如果对一种函数先验了解不多，而且对称性不知，则运用会出现失误，所以具有一定局限性。
     模拟$\int_{0}^{1}exp(-x^2/2)/\sqrt{2pi}dx$的代码如下：

     ^?View Code RUBY

     1 2 3

     a=runif(10^6) y=exp(-a^2/2)/sqrt(2*pi)+exp(-(1-a)^2/2)/sqrt(2*pi) b=sum(y)/(2*10^6)

     结果为 0.3413487，精确度也得到了提高。

  3. 重要性抽样：核心就是构造一个新的积分

     $$I(f)=\int_{0}^{1}f(x)dx=\int_{0}^{1}\frac{f(x)}{p(x)}p(x)dx$$

     其中$p(x)$为$[0,1]$的概率密度，$\int_{0}^{1}p(x)dx=1$,$p(x)\geq 0$。这样得到类似于最初最简单的那种均值法的新的形式积分近似$I(f)\approx 1/N\sum_{i=1}^{1}f(Y_i)/p(Y_i)$，$Y_i$为遵循p(y)的独立同分布随机变量(i.i.d),而之所以这样构造，可以根据这幅图来理解，当$\int_{0}^{a}+\int_{b}^{1}f(x)dx$占$I(f)$ 很少比例(即图像的边缘部分)，尽可能的不让随机点去计算这些值，这样的“拟合”积分较好（面积）,否则很多的随机点在 $[a,b]$ 外，为计算小比例的东西做贡献，这样必定精度较差，所以我们要期望尽可能的值出现在峰值处。我们称其为“重要性”抽样。理论分析，
     

     $$Var_X(f)=\int_{0}^{1}(f-i(f))^{2}dx=\int_{0}^{1}f^{2}dx-I^{2}(f)$$

     ,
     

     $$Var_Y(f/p)=\int_{0}^{1}(f/p)^{2}pdy-I^{2}(f)=\int_{0}^{1}f^2/pdy-I^{2}(f)$$

     .
     只要选取适当的 $p(y)$ 就可以使方差减少。
     当然不可能是$Var(f/p)=0$，这在图的模拟便可知，$p(y)$不可能与$f(x)$相合的同时，还有$\int_{0}^{1}p(x)dx=1$。但是$p$难以选择，需要对问题先验有了一定的了解程度才能适当选取。比如对于正态分布求$\int_{0}^{1}exp(-x^2/2)/\sqrt(2pi)dx$，我不太好选$p$，就没有做相关的模拟了。
  4. 控制变量法。
     核心也是构造，数学式子$I_N(f)=i/N\sum_{i=1}^{N}(f-g)(X_i)+I(g)$.
     需要$I(g)$已知，但是这是同重要性抽样一样，都不知道，也需要对问题先验有了一定的了解程度，才能好的选择，使$I(g)$已知,并且有$Var(f-g)\leq Var(f)$。相关模拟我也没有再做

总之，正如书中所言，减小误差都只是具有指导思想的意义，具体问题，往往需要我们对被积函数有先验的了解，才能够减小误差，提高精度。确实如此，我的例子都是比较特殊或者简单的，如果对于复杂的问题，而且对该函数比较怪异少见，那么在以上的方法中寻找什么$p(x)$，都是一项非常艰巨或者不可能完成的任务。因此，我下一结论，奇技淫巧之所以巧无外乎我们有了个千里眼和顺风耳，瞎碰瞎撞弄出一个巧妙的解法那是很难的，大部分都是前人解法突然激发灵感得来（此时此刻竟想起了贝叶斯来...)。所以对于专门解决高维积分强有力算法——Metropolis算法我是敬佩有加，这一算法2000年惊被评为20世纪十大算法之一。后面产生的什么模拟退火，拟$Monte Carlo$都与他有莫大的渊源。而伟大的东西终究充满了高深而深刻的东西，我目前楞死没有看太明白，已经被符号给堆死了。但是我想过段时间在查阅相关资料应该可以弄明白吧。

最后再扯句从中学到了的随机向量的构造——正态分布随机变量——Box-Muller法
以前学的是用n个均匀分布随机数由中心极限定律构造，但是那个要得到较好的精度需要N很大，故需要精确的还是用上面的那个方法更比较有保证，咱看公式也心安嘛。