Monte Carlo方法之一:积分方法

新用户注册 | 用户登陆 | 刷新
论坛嘉宾: sage

权权


发表文章数: 20
内力值: 94/94
贡献度: 277
人气: 18

学术成员

Monte Carlo方法之一:积分方法 [文章类型: 原创]

先说题外话,回复一下轩轩的帖子,蒙特卡洛(Monte Carlo)是摩纳哥公国一个城镇,位于地中海沿岸,以其赌场和豪华酒店而闻名,所以就有了以随机方法应用于数值计算的一类方法,被称为Monte Carlo.

子曰:诗三百,一言以蔽之,曰:思无邪。

物理学总是试图把纷繁的现象纳入单纯统一的纲领,至少可以减少不必要的记忆量。于是就总是想说“一言以蔽之”,怎么样怎么样。。。但说这样的话是要小心的,大多数情况下是领悟还不够透彻,还未到达把书从厚读薄的境界,以至于忽略了这样那样隐含的假设。

不过,如果不过分追求严格,只是想记忆起来方便一些,那多说几句这样的话倒也无妨。所谓Monte Carlo方法,一言以蔽之,曰:算积分。

数值计算总是与离散化形影不离的,数值积分也是一样。一般的计算方法是把积分区域划分为均匀的网格(如同Riemann积分的定义一样),用离散格点上的求和来近似积分值。这种方法在计算高维积分时会遇到困难,如果需要计算一个N维积分,于是我们有了N个坐标轴,每个坐标轴都均分成格点,比如说10个格点,那么总的数据点的个数是10的N次方,依赖于N指数式增长;估计一下,比如我们要算一个60维的积分,N=60,假设我们的计算机每秒钟能够做10的7次方次运算,那总的运算时间是10的53次方秒,这大约等于10的39次方个哈勃时间。

而高维积分在物理中一点都不罕见,经典统计力学的配分函数是一个在3N维位形空间中的积分(粒子数为N),量子多体问题中的可观测量的平均值是由对多体波函数的积分得到的,而量子力学中的路径积分实质上也是由高维积分来定义的。

而在物理学发展的现阶段,计算高维积分唯一的普遍方法,是Monte Carlo.

假设我们需要计算如下形式的积分:

I = ∫ f(x)p(x)dx

其中,x代表高维空间中的一个点;函数p(x)具有以下性质,

p(x) > 0 , ∫ p(x)dx = 1

那么这个恒正而且归一的函数将被诠释为概率密度分布函数,而我们的积分可以被诠释为函数f(x)在概率分布p(x)下的平均值(或者说期望值),它通过对这一概率分布取样来计算:

I ≈ (1/N) Σ f(x_i) = <f>

这里N为数据样本点的个数,数据点序列{x_i}是按照概率分布p(x)来抽样的,求和"Σ"对数据点进行。


上述变换就是Monte Carlo积分的基本精神,因为需要用到随机抽样,必然伴随统计误差。因为x_i是按照概率密度p(x)分布的随机变量,f(x_i)也是随机变量,而中心极限定理告诉我们,一组独立随机变量之和的概率分布是高斯,其方差等于每一项随机变量的方差之和,如果记随机变量f(x)的方差为σ(f),那么

σ(I) = (1/N^2) Σ σ(f(x_i)) = (1/N^2)∙Nσ(f) = (1/N)σ(f)

标准差(s)是方差的平方根,所以从上式我们可以得出,用随机抽样来近似计算积分所引起的误差与数据点个数的平方根成反比:

s(I) = (1/√N)∙s(f)

这样的收敛速度算是比较慢的,与通常的数值积分方法比较一下就可以知道。

设积分维数为D,数据点总数为N,那么每一个坐标轴上均匀分布着N^(1/D)个数据点,也就是说格点之间的距离 h~N^(-1/D),比如说一般用梯形近似面积的积分方法(trapezoidal),那么每个格点引起的误差正比于

h^(D+2) ~ N^(-1)∙N^(-2/D)

而所有格点累计起来的误差为

s ~ N∙h^(D+2) ~ N^(-2/D)

比较一下Monte Carlo方法:s ~ 1/√N 可以发现一个有趣的事实,当积分区域的维数小于4时,用均匀网格计算积分的方法收敛速度更快,当维数大于4时,Monte Carlo方法的收敛速度更快,而当维数等于4时,两种方法不分高下,(插一句废话,4=3+1哦). 而计算高维积分时,Monte Carlo 方法是较优的选择。

