当前位置:首页 >> 数学 >>

第五章蒙特卡罗方法(二)


蒙特卡罗方法(二)

5.2 蒙特卡罗积分
问题 物理学中常对一些自由度很大的系统感兴趣,例如 ? 多个电子组成的原子 ? 多个原子组成的气体或固体

对这种系统的描述常常可以归结为很高维数的积分。
例子 由 N 个原子组成的气体,这些原子在温度 1/β 下通过一个偶 对位势 v 相互作用,其经典配分函数为一 3N 维积分

通常积分方法面临的困难——计算量太大
我们很粗略的进行离散化,设每一维坐标只取 10 个不同的 值,于是被积函数必须在 103N 个点上求值。

对于一个并不大的 N 值,如 N=20,和一台每秒能求大约

1010个值的非常快速的计算机,这就需要大约 1050 秒的时
间——比宇宙的年龄的1032倍还多! 通常积分方法的适用范围 只有 N 非常小时,我们才可用前面讲过的求积公式。

蒙特卡罗积分的基本策略
从一维积分 说起

在区间 [a, b] 内随机的
撒 N 个点, 则 f(x) 在该区 间内的积分可写为

这是用蒙特卡罗方法计算定积分的基本策略

蒙特卡罗积分的误差分析
收敛性由大数定理来保证 设函数 f 在区间 [0,1] 上可积,我们在区间 [0,1] 上以均匀的 概率密度随机地取 N 个数 xi ,对每个 xi 计算出函数值 f(xi ) ,

大数法则告诉我们,这些函数值之和除以 N 的值将收敛于函
数 f 的期望值,也就是说下式是一个计算定积分的合理方法

收敛速度可从方差看出
1 2 1 ?1 2 ?I ? ? f ? ? N N? ?N

?
i ?1

N

?1 fi ? ? ?N
2

?
i ?1

N

? fi ? ?

2

? ? ? ?

其中 σf 2 是 f 的方差, 即 f 在积分区间上偏离其平均值的

程度的量度。

此式揭示了Monte Carlo 方法求积方法的两个重要的方面 1. 积分估值中的不确定度 σI 按 N-1/2 减小,所以用的点 越多,我们就会得到精确的答案,但是其收敛非常慢。

2. σf 越小(也就是说,若 f 尽可能平滑),精度越高。

考虑如下三种情况:
1. 一个极端情况是 f 一个常数,这时我们只需一点上的值就 可定出它的平均值。 2. 当 f 很平坦时,只需不多的点就能保证精度。 3. 另一极端情况是, f 除了在 x 的定义域有一个很窄的峰之 外其余都是零。若 xi 以等概率位于a 和 b 之间的任何地方, 那么除了少数几个之外它们绝大多数都将位于 f 的尖锋之外,

而只有这少数的几个的 fi 才不是零。这将导致对积分的一个
很差的估值。

例子 用Monte Carlo 方法计算积分
W(x)=1 N

dx ? ? ?0 1 ? x2 4 ? 0.78540
1

W(x)=(4-2x) /3

I
10 20 50 100 200 500 1000 2000 5000 0. 81491 0.73535 0.79606 0.79513 0.78677 0.78242 0.78809 0.78790 0.78963

?I
0.04638 0.03392 0.02259 0.01632 0.01108 0.00719 0.00508 0.00363 0.00227

I
0.79982 0.79071 0.78472 0.78838 0.78529 0.78428 0.78524 0.78648 0.78530

?I
0.00418 0.00392 0.00258 0.00194 0.00140 0.00091 0.00064 0.00045 0.00028

与传统积分方法做比较
求解一维定积分时,蒙特卡罗方法的效率其实是很差的。 标准的梯形积分方法的效率达到 N-2,而蒙特卡罗方法的效 率只有 N-1/2。 蒙特卡洛罗方法的误差不依赖于维数。梯形法是依赖于维数 的,总误差为 虽然在低维积分是不如常规积分方法,一般的,当维数大于 4 维时,用蒙特卡罗方法更好。 在维数很高的情况,只有蒙特卡洛罗方法可以用。

有待解决的问题
函数在积分区间内起伏很大时,积分误差也相应的很大。
物理中重要的积分,被积函数偏差都很大,例如统计力学中

的配分函数。
波尔兹曼因子在相空间的绝大多数区域都近似为 0。

解决办法——重要性抽样
Monte Carlo求积中的误差正比于被积函数的标准偏差。如果

