蒙特卡洛方法与定积分计算
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
蒙特卡洛方法与定积分计算
By 邓一硕 @ 2010/03/08
关键词:Monte-Carlo, 定积分, 模拟, 蒙特卡洛分类:统计计算
作者信息:来自中央财经大学;统计学专业。
版权声明:本文版权归原作者所有,未经许可不得转载。原文可能随时需要修改纰漏,全文复制转载会带来不必要的误导,若您想推荐给朋友阅读,敬请以负责的态度提供原文链接;点此查看如何在学术刊物中引用本文
本文讲述一下蒙特卡洛模拟方法与定积分计算,首先从一个题目开始:设,用蒙特卡洛模拟法求定积分的值。
随机投点法
设服从正方形上的均匀分布,则可知分别服从[0,1]上的均匀分布,且相互独立。记事件,则的概率为
即定积分的值就是事件出现的频率。同时,由伯努利大数定律,我们可以用重复试验中出现的频率作为的估计值。即将看成是正方形
内的随机投点,用随机点落在区域中的频率作为定积分的近似值。这种方法就叫随机投点法,具体做法如下:
图1 随机投点法示意图
1、首先产生服从上的均匀分布的个随机数(为随机投点个数,可以取很大,如)并将其配对。
2、对这对数据,记录满足不等式的个数,这就是事件发生的频数,由此可得事件发生的频率,则。
举一实例,譬如要计算,模拟次数时,R代码如下:n=10^4;
x=runif(n);
y=runif(n);
f=function(x)
{
exp(-x^2/2)/sqrt(2*pi)
}
mu_n=sum(y J=mu_n/n; J 模拟次数时,令,其余不变。 定积分的精确值和模拟值如下: 表1 0.3413447 0.342 0.344 0.34187 0.341539 0.341302 注:精确值用integrate(f,0,1)求得 扩展 如果你很细心,你会发现这个方法目前只适用于积分区间,且积分函数在区间上的取值也位于内的情况。那么,对于一般区间上 的定积分呢?一个很明显的思路,如果我们可以将与建立代数关系就可以了。 首先,做线性变换,令,此时, ,。 进一步如果在区间上有,令 , 则。此时,可以得到。 其中,,。 这说明,用随机投点法计算定积分方法具有普遍意义。 举一个实例,求定积分。 显然,,由于 在上时单调减函数,所以, 。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); 定积分的精确值和模拟值如下: 表2 0.02274985 0.02332792 0.02311736 0.02262659 0.02284152 0.02278524 注:精确值用integrate(g,2,5)求得) 平均值法 蒙特卡洛模拟法计算定积分时还有另一种方法,叫平均值法。这个原理也很简单:设随机变量服从上的均匀分布,则的数学期望为 所以估计的值就是估计的数学期望值。由辛钦大数定律,可以用 的观察值的均值取估计的数学期望。具体做法: 先用计算机产生个服从上均匀分布的随机数:。 对每一个,计算。 计算。 譬如,计算,R中的代码为 n=10^4; x=runif(n); f=function(x) { exp(-x^2/2)/sqrt(2*pi) } J=mean(f(x)); 其精确值和模拟值如下: 表3 0.3413447 0.3405831 0.3410739 0.3414443 0.3414066 0.3413366 平均值法与随机投点法类似可以进行扩展,这里不再赘述。 结论 用蒙特卡洛模拟法计算定积分具有普遍意义。蒙特卡洛模拟方法为我们提供了一个看待世界的新思路,即一个不具随机性的事件可以通过一定的方法用随机事件来模拟或逼近。