蒙特卡洛方法与定积分计算
沙发~~
蒙特卡洛法是现代统计计算的基石,是贝叶斯获第二春的精髓所在,只是可惜,遇到高维向量,蒙特卡洛法就死翘翘了。
蒙特卡洛法是现代统计计算的基石,是贝叶斯获第二春的精髓所在,只是可惜,遇到高维向量,蒙特卡洛法就死翘翘了。
以上的程序是二维的,如果是三维的,比如这样一个积分,相当于求一个曲面y+z=1在x,y,z都在[0,1]间上积分,那么可以明显肉眼测得,这个积分值为0.5,如果用蒙特卡洛的思路,那么:
这样的J也是在0.5附近的,我模拟的结果是0.500552
其中的n1设成1000000,所以,我觉得有些情况,这种
方法如果可以降维的话,也可以解决一些高维度问题吧!
PS:自己乱想的,献丑了!
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 Jn1是自己输入的一个很大的数,
这样的J也是在0.5附近的,我模拟的结果是0.500552
其中的n1设成1000000,所以,我觉得有些情况,这种
方法如果可以降维的话,也可以解决一些高维度问题吧!
PS:自己乱想的,献丑了!
咦,我先试试看这次的内容可否提交,怎么我刚才提交的东西没有了呢?
有的,我看得到。
[未知用户] 刚才被误判为垃圾评论了。已经恢复回来。
其实你可以做两个简单的动画图片插进来,这样会让文章更生动。你写的这两种方法我都已经在animation包中实现了:
## 需要最新版本的animation包(>=1.0-10) library(animation) ani.options(interval = 0.2, nmax = 100) ## 随机投点 MC.hitormiss() ## 取随机数,算函数实现,算样本均值 MC.samplemean()
[未知用户] 高维的情况下马上就遇到维数灾难了,即使样本点非常多,在高维空间下也是很稀疏的,所以蒙特卡洛不会太精确。
第一个例子中精确值是用
pnorm(1) - pnorm(0)
算出来的么?如果是,可以注明一下,因为这个积分似乎没有显式表达式(正态分布的密度嘛)。[未知用户] 还没有学习animation包。求精确值用的代码是
integrate(f,0,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里面
特定的包更新啊?
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里面
特定的包更新啊?
[未知用户] 索得斯呢……那么其实这个精确值也是近似值了,只不过近似方法不一样而已了。
[未知用户]
pnorm(1)-pnorm(0)
和integrate(f,0,1)
结果一样都是0.3413447,都不是精确值。想更精确的话就用级数展开近似好了。高维下计算积分除了MC还没有别的太好的方法。一个例子就是EM算法,积分维数经常相当的高,没办法,现在都是用的MC-EM的做法,先抽,再叠代优化似然函数,再抽,再叠代优化似然函数。我叫这种方法——抽风(疯)流
[未知用户] pnorm应该是用一个很特殊的函数erf算得,具体的算法我也不知道。但是把一维积分转化成pnorm,把高维积分转化为高维正态分布函数是一个不错的思路,R里面的mvtnorm包里有pmvnorm函数。pnorm和pmvnorm都比别的近似方法快很多。如果无法转化,可以用Gauss-Hermit Quadrature或者拉普拉斯展开的方法求得近似的解,但是这两种方法应该都无法应付太高维的积分问题,估计2,3维顶天了
[未知用户] 所谓维数灾难一般也就是指高维空间下数据的稀疏性。这个很容易理解,比如有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],用蒙特卡洛积分算算模拟的值,并和真实值比较:
至于animation包的问题,你最好是装最新版的R。不然的话就在加载animation包之前运行update.packages(),否则包的目录被R进程占用了就不能再安装写入了。
如果用Windows,装R的时候我有个习惯,就是把安装过程中的安装目录的版本号去掉,直接装到C:\Program Files\R\目录下,而不带R-2.8.1之类的子目录。这样每次都安装到同一个目录下,卸载R的时候不会卸载以前安装的附加包;在重新安装R之后,只需要update.packages()就可以更新所有可更新的包了。
在蒙特卡洛积分这里,如果是高维积分,那么抽取随机数的时候就很难照顾到每个变量,也就是说每个变量根本抽不到几个数字样本量就已经很大了,这样由于随机性很差(试想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()就可以更新所有可更新的包了。
[未知用户] 倒也用不着这么“抠”啦,反正只要涉及计算机的东西几乎全都是近似……
[未知用户] 你说的erf是误差函数(又叫高斯误差函数),是个非基本函数,是用级数算的。
我再提供一份蒙特卡洛方法的学习材料吧:http://www.ceremade.dauphine.fr/~xian/coursMC.pdf 可算是讲义,很详尽。
[未知用户] 多谢小谢了,这个例子举得很好,我一看就懂了你的意思了。另外,
多维情况的那个例子中,有个命令set.seed(123),我发现如果
我不执行这句程序,也可以得到结果,请教了几个同学,最后我们认为
set.seed的作用是不是可以防止重复运行intProd时候,产生相同的随
机数?以及里面参数123是如何设定的啊?是不是需要遵从什么规则啊?
对了,我还想请教:我发现很多时候library和require使用上没什么差别,
看了help后也发现不了什么差别,难道真的没有差别?为什么R不用一个而是
存在两个一样的命令呢?
PS:我初学,问题多多。劳烦小谢了。
关于安装R包的问题,非常感谢小谢提供的意见,我就按着你说的做,太好
了,方便了很多,谢谢。
多维情况的那个例子中,有个命令set.seed(123),我发现如果
我不执行这句程序,也可以得到结果,请教了几个同学,最后我们认为
set.seed的作用是不是可以防止重复运行intProd时候,产生相同的随
机数?以及里面参数123是如何设定的啊?是不是需要遵从什么规则啊?
对了,我还想请教:我发现很多时候library和require使用上没什么差别,
看了help后也发现不了什么差别,难道真的没有差别?为什么R不用一个而是
存在两个一样的命令呢?
PS:我初学,问题多多。劳烦小谢了。
关于安装R包的问题,非常感谢小谢提供的意见,我就按着你说的做,太好
了,方便了很多,谢谢。