可以设计出一种方案减小被积函数的标准偏差, 则Monte
Carlo 方法的效率也会得到改进。

设想有一个归一化的正的权函数 w (x)

积分

可以写成

如果现在我们作一个变量替换,从 x 变到



dy ? w ? x ? ; y ? x ? 0 ? ? 0; y ? x ? 1? ? 1 dx

积分变为

Monte Carlo求积方法:即在 y 坐标的区间 [0,1] 上取一组
均匀分布的随机样点,对相应的 f / w 的值求平均:

如果我们挑选一个 w 使它的行为和 f 相近,那么就可以使新的 被积函数 f / w 很平滑,从而减小 Monte Carlo估值的方差。 问题的关键在于能不能找到一个合适的 w,并且能够从y=w(x) 反推出函数 x(y) 。

dx ? ? ? 0.78540 ,我们选择归一化的权函数 对前例 ?0 2 1? x 4
1

其在积分区间上是正定的, 单调下降( f 也是)。由于在 x=0 和

x=1 两个端点上都有 f / w = 3 / 4, w 对 f 的行为近似得很好。
这样,新的积分变量为
1 y ? x ?4 ? x? 3

其反函数为

x ? 2 ? ?4 ? 3y?

1/2

计算精度相对前面 w(x)=1 的情况有明显的改进。

可以证明样点在 y 坐标的区间 [0,1] 内均匀分布隐含着样点在

x 坐标分布的概率密度是 dy/dx=w(x) 。
这意味着样点集中在最“重要”的 x 值附近,此处 w 大(因 而 f 也有希望大,对积分的贡献也大)。 在那些 w 和 f 值小的“不重要的”x 值附近,只有很少的 样点。 在对积分贡献大的区域多撒点,其余地方少撒点。 这种抽样方法称之为重要性抽样

重要性抽样小结
重要性抽样的两种理解方式:

1. 做积分变量的变换

新的积分变量 y 在区间 [0, 1] 均匀随机取样 2. 积分变量x 按照概率函数 w(x) 在 [0, 1] 随机取样

特定分布的随机变量的产生

Monte Carlo求积包括两种基本运算:产生以规定的概率密 度 w(x) (它可以是 1 )在积分区域内随机分布的坐标点, 然后在这些坐标点上计算函数 f / w 的值。

均匀分布的随机数是容易实现的,如何通过均匀分布的随机
数得到特定分布的随机数?

指数分布随机数的生成

我们已经知道如何生成 [0, 1] 之间均匀分布的随机数 y,现在
想寻找一个变换 x=x(y),将均匀随机数 y 变成想要的服从概 率分布 p 的随机数

因为


为简单起见 ,令 得

高斯分布随机数的生成 物理学中另一个有用的分布是高斯分布

用一个均匀分布的随机数 随机数

和一个指数分布 得到一个二维高斯分布随机



变换到极坐标

分布可写为

进一步简化 其中 生成[0,∞ ]范围内指数分布的随机数u和[0, 2π]范围内均匀分 布的随机数 θ ,然后通过下面的变换得到高斯分布随机数

任意分布的随机数的生成 概率分布函数w(x) 较简单时,可以通过解析的积分变换将 均匀分布的随机数转化为按照 w(x) 分布的随机数。 w(x) 较复杂时,只能通过数值方法来获得相应分布的随机数序

列。
设y( j) 为 [0, 1] 区间的一系列等间隔的值 y( j) ≡ j / M, j=0,

1, …, M。 y( j) 和 x ( j)通过下式相联系

现在问题当然就在于如何产生 x ( j) ,这可以通过把微分方 程 dy/dx= w(x) 的积分简单的离散化而得出。



我们得到方便的递推关系

这个关系可以同启动值 x(0) =0 一起使用。
从集合 0, 1, …, M 中以等概率选取一个整数j,对应的 x ( j) 就

是一个按w(x) 分布的随机数。

von Neumann 舍选法 舍弃

设我们对产生在 0 和 1 之间按 w(x)

分布的 x 有兴趣,并令 w '(x)是一个
正函数,并且在积分区域上 w '(x) > 接受

w(x)。

O

X→

1

选择两个随机变数 xi 和 η ,前者按 w’分布,而后者在 0

和1 之间均匀分布, 若 η 小于比值 w(xi )/ w’(xi ) ,就接受这
个 xi 值,否则就舍弃这个 xi 值。重复这一过程,就能得到 按照 w 分布的随机数。


赞助商链接
相关文章:
更多相关标签: