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

By 邓一硕 @ 2010/03/08
关键词：Monte-Carlo, 定积分, 模拟, 蒙特卡洛分类：统计计算
作者信息：来自中央财经大学;统计学专业。
版权声明：本文版权归原作者所有，未经许可不得转载。原文可能随时需要修改纰漏，全文复制转载会带来不必要的误导，若您想推荐给朋友阅读，敬请以负责的态度提供原文链接；点此查看如何在学术刊物中引用本文

常规引用方式
邓一硕. 蒙特卡洛方法与定积分计算. 统计之都, 2010.03. URL: http://cos.name/2010/03/monte-carlo-method-to-compute-integration/.
BibTeX引用

@ARTICLE{,
  AUTHOR = {邓一硕},
  TITLE = {蒙特卡洛方法与定积分计算},
  JOURNAL = {统计之都},
  YEAR = {2010},
  month = {03},
  URL = {http://cos.name/2010/03/monte-carlo-method-to-compute-integration/},
}

本文讲述一下蒙特卡洛模拟方法与定积分计算，首先从一个题目开始：设 $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

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

结论

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

# 程序写得稍微有点晦涩，大概意思： # 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

## 需要最新版本的animation包（>=1.0-10） library(animation) ani.options(interval = 0.2, nmax = 100) ## 随机投点 MC.hitormiss() ## 取随机数，算函数实现，算样本均值 MC.samplemean()

有边界区间上的核密度估计

一、一个例子
核密度估计应该是大家常用的一种非参数密度估计方法，从某种程度上来说它的性质比直方图更好，可以替代直方图来展示数据的密度分布。但是相信大家会经常遇到一个问题，那就是有些数据是严格大于或等…阅读全文 »

用R也能做精算—actuar包学习笔记（一）

本文是对R中精算学专用包actuar使用的一个简单教程。actuar项目开始于2005年，在2006年2月首次提供公开下载，其目的就是将一些常用的精算函数引入R系统。目前，提供的函数主要涉及风险理论，…阅读全文 »

分月存档

huang shuai on 2010/03/10 at 15:34

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

焦静 on 2010/03/10 at 19:06

以上的程序是二维的，如果是三维的，比如这样一个积分，相当于求一个曲面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：自己乱想的，献丑了！

谢益辉 on 2010/03/11 at 05:19

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

回复
- 焦静 on 2010/03/11 at 11:16
  
  哦，好像还是不是很了解，我再研究研究吧，高维用降维思路怎么就不可以了啊，哎，本人愚笨，要不小谢给我举个例子吧！另外，我最近才算是踏踏实实的准备从新扫盲学习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里面
  特定的包更新啊？
- 谢益辉 on 2010/03/11 at 12:48
  所谓维数灾难一般也就是指高维空间下数据的稀疏性。这个很容易理解，比如有10000个样本点均匀分布在样本空间里，如果数据的维数是2，那么正方形的每条边平均有100个样本点；如果数据维数是4，那么一个4维“正方体”每条边上平均只能分布10个数据点；如果维数更高，那么每一维都只能分布到更少的数据点。
  
  在蒙特卡洛积分这里，如果是高维积分，那么抽取随机数的时候就很难照顾到每个变量，也就是说每个变量根本抽不到几个数字样本量就已经很大了，这样由于随机性很差（试想rnorm(2)能算是正态分布的代表么？），积分的结果可能就不靠谱。不知道我这种解释是不是对的。
  
  用程序来说，计算一个n重积分 $\int_0^1\int_0^1\cdots\int_0^1x_1x_2\cdots x_ndx_1dx_2\cdots dx_n$ ，这个积分的结果我们知道，就是 $1/2^n$ ，用蒙特卡洛积分算算模拟的值，并和真实值比较：
```
# 程序写得稍微有点晦涩，大概意思：
# 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()就可以更新所有可更新的包了。
- gaotao on 2010/03/12 at 00:35
  
  冒个泡，同意谢老大的观点。因为随机投点法比较粗糙，而且用的随机数应该是伪随机的，遇到高维时，单纯的加大投点的随机数，意义不是很大，当然还因为在高维下数据太稀疏的缘故，所以随机投点法虽常用但是对于需要比较精确一点的解时误差还是挺大的。
  而对于平均值法，现在改进的方法倒是还蛮多的。在平均值法中的均方误差等于E|e_N|^2=Var(f)/N，Var(f)为f(X)的方差，所以改进时需要增大N和减小Var(f)，因此抽样方法的选择就是一个很重要的技巧，这里就不再过多探讨，可以参阅相关计算方法的书。用这个方法的话，根据上面的公式可知MC收敛阶不依赖维数（半阶收敛性），因此用较大型的计算机便可以得到较精确的解，而不用随机投点法了。
  但是当有更高效率和精度的确定性算法时，还是不会用MC法的，他是往往没有更好的算法时才采用，比如上千上万维的积分…
- 胡荣兴 on 2010/03/30 at 20:55
  
  很久之前用过蒙特卡罗方法，现在都生疏了。蒙特卡罗的好处就是计算简化，不用去考虑那么复杂的积分或微分甚至偏微分。
  而对维数不敏感也是其一个优点之一，在谢兄这倒成了其缺点了。有两位经济学家用该方法解作过1900多维计算。
- 谢益辉 on 2010/03/30 at 22:52
  
  那要看怎么理解“对维数不敏感”了，如果只是“不限制维数”的话，那这算不上优点，因为高维情况下需要的样本点会急剧增加，否则积分可能不精确（如前例）。当然，蒙特卡洛方法包括思想类似的一大类方法，也许存在某些方法能克服维数问题，我不了解。

焦静 on 2010/03/10 at 19:08

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

邓一硕 on 2010/03/10 at 23:24

有的，我看得到。

谢益辉 on 2010/03/10 at 23:52

刚才被误判为垃圾评论了。已经恢复回来。

回复
- 焦静 on 2010/03/11 at 21:03
  
  多谢小谢了，这个例子举得很好，我一看就懂了你的意思了。另外，
  多维情况的那个例子中，有个命令set.seed(123)，我发现如果
  我不执行这句程序，也可以得到结果，请教了几个同学，最后我们认为
  set.seed的作用是不是可以防止重复运行intProd时候，产生相同的随
  机数？以及里面参数123是如何设定的啊？是不是需要遵从什么规则啊？
  
  对了，我还想请教：我发现很多时候library和require使用上没什么差别，
  看了help后也发现不了什么差别，难道真的没有差别？为什么R不用一个而是
  存在两个一样的命令呢？
  PS：我初学，问题多多。劳烦小谢了。
  关于安装R包的问题，非常感谢小谢提供的意见，我就按着你说的做，太好
  了，方便了很多，谢谢。
- 邱怡轩 on 2010/03/14 at 09:07
  
  set.seed()不是为了防止产生相同的随机数，相反，正是为了产生相同的随机数，因为这样一来，每次运行的结果就是一样的，这使得结果可重复被验证。参数的选择很随意，比如选个幸运数字什么的。
  library()和require()在作用上没什么大差别，但require()会返回一个逻辑值，加载成功就返回TRUE，否则为FALSE。

谢益辉 on 2010/03/11 at 05:12

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

谢益辉 on 2010/03/11 at 06:05

第一个例子中精确值是用pnorm(1) - pnorm(0)算出来的么？如果是，可以注明一下，因为这个积分似乎没有显式表达式（正态分布的密度嘛）。

邓一硕 on 2010/03/11 at 10:56

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

回复
- 谢益辉 on 2010/03/11 at 11:30
  
  索得斯呢……那么其实这个精确值也是近似值了，只不过近似方法不一样而已了。
- 邓一硕 on 2010/03/11 at 11:39
  
  pnorm(1)-pnorm(0)和integrate(f,0,1)结果一样都是0.3413447，都不是精确值。想更精确的话就用级数展开近似好了。
- 谢益辉 on 2010/03/11 at 12:54
  
  倒也用不着这么“抠”啦，反正只要涉及计算机的东西几乎全都是近似……
colinisstudent on 2010/03/11 at 12:43

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

回复
- 邓一硕 on 2010/03/11 at 13:18
  
  你说的erf是误差函数（又叫高斯误差函数），是个非基本函数，是用级数算的。

colinisstudent on 2010/03/11 at 12:30

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

谢益辉 on 2010/03/11 at 13:42

我再提供一份蒙特卡洛方法的学习材料吧：http://www.ceremade.dauphine.fr/~xian/coursMC.pdf 可算是讲义，很详尽。

焦静 on 2010/03/12 at 17:27

大致看了看，非常好的一份材料，不仅研究了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

回复

魏太云 on 2010/03/12 at 14:47

还是期待量子计算机的横空出世吧，到时候可不惧怕这些大计算量问题了。

colinisstudent on 2010/03/14 at 03:41

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

回复

提高Monte Carlo模拟定积分计算的精度 - Just Do It! on 2010/03/14 at 02:07

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

huang shuai on 2010/04/25 at 15:55

http://www.stat.columbia.edu/~cook/movabletype/mlm/RMHMC_MG_BC_SC_REV_08_04_10.pdf

刚在GELMAN博客上看到这篇文章，非常兴奋，这个东东将是蒙特卡洛计算的一个里程碑！

P.S. 看到“有两位经济学家用该方法解作过1900多维计算。” 请问可有文献啊。这个维数有点吓人，如果他们成功了，那简直太值得一学了。。。有文献就拿出来共享嘛。。

贡献我的力量

最新评论

随机文章

标签云

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

相关文章

28 Responses to “ 蒙特卡洛方法与定积分计算 ”

Leave a Reply

搜索

推荐阅读

有边界区间上的核密度估计

相关文章

用R也能做精算—actuar包学习笔记（一）

相关文章

本月热评

中文论坛新帖

维基百科动态

英文论坛新帖

分月存档