蒙特卡罗模拟方法

19
Inverse Transform Method
设随机变量 X 的分布函数为 F ( x) ,定义: F 1 ( y ) = inf{x : F ( x) ≥ y}, 0 ≤ y ≤1 即F 1 ( y )是F ( x)的反函数. 定理 1 :设随机变量 U 服从 (0,1) 上的均匀分布,
则 X = F 1 (U ) 的分布函数为 F ( x) .
24
标准正态分布随机数的生成
正态分布是概率统计中最重要的分布,在此 我们着重讨论如何生成标准正态分布随机数. 引理: 设 U1 , U 2 是独立同分布的 U (0,1) 变量,令:
Z1 = (2 ln U1 )
1/ 2
cos(2π U 2 )
Z 2 = (2 ln U1 )1/ 2 sin(2π U 2 ) 则 Z1 与 Z 2 独立,且均服从标准正态分布. 证明:由随机向量的变换公式易得,略.
5
Monte Carlo数值积分的优点
与一般的数值积分方法比较,Monte Carlo方法 具有以下优点:
1. 一般的数值方法很难推广到高维积分的情形,而 Monte Carlo方法很容易推广到高维情形 Carlo方法很容易推广到高维情形
2. 一般的数值积分方法的收敛阶为 O(n 收敛阶为 O(n
16
由rand()函数生成的U[0,1]随机数
17
由rand函数生成的2维随机点
18
从U(0,1)到其它概率分布的随机数
U(0,1)的均匀分布的随机数,是生成其他概率 分布随机数的基础,下面我们主要介绍两种将 U(0,1)随机数转换为其他分布的随机数的方法.
1.逆变换方法 (Inverse Transform Method) 2.舍取方法 (Acceptance-Rejection Method)
25
Box-Muller 算法
1. 生成两个独立的均匀分布随机数 U1 , U 2
2.令 Z1 = (2 ln U1 )
1/ 2 1/ 2
cos(2π U 2 ) sin(2π U 2 )
Z 2 = (2 ln U1 )
则 Z1,Z 2为独立的正态随机数
26
逆变换方法(一)
我们无法通过具体的数学表达式计算正态分布函数 的逆函数,我们必须通过数值的方法逼近正态函数 下面我们介绍 Beasley-Springer-Moro 方法.
2147483399 2147483563
复杂一些的生成器(一)
1.Combining Generators:
考虑 J 个简单线性同余生成器: x j ,i +1 = a j x j ,i mod m j , u j ,i +1 = x j ,i +1 / m j j = 1, 2....J xi +1 = ∑ (1)
(1) 生成 (0,1)上 均匀分布的随机数U.
1
(2) 计算 X = F -1 (U ) ,则 X 为来自 F ( x) 分布的随机数.
21
几个具体例子(一)
例 1 : X ~ U (a, b) ,则其分布函数为 设 0 xa F ( x) = b a 1,
-1
x<a a≤ x<b x≥b
随机投针可以理解成针的中心 点与最近的平行线的距离X是均匀 地分布在区间 [0, a / 2] 上的r.v.,针 与平行线的夹角Φ是均匀地分布 在区间 [0, π ] 上的r.v.,且X与Φ相互独立,
l 于是针与平行线相交的充要条件为 X ≤ 2 sin Φ, l 即相交 A = {ω : X ≤ 2 sin Φ}.
1/ 2
2 / d
),而
由中心极限定理可以保证 Monte Carlo 方法的 ) .此收敛阶与维数无关,且在
6
高维时明显优于一般的数值方法.
随机模拟计算的基本思路
1.针对实际问题建立一个简单且便于实现的概率统计模 型,使所求的量(或解)恰好是该模型某个指标的概率 分布或者数字特征. 2.对模型中的随机变量建立抽样方法,在计算机上进行 模拟测试,抽取足够多的随机数,对有关事件进行统计 3.对模拟试验结果加以分析,给出所求解的估计及其精 度(方差)的估计 4.必要时,还应改进模型以降低估计方差和减少试验费 用,提高模拟计算的效率
j =1 J j 1
x j ,i +1 mod(m1 1) xi +1 > 0
ui +1 = xi +1 / m1
(m1 1) / m1 xi +1 = 0 (其中 m1 为所有 m j 中最大的一个)
14
复杂一些的生成器(二)
2.Multiple recursive generator
xi = (a1 xi 1 + a2 xi 2 + ...... + ak xi k ) mod m ui = xi / m 需要选取种子(xk 1 , xk 2 .......x0 )
7
随机数的生成
1.蒙特卡罗模拟的关键是生成优良的随机数. 2.在计算机实现中,我们是通过确定性的算法生成 随机数,所以这样生成的序列在本质上不是随机 的,只是很好的模仿了随机数的性质(如可以通过 统计检验).我们通常称之为伪随机数(pseudorandom numbers). 3.在模拟中,我们需要产生各种概率分布的随机数, 而大多数概率分布的随机数产生均基于均匀分布 U(0,1)的随机数.
证明: 由 F (U ) 的定义和均匀分布的分布函数可得: P( X ≤ x) = P( F 1 (U ) ≤ x) = P (U ≤ F ( x)) = F ( x)
20
1
Inverse Transform Method
由定理 1 ,要产生来自 F ( x) 的随机数,只要先 产生来自U (0,1) 随机数 u ,然后计算 F (u ) 即 可.具体步骤如下:
12
常用的线性同余生成器
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 13
10
一个简单的例子(续)
上面的例子中,第一个随机数生成器的周期 长度是 10,而后两个生成器的周期长度只有 它的一半.我们自然希望生成器的周期越长 越好,这样我们得到的分布就更接近于真实 的均匀分布.
在给定 m 的情况下,生成器的周期与 a 和 初值 x0 (种子)选择有关.
11
线性同余生成器 (Linear Congruential Generator )
时间(年) 1850 1855 1884 1925
针长 投针次数 相交次数 0.80 0.60 0.75 0.83 5000 3204 1030 3408 2532 1218 489 1808
π的估计值 3.15956 3.15665 3.15951 3.14159292
4
数值积分问题
计算积分 α = ∫ f ( x)dx = Ef ( X )
8
U(0,1)随机数的生成
一个简单的随机数生成器:
xi +1 = axi
mod m
ui +1 = xi +1 / m 其中 xi , a, m 均为整数, x0 可以任意选取.
随机数生成器的一般形式: xi +1 = f ( xi ) , ui +1 = g ( xi +1 )
9
源自文库
一个简单的例子
xi +1 = 6 xi mod11, u i+1 = xi +1 /11 (a = 6, m = 11 )
F ( y ) = a + (b a ) y , 0 ≤ y ≤ 1 生成 U (0,1) 随机数 U,则 a + (b - a )U 是来自 U (a, b) 的随机数.
22
几个具体例子(二)
例 2: 设 X ~ exp(θ ) 服从指数分布,则 X 的分布函数为: F ( x) = 1 e
当 x0 = 1 时,得到序列: 1,6,3,7,9,10,5,8,4,2,1,6,3...... 如果令 a = 3, x0 = 1 ,得到序列: 1,3,9,5, 4,1,3,9........ 如果令 a = 3, x0 = 2, 得到序列: 2, 6, 7,10,8, 2, 6.......
每个 xi j 有 m 种选择,所以向量 ( xi 1 , xi 2 ......xi -k ) 可以取 m k -1 个不同的值,所以这样的随机数生成器 的最大周期可以达到 m 生成器的周期.
k -1
1,大大提高了简单同余
15
算法实现
许多程序语言中都自带生成随机数的方法,如 c 中的 random() 函数,Matlab中的rand()函数等. 但这些生成器生成的随机数效果很不一样,比如 c 中的函数生成的随机数性质就比较差,如果用 c ,最好自己再编一个程序.Matlab 中的 rand() 函数,经过了很多优化.可以产生性质很好的随 机数,可以直接利用.
0 1
我们可以将此积分看成 f ( x) 的数学期望.其中 X ~ U [0,1] (均匀分布).于是可以将上式积分看 作是f(X)的数学期望.若{U k ,1 ≤ k ≤ n}为i.i.d .U k ~U[0,1]. 1 n 则可以取α n = ∑ f (U i )作为α的估计, n i =1 由大数定律,可以保证收敛性,即: α n → α with probability as n → ∞ 这表明可以用随机模拟的方法计算积分.
Monte Carlo Simulation Methods (蒙特卡罗模拟方法)
主要内容: 1.各种随机数的生成方法. 2.MCMC方法.
1
从Buffon 投针问题谈起
Buffon 投针问题:平面上画很多平行线,间距为 a.向此平面投掷长 为 l ( l < a) 的针, 求此针与任一平行线相交的概率 p.
一般形式: xi +1 = (axi + c) mod m ui +1 = xi +1 / m 1.可以证明 c 的选择对生成的随机数的均匀性影响
不大,所以为了提高计算速度,一般都令 c = 0 .
2. 线性同余器可以达到的最长周期为 m 1 ,我们 可以通过适当的选择 m 和 a ,使无论选取怎样的 初值 x0 都可以达到最大周期(一般选取 m 为质数)
由于标准正态分布的对称性:Φ (1-u ) = Φ (u )
2
Buffon 投针问题
于是有: l p = P( X ≤ sin Φ ) = ∫ 2 0
π
l sin 2
∫
0
2 2l dxd = πa aπ
π=
2l ap
若我们独立重复地作 n 次投针试验,记 n ( A) 为 A 发生的次数. f n ( A) 为 A
在 n 次中出现的频率.假如我们取 f n ( A) 作为 p = P ( A) 的估计,即 p = f n ( A) .
然后取 π =
2l a .s . 作为 π 的估计.根据大数定律,当 n → ∞ 时, p = f n ( A) → p. af n ( A)
2l P π .这样可以用随机试验的方法求得 π 的估计.历史上 → af n ( A)
3
从而有 π =
有如下的试验结果.
试验者 Wolf Smith Fox Lazzarini
j =1 i
F (ci ) = qi ,我们按照如下步骤生成随机数:
(1) 生成(0,1)上均匀分布随机数 U ;
(2) 寻找 k 满足 qk 1 < U ≤ qk , 1 ≤ k ≤ n; (3) 令 X = cK,显然P(X=ck )=P(q k-1 <U ≤ q k )=p k 则 X 的分布函数是 F ( x).
x /θ
, x≥0
通过计算得 F 1 ( y ) = θ log(1 y ) ,则: X = θ log(1 U ) 服从指数分布(其中 U 服从均匀分布) 又因为 1-U 和 U 有着同样的分布,所以也可以取: X = θ log(U )
23
几个具体例子(三)
例 3 (离散分布随机数的生成) 设 X ~ F ( x) 为取值 c1 , c2 ......cn 的离散随机变量 且 P ( X = ci ) = pi ,令 q0 = 0, qi = ∑ p j ,则有