蒙特卡罗算法

31
4.必要时，还应改进模型以降低估计方差和减少试验费 用，提高模拟计算的效率
11
回顾几种连续型分布
1.均匀分布U(a,b) 其概率密度函数为
1 ,a x b f ( x) b a 0 , 其他
有
d c P{c x d } , 其中(c, d ) (a, b) ba
Monte Carlo Simulation Methods 蒙特卡罗模拟方法
1
主要内容
M-C方法概述
随机数的生成 模拟训练
2
一.M-C方法概述
蒙特卡洛(Monte Carlo)方法，或称计算机 随机模拟方法，是一种基于“随机数”的计 算方法。这一方法源于美国在第二次世界大 战中研制原子弹的“曼哈顿计划”。该计划 的主持人之一、数学家冯· 诺伊曼用驰名世 界的赌城—摩纳哥的Monte Carlo—来命名 这种方法，为它蒙上了一层神秘色彩。
xi 1 6xi mod11, ui+1 xi 1 /11 （a 6, m 11 ）
如果令 a 3, x0 1 ,得到序列: 1,3,9,5, 4,1,3,9........ 如果令 a 3, x0 2, 得到序列： 2, 6, 7,10,8, 2, 6.......
一般地，

A
f ( x)d (measure of A) (average of f over n random points in A)
9
Monte Carlo数值积分的优点
与一般的数值积分方法比较，Monte Carlo方法 具有以下优点：
1. 一般的数值方法很难推广到高维积分的情形，而 Monte Carlo方法很容易推广到高维情形
一般形式： xi 1 (axi c) mod m ui 1 xi 1 / m
1. c是非负整数.通过适当选取参数c可以改善 随机数的统计性质(独立性,均匀性).
2. 线性同余器可以达到的最长周期为 m 1 ，我们 可以通过适当的选择 m 和 a ，使无论选取怎样的 初值 x0 都可以达到最大周期(一般选取 m 为质数)
21
常用的线性同余生成器
Modulus m 2^31-1 =2147483647 Multiplier a 16807 39373 742938285 950706376 1226874159 40692 40014 Reference Lewis, Goodman, and Miller L’Ecuyer Fishman and Moore Fishman and Moore Fishman and Moore L’Ecuyer L’Ecuyer 22
针长 投针次数 相交次数
π的估计值
Wolf
1850
0.80
5000
2532
3.15956
Smith
1855
0.60
3204
1218
3.15665
Fox
1884
0.75
1030
489
3.15951
Lazzarini
1925
0.83
3408
1808
3.14159292
7
数值积分问题
选取(0, 1)中随机数序列x1， x2， x3， …… xn。则
24
从U(0,1)到其它概率分布的随机数
1.离散型随机数的模拟
2.连续型随机数的模拟
3.正态随机数的模拟
25
1.离散型随机数的模拟
设随机变量 X 的分布律为 令
P{X xi } pi
n i 1
(i 1,2,)
(n 1, 2 ,)
P(0) 0, P(n) pi ,
每个 xi j 有 m 种选择，所以向量 ( xi 1 , xi 2 ......xi -k ) 可以取 mk -1 个不同的值，所以这样的随机数生成器 的最大周期可以达到 m 生成器的周期。
k -1
1，大大提高了简单同余
23
算法实现
许多程序语言中都自带生成随机数的方法，如 c 中的 random() 函数，Matlab中的rand()函数等。 但这些生成器生成的随机数效果很不一样，比如 c 中的函数生成的随机数性质就比较差，如果用 c ，最好自己再编一个程序。Matlab 中的 rand() 函数，经过了很多优化。可以产生性质很好的随 机数，可以直接利用。
据此,可得产生 X 的随机数的具体过程为:每产 生一个(0,1)区间上均匀分布随机数U ,若
P(n 1) U P(n)
则令 X 取值 xn .
27
例1:
离散型随机变量X有如下分布律: X 0 1 2 P(x) 0.3 0.3 0.4 设 U1 ,U 2 ,,U N 是(0,1)上均匀分布的随机数,令
我们可产生一系列随机数 1, 2 , 3 , 4 , 5 , 6 ,....... 可简单取3个随机数构成一个随机点，即
(1, 2 , 3 ), (4 , 5 , 6 ),.......
相应地，
1 b 1 n a f ( x)dx n f ( xi ) ba i 1 d b 1 1 n c a f ( x, y)dxdy n f ( xi , yi ) (d c)(b a) i 1
5
于是有： l p P( X sin ) 2 0

l sin 2

