蒙特卡洛方法与定积分计算

本文讲述一下蒙特卡洛模拟方法与定积分计算,首先从一个题目开始:设$0\leq f(x) \leq 1$,用蒙特卡洛模拟法求定积分$J=\int_{0}^{1}f(x)dx$的值。

随机投点法

设$(X,Y)$服从正方形 $\{0\leq x \leq 1,0\leq y\leq 1\}$上的均匀分布,则可知 $X,Y$分别服从[0,1]上的均匀分布,且$X,Y$相互独立。记事件$A=\{Y\leq f(X)\}$,则$A$的概率为

$
P(A)=P(Y\leq f(X))=\int_{0}^{1}\int_{0}^{f(x)}dydx=\int_{0}^{1}f(x)dx=J
$

即定积分$J$的值 就是事件$A$出现的频率。同时,由伯努利大数定律,我们可以用重复试验中$A$出现的频率作为 $p$的估计值。即将$(X,Y)$看成是正方形$\{0\leq x \leq 1,0\leq y \leq 1\}$内的随机投点,用随机点落在区域${y\leq f(x)}$中的频率作为定积分的近似值。这种方法就叫随机投点法,具体做法如下:

图1 随机投点法示意图

1、首先产生服从 $[0,1]$上的均匀分布的$2n$个随机数( $n $为随机投点个数,可以取很大,如 $n=10^4$ )并将其配对。

2、对这$n$对数据 $(x_i,y_i),i=1,2,…,n$ ,记录满足不等式$y_i\leq f(x_i)$的个数,这就是事件 $A$ 发生的频数$\mu_n$,由此可得事件$A$发生的频率$\frac{\mu_n}{n}$,则$J\approx \frac{\mu_n} {n}$。

举一实例,譬如要计算$\int_{0}^{1}e^{-x^2/2}/\sqrt{2\pi}dx$,模拟次数$n=10^4$时,R代码如下:

n=10^4;
 x=runif(n);
 y=runif(n);
 f=function(x)
 {
 exp(-x^2/2)/sqrt(2*pi)
 }
 mu_n=sum(y<f(x));
 J=mu_n/n;
 J

模拟次数 $n=10^5$ 时,令$n=10^5$,其余不变。

定积分$\int_{0}^{1}e^{-x^2/2}/\sqrt{2\pi}dx$的精确值和模拟值如下:

表1

精确值 $10^3$ $10^4$ $10^5$ $10^6$ $10^7$
0.3413447 0.342 0.344 0.34187 0.341539 0.341302
注:精确值用integrate(f,0,1)求得

扩展

如果你很细心,你会发现这个方法目前只适用于积分区间$[0,1]$ ,且积分函数 $f(x)$ 在区间$[0,1]$上的取值也位于 $[0,1]$ 内的情况。那么,对于一般区间 $[a,b]$ 上的定积分$J’=\int_{a}^{b}g(x)dx$ 呢?一个很明显的思路,如果我们可以将 $J’$ 与$J$ 建立代数关系就可以了。

首先,做线性变换,令 $y=(x-a)/(b-a)$ ,此时,

$x=(b-a)y+a$, $J’=(b-a)\int_{0}^{1}g[(b-a)y+a]dy$。

进一步如果在区间$[a,b]$上有$c\leq g(x) \leq d $ ,令

$f(y)=\frac{1}{d-c}{g(x)-c}=\frac{1}{d-c}{g[a+(b-a)y]-c} $,

则$0\leq f(y) \leq 1$。此时,可以得到$J’=\int_{a}^{b}g(x)dx=S_0J+c(b-a) $。

其中,$S_0=(b-a)(d-c) $, $J=\int_{0}^{1}f(y)dy $。

这说明,用随机投点法计算定积分方法具有普遍意义。

举一个实例,求定积分 $J’=\int_{2}^{5}e^{-x^2/2}/\sqrt{2\pi}dx $。

显然$a=2,b=5$,$c=min\{g(x)|2\leq x \leq 5\},d=max\{g(x)|2\leq x \leq 5\}$,由于$g(x)=e^{-x^2/2}/\sqrt{2\pi}$在 $[2,5]$上时单调减函数,所以$c=g(5),d=g(2)$,$S_0=(b-a)(d-c)$。R中代码为