不妨回顾一下布丰投针实验来结束本篇,在分布着等距平行木纹的地板上投针,要求针的长度小于木纹之间的距离,几何概形计算结果表明,针与一条木纹相交的概率可以用针的长度、木纹间距和圆周率π表示。而用几何概形计算概率实质上归结为面积的计算,也就是积分的计算,布丰投针实验可以说是用随机抽样计算积分的始祖。

如今在计算机上,可以很方便地进行简化版的部分投针实验,比若设一个边长为2的正方形,随机地在正方形区域内生成点,并记录落在正方形的内切圆中的点的个数与总点数的比例,就可以逼近圆周率。

发表时间: 2007-07-22, 23:55:40 个人资料

轩轩


发表文章数: 59
内力值: 126/126
贡献度: 465
人气: 118

学术成员

Re: Monte Carlo方法之一:积分方法 [文章类型: 原创]

假设我们需要计算如下形式的积分:

I = ∫ f(x)p(x)dx


假如p(x)不归一呢?_____________我没有仔细看全文:)


这种离散的方法在我印象中和离散的空间(loop量子引力)还可能扯上关系,我记得有一个人,用Monte Carlo方法计算我们空间的维度,在小尺度上是2.x. 反正不是一个整数. 这看上去似乎有分形有关系,也许自然是牛比的



还有一个是题外话了

格点是很有趣的.
1.椭圆曲线不能通过无穷多个格点
2.圆内格点计数,是一个很好的数论问题. 陈景润他们喜欢的问题_____当然这估计也可以用Monte Carlo方法计算.因为Monte Carlo方法想来是牛比的

i will love you till null infinity
<<相对论通俗演义>>

发表时间: 2007-07-23, 00:23:58 个人资料

权权


发表文章数: 20
内力值: 94/94
贡献度: 277
人气: 18

学术成员

Re: Monte Carlo方法之一:积分方法 [文章类型: 原创]

是这样的,任意一个积分都可以写成

I = ∫ f(x)p(x)dx

的形式,而且有无穷多种写法。比如要算某积分形如:

I = ∫ g(x)dx

只要找任意一个概率分布函数:

p(x) > 0 , ∫ p(x)dx = 1

然后把积分写成下面的形式:

I = ∫ [g(x)/p(x)]∙p(x)dx

然后令

f(x) = g(x)/p(x)

就可以了。问题是如何选取适当的概率分布p(x)才能是最有效率,一般的原则是要使得

f(x) = g(x)/p(x)

的形状越平坦越好,也就是说被积函数g(x)比较大的地方p(x)也比较大,反之亦然,这样很明显的结果就是f(x)的涨落会比较小。所以一般会选取一些解析形式简单又能大致模拟被积函数形状的概率分布p(x)来作上述变换。

发表时间: 2007-07-23, 00:35:14 个人资料

轩轩


发表文章数: 59
内力值: 126/126
贡献度: 465
人气: 118

学术成员

Re: Monte Carlo方法之一:积分方法 [文章类型: 原创]

问题是如何选取适当的概率分布p(x)才能是最有效率,


------------------这个看来还是比较人为的,比如,poisson分布,或者高斯分布?只有一个概率最高处,但被积分的函数的峰值可能有很多处:)



第2个问题

matlab mathematica全可以做积分 他们的原理是Monte carlo吗? 或者,它们之间的区别在哪里? 高维??

i will love you till null infinity
<<相对论通俗演义>>

发表时间: 2007-07-23, 00:59:38 个人资料

权权


发表文章数: 20
内力值: 94/94
贡献度: 277
人气: 18

学术成员

Re: Monte Carlo方法之一:积分方法 [文章类型: 原创]

概率分布函数也不是只有高斯、泊松那样的,也是可以取得奇形怪状的,只要满足归一性和恒正就可以了。

mathematica里算积分估计就是用传统的网格算法吧,当然会有各种各样的改进,因为用mathematica一般都是算低维积分,用Monte Carlo不合算的。

发表时间: 2007-07-23, 01:14:16 个人资料

轩轩


发表文章数: 59
内力值: 126/126
贡献度: 465
人气: 118

学术成员

Re: Monte Carlo方法之一:积分方法 [文章类型: 原创]

还有一个弱问题:)_____就直接问吧:)


Monte Carlo在目前有没有现成的软件?

还是说,要自己编写程序,然后再进行计算.

i will love you till null infinity
<<相对论通俗演义>>

发表时间: 2007-07-23, 07:44:15 个人资料

权权


发表文章数: 20
内力值: 94/94
贡献度: 277
人气: 18

学术成员

Re: Monte Carlo方法之一:积分方法 [文章类型: 原创]

有很多现成的code,但并没有做成用户界面很友好的软件。

发表时间: 2007-07-23, 17:35:06 个人资料
您尚未登陆 | 用户登陆