0
2 2l dxd a a
2l ap
若我们独立重复地作 n 次投针试验，记 n ( A) 为 A 发生的次数。 fn ( A) 为 A
ˆ 在 n 次中出现的频率。假如我们取 fn ( A) 作为 p P( A) 的估计，即 p f n ( A) 。
3
基本思想很早以前就被人们所发现和 利用。17世纪，人们就知道用事件发 生的“频率”来决定事件的“概率”。 19世纪人们用投针试验的方法来决定π。 高速计算机的出现，使得用数学方法 在计算机上大量模拟这样的试验成为 可能。
4
从Buffon(蒲丰)投针问题谈起
Buffon 投针问题:平面上画很多平行线，间距为 a．向此平面投掷长为 l ( l < a) 的 针,此针与任一平行线相交的概率 p。
15
3. 指数分布 指数分布随机变量X的概率密度为
e x , x 0 f ( x) 0, x 0
指数分布常用来描述寿命问题.
16
二.随机数的生成
1.蒙特卡罗模拟的关键是生成优良的随机数。 2.在计算机实现中,我们是通过确定性的算法生成 随机数,所以这样生成的序列在本质上不是随机 的,只是很好的模仿了随机数的性质(如可以通过 统计检验)。我们通常称之为伪随机数(pseudorandom numbers)。 3.在模拟中，我们需要产生各种概率分布的随机数， 而大多数概率分布的随机数产生均基于均匀分布 U(0,1)的随机数。
12
均匀分布随机变量X的取值具有”均匀性” 均匀性特点:均匀分布随机变量X落在(a,b)内 任意子区间的概率只与子区间的长度有关,而与 子区间的位置无关. 可假设有这种特性的随机变量服从均匀分布.
13
2. 正态分布 正态分布随机变量X的概率密度函数是
1 1 x 2 f ( x) exp[ ( ) ], x R 2 2

1
0
1 n f ( x)dx E[ f ( xi )] f ( xi ) n i 1
1 误差约 ，它并不能和一些高级的数值积分算法比拟， n
但对多维情况，MC方法却很有吸引力。
8

1 1 1
0 0 0
1 n f ( x, y, z )dxdydz f ( xi , yi , zi ) n i 1
2147483399 2147483563
复杂一些的生成器
Multiple recursive generator
xi (a1 xi 1 a2 xi 2 ...... ak xi k ) mod m ui xi / m 需要选取种子(xk 1 , xk 2 .......x0 )
正态分布由两个参数 和 唯一确定.其中 是X 的均值(数学期望): =E(X),它确定了概率曲线的 中心位置,而 是X的标准差: D(X ) ,它确定了 概率曲线的”宽窄”程度.
14
在许多实际问题中,有一类随机变量可以表 示成为许多相互独立的随机变量之和,而其中每 个随机变量对总和只起微小的影响,这类随机变 量往往服从或近似服从正态分布.在实际应用中, 如果我们分析到一个随机变量受到较多独立的 微小因素的叠加影响,就可以用正态分布来模拟 这个变量.如:工厂产品的测量尺寸,农作物的收 获量,某地区成年人的身高,体重等可看成服从正 态分布的随机变量.
29
Matlab程序
Function r=rnd-beta(lamda) %模拟指数分布 %lamda表示指数分布的参数 r=-log(rand)/lamda; return
30
Matlab程序
Function y= rnd(mean, segema) %模拟均值为mean,方差为segma的正态分布 r=rand(1,12); x=sum(r)-6; y=segma*x+mean; return
ˆ 然后取 2l a.s. ˆ 作为 的估计。根据大数定律，当 n 时， p fn ( A) p. af n ( A) 2l P 。这样可以用随机试验的方法求得 的估计。历史上 af n ( A)
ˆ 从而有
有如下的试验结果。
6
试验者
时间(年)
17
乘同余法：
U(0,1)随机数的生成
xi 1 axi mod m ui 1 xi 1 / m
其中 xi , a, m 均为整数, x0 可以任意选取。
x0 称为种子,a 是乘因子,m是模数
18
一个简单的例子
当 x0 1 时,得到序列: 1,6,3,7,9,10,5,8,4,2,1,6,3......
0, 0 U i 0.3 xi 1, 0.3 U i 0.6 2, 0.6 U i
则 x1 , x2 ,, x N 是具有X分布律的随机数.
28
Matlab程序
Function r=rnd-u(a,b) %产生在[a,b]间均匀分布的随机数 r=a+(b-a)*rand; return
2. 一般的数值积分方法的收敛阶为 O(n 2 / d )，而 由中心极限定理可以保证 Monte Carlo 方法的 收敛阶为 O(n
1/ 2
) 。此收敛阶与维数无关，且在
10
高维时明显优于一般的数值方法。
M-C的基本思路
1.针对实际问题建立一个简单且便于实现的概率统计模 型，使所求的量（或解）恰好是该模型某个指标的概率 分布或者数字特征。 2.对模型中的随机变量建立抽样方法，在计算机上进行 模拟测试，抽取足够多的随机数，对有关事件进行统计 3.对模拟试验结果加以分析，给出所求解的估计及其精 度(方差)的估计
随机投针可以理解成针的中心 点与最近的平行线的距离X是均匀 地分布在区间 [0, a / 2] 上的r.v.，针 与平行线的夹角是均匀地分布 在区间 [0, ] 上的r.v.，且X与相互独立，
l 于是针与平行线相交的充要条件为 X 2 sin , l 即相交 A { : X 2 sin }.
将 {P(n)} 作为区间(0,1)的分点.若随机变量 U ~ U (0,1) ,有
P{P(n 1) U P(n)} P(n) P(n 1) pn , (n 1,2, )
26
令 {P(n 1) U P(n)} {X xn }
则有
P{X xn } pn
19
一个简单的例子(续)
上面的例子中，第一个随机数生成器的周期 长度是 10，而后两个生成器的周期长度只有 它的一半。我们自然希望生成器的周期越长 越好，这样我们得到的分布就更接近于真实 的均匀分布。
Hale Waihona Puke Baidu
在给定 m 的情况下，生成器的周期与 a 和 初值 x0 (种子)选择有关。
20
线性同余生成器(混合同余法) (Linear Congruential Generator )