a=2;
 b=5;
 g=function(x)
 {
 exp(-x^2/2)/sqrt(2*pi);
 }
 f=function(y)
 {
 (g(a+(b-a)*y)-c)/(d-c);
 }
 c=g(5);d=g(2);s_0=(b-a)*(d-c);
 n=10^4;
 x=runif(n);y=runif(n);
 mu_n=sum(y<=f(x));
 J=mu_n/n;
 J_0=s_0*J+c*(b-a);

定积分 $J’=\int_{2}^{5}e^{-x^2/2}/\sqrt{2\pi}dx$ 的精确值和模拟值如下:

表2

真实值 $10^3$ $10^4$ $10^5$ $10^6$ $10^7$
0.02274985 0.02332792 0.02311736 0.02262659 0.02284152 0.02278524
注:精确值用integrate(g,2,5)求得)

平均值法

蒙特卡洛模拟法计算定积分时还有另一种方法,叫平均值法。这个原理也很简单:设随机变量 $X$ 服从$[0,1]$上的均匀分布,则$Y=f(X)$的数学期望为

$E(f(X))=\int_{0}^{1}f(x)dx=J $

所以估计$J$的值就是估计$f(X)$的数学期望值。由辛钦大数定律,可以用$f(X)$的观察值的均值取估计$f(X)$的数学期望。具体做法:

先用计算机产生$n$个服从$[0,1]$上均匀分布的随机数:$x_i,i=1,2,…,n$。

对每一个$x_i$ ,计算$f(x_i)$。

计算$J\approx \frac{1}{n}\sum_{i=1}^{n}f(x_i)$。

譬如,计算 $J=\int_{0}^{1}e^{-x^2/2}/\sqrt{2\pi}dx $,R中的代码为

n=10^4;
 x=runif(n);
 f=function(x)
 {
 exp(-x^2/2)/sqrt(2*pi)
 }
 J=mean(f(x));

其精确值和模拟值如下:
表3

真实值 $10^3$ $10^4$ $10^5$ $10^6$ $10^7$
0.3413447 0.3405831 0.3410739 0.3414443 0.3414066 0.3413366

平均值法与随机投点法类似可以进行扩展,这里不再赘述。

结论

用蒙特卡洛模拟法计算定积分具有普遍意义。蒙特卡洛模拟方法为我们提供了一个看待世界的新思路,即一个不具随机性的事件可以通过一定的方法用随机事件来模拟或逼近。

关于邓一硕

专注于R语言在金融投资分析和计量经济学中的应用;《R in a Nut Shell》、《R Graphics Cookbook》、《Introductory Statistics with R》等书的译者。

