Monte Carlo(蒙特卡洛方法)

蒙特卡罗模拟是一种随机模型．
1.抛100次硬币得到51个正面，并且接下来的10 次(即使不太可能刚巧10次)全为正面的情况是可 能出现的，这样，用110次的结果进行估计实际 上比用100次要差．
2.同样抛100次的结果也不近相同。
做任何的蒙特卡罗模拟，都要用到随机数．
要记住，对于根据模拟结果的预测寄予太多的 信任是有危险的，特别是在模拟中包含的假设没有 清楚表明的时候．还有，由于用了大量的数据和庞 大的计算，再加上非专业人员理解模拟模型和计算 机输出相对容易，所以常会导致对模拟结果的过分 相信．
一般形式： xi 1 (axi c) mod m ui 1 xi 1 / m
1. c是非负整数.通过适当选取参数c可以改善 随机数的统计性质(独立性,均匀性).
2. 线性同余器可以达到的最长周期为 m 1 ，我们 可以通过适当的选择 m 和 a ，使无论选取怎样的 初值 x0 都可以达到最大周期(一般选取 m 为质数)
蒙特卡洛(Monte Carlo)方法，或称计算机随机 模拟方法，是一种基于“随机数”的计算方法。 这一方法源于美国在第二次世界大战中研制原 子弹的“曼哈顿计划”。该计划的主持人之一、 数学家冯· 诺伊曼用驰名世界的赌城—摩纳哥的 Monte Carlo—来命名这种方法，为它蒙上了一层 神秘色彩。
P{P(n 1) U P(n)} P(n) P(n 1) pn , (n 1,2, )
令 {P(n 1) U P(n)} {X x } n
则有
P{X xn } pn
据此,可得产生 X的随机数的具体过程为:每产 生一个(0,1)区间上均匀分布随机数 U,若
例 2： 设 X ~ exp( ) 服从指数分布，则 X 的分布函数为： F ( x) 1 e x / , x 0 通过计算得 F ( y ) log(1 y ) ，则： X log(1 U )
1
服从指数分布(其中 U 服从均匀分布)
又因为 1－U 和 U 有着同样的分布，所以也可以取： X log(U )
1 确定行为的模拟
例：曲线下的面积
本节以曲线下的面 积为例说明蒙特卡罗 模拟在确定行为建模 中的应用．
下面的算法给出了用蒙特卡罗方法求曲线下面积 的计算机模拟的计算格式．
在给定区间上曲线y=cosx下面积的真值是2．注意到即使对 于产生的相当多的点数，误差也是可观的．对单变量函数，一般 说来，蒙特卡罗方法无法与在数值分析中学到的积分方法相比， 没有误差界以及难以求出函数的上界M也是它的缺点．然而，蒙 特卡罗方法可以推广到多变量函数，在那里它变得更加实用．
二.随机数的生成
1.蒙特卡罗模拟的关键是生成优良的随机数。 2.在计算机实现中,我们是通过确定性的算法生成 随机数,所以这样生成的序列在本质上不是随机 的,只是很好的模仿了随机数的性质(如可以通过 统计检验)。我们通常称之为伪随机数(pseudorandom numbers)。 3.在模拟中，我们需要产生各种概率分布的随机数， 而大多数概率分布的随机数产生均基于均匀分布 U(0,1)的随机数。
1. 生成 g ( x) 的样本 X ; 2. 生成U~U (0,1), 且U与X 独立； 3. 如果 U f ( X ) / cg ( X ) ，则取Y=X,返回步骤1， 否则舍去X,返回步骤1.
从Buffon(蒲丰)投针问题谈起
Buffon 投针问题:平面上画很多平行线，间距为 a．向此平面投掷长为 l ( l < a) 的 针,此针与任一平行线相交的概率 p。
随机投针可以理解成针的中心 点与最近的平行线的距离X是均匀 地分布在区间 [0, a / 2] 上的r.v.，针 与平行线的夹角是均匀地分布 在区间 [0, ] 上的r.v.，且X与相互独立，
一个简单的例子(续)
上面的例子中，第一个随机数生成器的周期 长度是 10，而后两个生成器的周期长度只有 它的一半。我们自然希望生成器的周期越长 越好，这样我们得到的分布就更接近于真实 的均匀分布。
在给定 m 的情况下，生成器的周期与 a 和 初值 x0 (种子)选择有关。
线性同余生成器(混合同余法) (Linear Congruential Generator )
(1) 生成 (0,1)上 均匀分布的随机数U。
(2) 计算 X F (U ) ,则 X 为来自 F ( x)
-1
分布的随机数.
例 1： 设 X ~ U (a, b) ，则其分布函数为 0 xa F ( x) b a 1, xa a xb xb
F -1 ( y ) a (b a) y ， 0 y 1 生成 U (0,1) 随机数 U，则 a (b - a)U 是来自 U (a, b) 的随机数。
P(n 1) U P(n)
则令 X取值
xn.
例1:
离散型随机变量X有如下分布律: X 0 1 2 P(x) 0.3 0.3 0.4 设 U1 ,U 2 ,,U 是 (0,1)上均匀分布的随机数,令 N
0, 0 U i 0.3 xi 1, 0.3 U i 0.6 2, 0.6 U i
所以在某些情况下，对对象的行为进行直接观测
或重复试验可能是不可行的。
在对象的行为不能做分析性的解释，或数据无法直 接收集的情况下，建模者可以用某种方式间接地模拟 其行为，试验所研究的供选择的各种方案，以估计它 们怎样影响对象的行为，然后收集数据来确定哪种方 案是最好的． 例如，为了得到一艘拟建造的潜艇受到的阻力， 造一个原型是不可行的，我们可以按比例建一个模型， 去模拟实际的潜艇的行为． 这里将研究另外一种形式的模拟——蒙特卡罗 (MonteCarlo)模拟，一般是借助于计算机完成的．
每个 xi j 有 m 种选择，所以向量 ( xi 1 , xi 2 ......xi -k ) 可以取 mk -1 个不同的值，所以这样的随机数生成器 的最大周期可以达到 m 生成器的周期。
k -1
1，大大提高了简单同余
算法实现
许多程序语言中都自带生成随机数的方法，如 c 中的 random() 函数，Matlab中的rand()函数等。 但这些生成器生成的随机数效果很不一样，比如 c 中的函数生成的随机数性质就比较差，如果用 c ，最好自己再编一个程序。Matlab 中的 rand() 函数，经过了很多优化。可以产生性质很好的随 机数，可以直接利用。
从U(0,1)到其它概率分布的随机数
1.离散型随机数的模拟
2.连续型随机数的模拟
3.正态随机数的模拟
1.离散型随机数的模拟
设随机变量 X 的分布律为 令
P{X xi } pi
n i 1
(i 1,2,)
(n 1, 2 ,)
P(0) 0, P(n) pi ,
Βιβλιοθήκη Baidu
将 {P(n)} 作为区间(0,1)的分点.若随机变量 U ~ U (0,1) ,有
ˆ 从而有
有如下的试验结果。
试验者
时间(年)
针长 投针次数 相交次数
π的估计值
Wolf
Smith
1850
1855
0.80
0.60
5000
3204
2532
1218
3.15956
3.15665
Fox
Lazzarini
1884
1925
0.75
0.83
1030
3408
489
1808
3.15951
3.14159292
则 x1 , x2 ,, x N 是具有X分布律的随机数.
2.连续型随机数的模拟
a.逆变换方法(常用) (Inverse Transform Method) b.舍取方法 (Acceptance-Rejection Method)
定理: 设随机变量Y的分布函数F(y)是连续函数, 而U是在(0,1)上均匀分布的随机变量, 令 X F 1 (U ) , 则Y与X有相同的分布.
证明： 由 F 1 (U ) 的定义和均匀分布的分布函数可得： P ( X x) P ( F 1 (U ) x) P (U F ( x )) F ( x )
由定理 1 ，要产生来自 F ( x) 的随机数，只要先 产生来自U (0,1) 随机数 u ，然后计算 F 1 (u ) 即 可。具体步骤如下：
能为那些要到达特定楼层的乘客提供最好的服务。 然而这种做法可能是难以接受的，因为在收集统计 数据时要再三惊扰乘客，并且电梯运行模式的不断变 化也会使乘客感到迷惑．
与此有关的另一个问题是大城市交通控制系统可 供选择的运行模式的检验，为了做试验而不停地改变 单行道的交通方向和配置交通信号将是不现实的．
舍取方法
舍取方法(Acceptance-Rejection )方法最早由 Von Neumann提出，现在已经广泛应用于各 种随机数的生成。 基本思路： 通过一个容易生成的概率分布 g 和一个取舍 准则生成另一个与 g 相近的概率分布 f 。
具体步骤：
假设 f ( x) 和 g ( x) 均为集合 上的概率密度函数。 且满足： f ( x) cg ( x) c 1, x 。易于从g(x)中 生成样本X.用如下步骤生成有Y:
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 )
l 于是针与平行线相交的充要条件为 X 2 sin , l 即相交 A { : X 2 sin }.
于是有： 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
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.......
常用的线性同余生成器
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
ˆ f n ( A) 。 在 n 次中出现的频率。假如我们取 fn ( A) 作为 p P( A) 的估计，即 p
ˆ 然后取 2l a.s. ˆ fn ( A) 作为 的估计。根据大数定律，当 n 时， p p. af n ( A) 2l P 。这样可以用随机试验的方法求得 的估计。历史上 af n ( A)
U(0,1)随机数的生成
乘同余法：
xi 1 axi
mod m
ui 1 xi 1 / m 其中 xi , a, m 均为整数, x0 可以任意选取。
x0称为种子,a 是乘因子,m是模数
一个简单的例子
当 x0 1 时,得到序列: 1,6,3,7,9,10,5,8,4,2,1,6,3......
Monte Carlo Simulation Methods (蒙特卡罗模拟方法)
一. 概述与思想 二. 随机数的生成. 三. 实例---港口模型 四. 作业
一 概述与思想
引例：电梯系统． 我们可以提出若干供选择的电梯运行模式，如设定 停偶数层、奇数层的电梯或直达电梯．理论上，对每种
供选择的模式都能够做若干次试验，以确定哪一种模式