蒙特卡洛方法与定积分计算》有38个想法

  1. 沙发~~
    蒙特卡洛法是现代统计计算的基石,是贝叶斯获第二春的精髓所在,只是可惜,遇到高维向量,蒙特卡洛法就死翘翘了。

  2. 以上的程序是二维的,如果是三维的,比如这样一个积分,相当于求一个曲面y+z=1在x,y,z都在[0,1]间上积分,那么可以明显肉眼测得,这个积分值为0.5,如果用蒙特卡洛的思路,那么:

    n1 < - scan()
    n <- as.numeric(n1)
    x <- runif(n)
    y <- runif(n)
    z <- runif(n)
    z1 <- function(y) {
        1 - y
    }
    num <- sum(z < z1(y))
    J <- num/n
    J 
    

    n1是自己输入的一个很大的数,
    这样的J也是在0.5附近的,我模拟的结果是0.500552
    其中的n1设成1000000,所以,我觉得有些情况,这种
    方法如果可以降维的话,也可以解决一些高维度问题吧!
    PS:自己乱想的,献丑了!

    1. 高维的情况下马上就遇到维数灾难了,即使样本点非常多,在高维空间下也是很稀疏的,所以蒙特卡洛不会太精确。

      1. 哦,好像还是不是很了解,我再研究研究吧,高维用降维思路怎么就不可以了啊,哎,本人愚笨,要不小谢给我举个例子吧!另外,我最近才算是踏踏实实的准备从新扫盲学习R,但执行你上面的程序时出现了问题,百思不得其解啊!特来请教这个低级问题:我用2.8.1版本的R去加载animation版本1.0-4后,我就用了一个命令update.packages(),这个命令不是能将所有安装包都更新么,但是后面继续执行:
        library(animation)
        ani.options(interval = 0.2, nmax = 100)
        MC.hitormiss()
        后,仍然报错为:错误: 没有”MC.hitormiss”这个函数
        我直接用update.packages(“animation”)的命令也总是报错为:

        Warning message:
        In list.files(lib) : list.files: ‘animation’目录不可读

        而我同学直接安装后,就可以直接执行以上命令,他的R是2.10.1的,
        所以我们推测可以是我的R版本太低的缘故,是这样么?如果是,那么
        我还想请教小谢,怎么可以在不改变R版本的情况下,将library里面
        特定的包更新啊?

      2. 所谓维数灾难一般也就是指高维空间下数据的稀疏性。这个很容易理解,比如有10000个样本点均匀分布在样本空间里,如果数据的维数是2,那么正方形的每条边平均有100个样本点;如果数据维数是4,那么一个4维“正方体”每条边上平均只能分布10个数据点;如果维数更高,那么每一维都只能分布到更少的数据点。

        在蒙特卡洛积分这里,如果是高维积分,那么抽取随机数的时候就很难照顾到每个变量,也就是说每个变量根本抽不到几个数字样本量就已经很大了,这样由于随机性很差(试想rnorm(2)能算是正态分布的代表么?),积分的结果可能就不靠谱。不知道我这种解释是不是对的。

        用程序来说,计算一个n重积分[latex]\int_0^1\int_0^1\cdots\int_0^1x_1x_2\cdots x_ndx_1dx_2\cdots dx_n[/latex],这个积分的结果我们知道,就是[latex]1/2^n[/latex],用蒙特卡洛积分算算模拟的值,并和真实值比较:

        # 程序写得稍微有点晦涩,大概意思:
        # 10个均匀分布随机数相乘(x1x2...xn),看看比1个均匀分布
        # (y)小的概率,考虑到数值精度问题,把乘法改成了exp(log(加法))
        # ndim:数据维数;nrep:抽样次数
        intProd = function(ndim, nrep = 1e+05) {
            c(MC.value = mean(runif(nrep, 0, 1) < = exp(replicate(nrep, 
                sum(log(runif(ndim)))))), TRUE.value = 1/2^ndim)
        }
        set.seed(123)
        res1 = replicate(10, intProd(15))
        res2 = replicate(10, intProd(1))
        # 模拟值/真实值
        > res1[1, ]/res1[2, ]
         [1] 1.96608 0.65536 0.98304 1.31072 1.63840
         [6] 0.65536 1.63840 1.63840 1.96608 0.32768
        > res2[1, ]/res2[2, ]
         [1] 1.00626 1.00260 1.00078 0.99894 0.99892
         [6] 1.00616 1.00214 0.99872 1.00050 1.00618
        

        可以看出,样本量固定为十万,当维数为15的时候,模拟结果和真实结果的倍数波动相对比较大,维数为1的时候,结果波动相对较小。

        至于animation包的问题,你最好是装最新版的R。不然的话就在加载animation包之前运行update.packages(),否则包的目录被R进程占用了就不能再安装写入了。

        如果用Windows,装R的时候我有个习惯,就是把安装过程中的安装目录的版本号去掉,直接装到C:\Program Files\R\目录下,而不带R-2.8.1之类的子目录。这样每次都安装到同一个目录下,卸载R的时候不会卸载以前安装的附加包;在重新安装R之后,只需要update.packages()就可以更新所有可更新的包了。

      3. 冒个泡,同意谢老大的观点。因为随机投点法比较粗糙,而且用的随机数应该是伪随机的,遇到高维时,单纯的加大投点的随机数,意义不是很大,当然还因为在高维下数据太稀疏的缘故,所以随机投点法虽常用但是对于需要比较精确一点的解时误差还是挺大的。
        而对于平均值法,现在改进的方法倒是还蛮多的。在平均值法中的均方误差等于E|e_N|^2=Var(f)/N,Var(f)为f(X)的方差,所以改进时需要增大N和减小Var(f),因此抽样方法的选择就是一个很重要的技巧,这里就不再过多探讨,可以参阅相关计算方法的书。用这个方法的话,根据上面的公式可知MC收敛阶不依赖维数(半阶收敛性),因此用较大型的计算机便可以得到较精确的解,而不用随机投点法了。
        但是当有更高效率和精度的确定性算法时,还是不会用MC法的,他是往往没有更好的算法时才采用,比如上千上万维的积分…

      4. 很久之前用过蒙特卡罗方法,现在都生疏了。蒙特卡罗的好处就是计算简化,不用去考虑那么复杂的积分或微分甚至偏微分。
        而对维数不敏感也是其一个优点之一,在谢兄这倒成了其缺点了。有两位经济学家用该方法解作过1900多维计算。

      5. 那要看怎么理解“对维数不敏感”了,如果只是“不限制维数”的话,那这算不上优点,因为高维情况下需要的样本点会急剧增加,否则积分可能不精确(如前例)。当然,蒙特卡洛方法包括思想类似的一大类方法,也许存在某些方法能克服维数问题,我不了解。

  3. 咦,我先试试看这次的内容可否提交,怎么我刚才提交的东西没有了呢?

      1. 多谢小谢了,这个例子举得很好,我一看就懂了你的意思了。另外,
        多维情况的那个例子中,有个命令set.seed(123),我发现如果
        我不执行这句程序,也可以得到结果,请教了几个同学,最后我们认为
        set.seed的作用是不是可以防止重复运行intProd时候,产生相同的随
        机数?以及里面参数123是如何设定的啊?是不是需要遵从什么规则啊?

        对了,我还想请教:我发现很多时候library和require使用上没什么差别,
        看了help后也发现不了什么差别,难道真的没有差别?为什么R不用一个而是
        存在两个一样的命令呢?
        PS:我初学,问题多多。劳烦小谢了。
        关于安装R包的问题,非常感谢小谢提供的意见,我就按着你说的做,太好
        了,方便了很多,谢谢。

      2. set.seed()不是为了防止产生相同的随机数,相反,正是为了产生相同的随机数,因为这样一来,每次运行的结果就是一样的,这使得结果可重复被验证。参数的选择很随意,比如选个幸运数字什么的。
        library()require()在作用上没什么大差别,但require()会返回一个逻辑值,加载成功就返回TRUE,否则为FALSE。

  4. 其实你可以做两个简单的动画图片插进来,这样会让文章更生动。你写的这两种方法我都已经在animation包中实现了:

    ## 需要最新版本的animation包(>=1.0-10)
    library(animation)
    ani.options(interval = 0.2, nmax = 100)
    
    ## 随机投点
    MC.hitormiss()
    
    ## 取随机数,算函数实现,算样本均值
    MC.samplemean()
    
  5. 第一个例子中精确值是用pnorm(1) - pnorm(0)算出来的么?如果是,可以注明一下,因为这个积分似乎没有显式表达式(正态分布的密度嘛)。

    1. 还没有学习animation包。求精确值用的代码是integrate(f,0,1)

      1. 索得斯呢……那么其实这个精确值也是近似值了,只不过近似方法不一样而已了。

      2. pnorm(1)-pnorm(0)integrate(f,0,1)结果一样都是0.3413447,都不是精确值。想更精确的话就用级数展开近似好了。

      3. 倒也用不着这么“抠”啦,反正只要涉及计算机的东西几乎全都是近似……

    2. pnorm应该是用一个很特殊的函数erf算得,具体的算法我也不知道。但是把一维积分转化成pnorm,把高维积分转化为高维正态分布函数是一个不错的思路,R里面的mvtnorm包里有pmvnorm函数。pnorm和pmvnorm都比别的近似方法快很多。如果无法转化,可以用Gauss-Hermit Quadrature或者拉普拉斯展开的方法求得近似的解,但是这两种方法应该都无法应付太高维的积分问题,估计2,3维顶天了

      1. 你说的erf是误差函数(又叫高斯误差函数),是个非基本函数,是用级数算的。

  6. 高维下计算积分除了MC还没有别的太好的方法。一个例子就是EM算法,积分维数经常相当的高,没办法,现在都是用的MC-EM的做法,先抽,再叠代优化似然函数,再抽,再叠代优化似然函数。我叫这种方法——抽风(疯)流

    1. 大致看了看,非常好的一份材料,不仅研究了Monte Carlo,还研究了Bayesian,非常值得大家看哦,不过缺点就在于这份PDF没有目录,所以我为这个材料做了一份目录,希望可以帮助大家学习:

      Monte Carlo Statistical Methods
      by Christian P.Robert

      context

      1 introduction
      1.1 Statistical Models 5
      1.2 Likelihood Methods 10
      1.3 Bayesian Methods 21
      1.4 Deterministic Numerical Methods 28
      1.5 simulation versus numerical analysis:
      when is it useful? 31

      2 Random Variable Generation 35
      2.1 Basic Methods 37
      2.2 Beyond uniform distributions 51

      3 Monte Carlo Intergration 76
      3.1 introduction 77
      3.2 Classical Monte Carlo Integration 80
      3.3 importance sampling 87
      3.4 Acceleration Methods 98

      4 Markov Chains 108
      4.1 Basic Notions 110
      4.2 Irreducibility 115
      4.3 Transience/Recurrence 123
      4.4 Invariant Measures 126
      4.5 Ergodicity and stationarity 130
      4.6 Limit Theorems 134

      5 Monte Carlo Optimization 139
      5.1 Introduction 140
      5.2 Stochastic Exploration 142
      5.3 Stochastic Approximation 173
      5.3.3 MCEM 195
      6 The Metropolis-Hastings Algorithm
      6.1 Markov Chain Monte Carlo 197
      6.2 The Metropolis-Hastings Algorithm 199
      6.3 A Collection of Metropolis-Hastings Algorithms 204
      6.4 Extensions 217

      7 The Gibbs Sampler 231
      7.1 General Principles 232
      7.1.5 Hierarchical models 253
      7.2 Data Augmentation 255
      7.3 Improper Priors 271

      8 Diagnosing Convergence 278
      8.1 Stopping the Chain 279
      8.2 Monitoring Stationarity Convergence 282
      8.3 Monitoring Average Convergence 290

      9 Implementation in Missing Data Models 317
      9.1 First examples 319
      9.2 Finite mixtures of distributions 340
      9.3 Extensions 354

    1. 我倒是觉得并不是所有的积分问题都值得硬上的,一些问题能够被简化,使得积分问题在相对比较低维的情况下进行。

  7. 非常喜欢MCMC,就是不太懂,希望大侠们多提供一些基础性的材料和讲义。

  8. 如果能加点图例就更好了,像hit and miss方法,看图一目了然。再就是可以考虑再加上Importance sampling, MCMC等方法,感觉两种方法还是略显单调。最后说下版面问题:1.最上面因为页面右面有标签看起来还挺舒服,但是往下走走,感觉右边空白太大了(莫非是为了非宽屏电脑考虑?) 再就是留言评论呈阶梯状排列使得后面发言的人评论区很小,像谢老大写的代码被挤到一个小狭长的空间里了,都要托滚动条才能看完整,不知到能否改进(能不能让右边距也右移呢?)

    1. 1、这个建议倒是提醒了我如何充分利用页面空间。往下看评论的时候页面右边确实显得有点空,此时不妨放一个浮动层,里面写点重要通知之类的,既能引起人的注意,又能改善空间布局。
      2、目前没办法,这种情况应该也不会经常出现。你可以双击代码,然后复制到别的地方看。

  9. 任何方法,在处理高维问题是,都会显得捉襟见肘。但是相比普通的数值积分,蒙特卡罗积分的收敛速度和维度无关。这是因为数值积分相当于在积分局域上均匀抽样,而蒙特卡罗积分则在“重要的区域”抽取多的样本,而不重要的区域抽取少的样本,这样的比较下,高下立判。当然了,现在的主要问题是如何从高维分布里抽样,这才是真正的难题!

  10. 大哥,【进一步如果在区间[a,b]上有c≤g(x)≤d ,令f(y)=1d−cg(x)−c=1d−cg[a+(b−a)y]−c,则0≤f(y)≤1。此时,可以得到J′=∫bag(x)dx=S0J+c(b−a)。】这里少了括号好吗!

W.C Zhu进行回复 取消回复

邮箱地址不会被公开。 必填项已用*标注