蒙特卡罗模拟方法

1、1Monte Carlo Simulation Methods(蒙特卡罗模拟方法)主要内容: 1.各种随机数的生成方法. 2.MCMC方法.2从Buffon 投针问题谈起 220, /2 0, sin,:sin.llaXAX随机投针可以理解成针的中心点与最近的平行线的距离X是均匀地分布在区间 上的r.v.，针与平行线的夹角 是均匀地分布在区间 上的r.v.，且X与 相互独立，于是针与平行线相交的充要条件为 即相交 3Buffon 投针问题2sin0022(sin)2lllpP Xdxdaa 于是有：2lap4试验者时间(年)针长投针次数相交次数的估计值Wolf18500.80500025323

2、.15956Smith18550.60320412183.15665Fox18840.7510304893.15951Lazzarini19250.83340818083.141592925数值积分问题1( )()( ) 0,1 . . .1() nniinf x dxEf Xf xXUknii df Unwith probability as n 10kk 计算积分 我们可以将此积分看成 的数学期望。其中 (均匀分布)。于是可以将上式积分看作是f(X)的数学期望.若U ,1为U U0,1.则可以取作为 的估计，由大数定律，可以保证收敛性，即：这表明可以用随机模拟的方法计算积分。6Monte

3、Carlo数值积分的优点与一般的数值积分方法比较，Monte Carlo方法具有以下优点：1. Monte Carlo一般的数值方法很难推广到高维积分的情形，而方法很容易推广到高维情形2/1/22. ()() dO nO n一般的数值积分方法的收敛阶为 ，而由中心极限定理可以保证 Monte Carlo 方法的收敛阶为 。此收敛阶与维数无关，且在高维时明显优于一般的数值方法。7随机模拟计算的基本思路1.针对实际问题建立一个简单且便于实现的概率统计模型，使所求的量（或解）恰好是该模型某个指标的概率分布或者数字特征。2.对模型中的随机变量建立抽样方法，在计算机上进行模拟测试，抽取足够多的随机数，对

4、有关事件进行统计3.对模拟试验结果加以分析，给出所求解的估计及其精度(方差)的估计4.必要时，还应改进模型以降低估计方差和减少试验费用，提高模拟计算的效率8随机数的生成1.蒙特卡罗模拟的关键是生成优良的随机数。2.在计算机实现中,我们是通过确定性的算法生成 随机数,所以这样生成的序列在本质上不是随机 的,只是很好的模仿了随机数的性质(如可以通过 统计检验)。我们通常称之为伪随机数(pseudo-random numbers)。3.在模拟中，我们需要产生各种概率分布的随机数，而大多数概率分布的随机数产生均基于均匀分布U(0,1)的随机数。9U(0,1)随机数的生成一个简单的随机数生成器：1101

5、 mod , , /iiiiixaxmuaxxmxm其中 均为整数, 可以任意选取。111 ( ) , () iiiixf xug x随机数生成器的一般形式：10一个简单的例子1i+116 mod11, u/11 6,11iiixxxam（）0 1 ,x 1,6,3,7,9,10,5,8,4当时 得到序列:,1,6,3.,2.003,1 ,1,3,9.3,2,2,1,3,9,5,42,6,7,10,8,6.axax如果令 得到序列:如果令 得到序列：11一个简单的例子(续)上面的例子中，第一个随机数生成器的周期长度是 10，而后两个生成器的周期长度只有它的一半。我们自然希望生成器的周期越长越好

6、，这样我们得到的分布就更接近于真实的均匀分布。0 (max在给定 的情况下，生成器的周期与和初值种子)选择有关。12线性同余生成器(Linear Congruential Generator )111 () mod /iiiixaxcmuxm一般形式：0 cc 1.可以证明 的选择对生成的随机数的均匀性影响不大，所以为了提高计算速度，一般都令 。02. 1 mmaxm线性同余器可以达到的最长周期为 ，我们可以通过适当的选择 和，使无论选取怎样的初值 都可以达到最大周期(一般选取 为质数)13常用的线性同余生成器Modulus mMultiplier aReference231-1=214748

7、364716807Lewis, Goodman, and Miller39373LEcuyer742938285Fishman and Moore950706376Fishman and Moore1226874159Fishman and Moore214748339940692LEcuyer214748356340014LEcuyer14复杂一些的生成器（一）1.Combining Generators:,1,1,111,11111111111 mod, / 1,2.( 1) mod(1)/ 0 (1)/ 0( )j ijj ijj ij ijJjij ijiiiijJxa xmuxmjJ

8、xxmuxmxmm xmm考虑个简单线性同余生成器：其中为所有中最大的一个15复杂一些的生成器（二）2.Multiple recursive generator1122102(.) mo,.)d/iiiki kikkixa xa xa xmuxxxxm需(要选取种子12-1-1 (,.) 1ijiii kkkxmxxxmm每个有种选择，所以向量 可以取个不同的值，所以这样的随机数生成器的最大周期可以达到 ，大大提高了简单同余生成器的周期。16算法实现许多程序语言中都自带生成随机数的方法，如 c 中的 random() 函数，Matlab中的rand()函数等。但这些生成器生成的随机数效果很不一

9、样，比如 c 中的函数生成的随机数性质就比较差，如果用 c ，最好自己再编一个程序。Matlab 中的 rand() 函数，经过了很多优化。可以产生性质很好的随机数，可以直接利用。17由rand()函数生成的U0,1随机数18由rand函数生成的2维随机点19从U(0,1)到其它概率分布的随机数U(0,1)的均匀分布的随机数，是生成其他概率分布随机数的基础，下面我们主要介绍两种将U(0,1)随机数转换为其他分布的随机数的方法。1.逆变换方法 (Inverse Transform Method)2.舍取方法 (Acceptance-Rejection Method)20Inverse Trans

10、form Method11 ( ) ( )inf :( ), 01( )( )XF xFyx F xyyFyF x设随机变量 的分布函数为，定义：即是的反函数。11 (0,1) () ( ) UXFUF x定理 ：设随机变量服从 上的均匀分布，则的分布函数为 。11( ) ()( )( )( )FUP XxP FUxP UF xF x 由 的定义和均匀分布的分布函明数可证：得：21Inverse Transform Method1 1 ( ) (0,1) ( ) F xUuFu由定理，要产生来自的随机数，只要先产生来自随机数，然后计算 即可。具体步骤如下：(1) (0,1)U上生成均匀分布的随

11、机数 。-1(2) , ( ) ( ) XXUFFx计算则为来自 分布的随机数.22几个具体例子（一）-11 ( , ) 0 ( ) 1, ( )() 01 (0,1) ( - ) ( , ) XU a bxaxaF xaxbbaxbFyaba yyUUab a UU a b例 ： 设，则其分布函数为 ，生成随机数，则是来自的随机数。23几个具体例子（二）/12 exp( ) ( )1, 0 ( )log(1) log(1) ( ) 1 xXXF xexFyyXUUUU 例 ：设服从指数分布，则的分布函数为：通过计算得，则：服从指数分布 其中服从均匀分布又因为 和有着同样的分布，所以也可以取：

12、 log( )XU 24几个具体例子（三）12013 () ( ) ,. () 0, ( ) niiiijjiiXF xc ccP XcpqqpF cq例 离散分布随机数的生成设为取值的离散随机变量且，令 ，则有，我们按照如下步骤生成随机数：(0,1) ;U(1) 生成上均匀分布随机数1 , 1;kkkqUqkn(2) 寻找满足 ( )KXcXF xkk-1kk(3) 令，显然P(X=c )=P(q 0.5,Y=-X,1;exp( 0.5(1) )U U UXUUXUUXUUX 由下述步骤生成Y：生成 上的均匀分布且相互独立的样本 ；令若 ，（生成指数分布的样且则取返回步骤若，且则取返回步骤若

13、，则舍去X,返本）回步骤139随机向量的抽样方法（一）12(,.) nXXX若随机向量 独立，则可用生成一维随机数的方法分别生成向量的每一维即可。121212121121121121(,.) (,.) ( ,.)( ) (|). (|,.)( ),(|). (|,.)nnnnnnnXXXXXXp x xxp x p xxp xx xxF xF xxF xx xx若随机向量 不独立。设概率密度函数为：为对应的分布函数40随机向量的抽样方法（二）1212111211212,. (0,1) ,. () (|,.) 2,3.,.( ,.)nniiinnU UUUXXXF xUF xx xxUinXXX

14、p x xx定理2:设 是独立同分布的变量，是方程组：的解，则 的分布为 1121221. (0,1)(,.),. 2. , . nnnUu uuxXXXxx从分布中独立抽取 解上面的方程组得到样本 生成的基本步骤：41生成多维正态随机数的方法（一）122121(,.) ( , ),11 ( )exp()()2(2 ) |()|nnnTijXXXXnNf xxx 若 服从 维正态分布 则密度函数为:其中 为已知均值向量，为已知协方差阵，为其行列式。12(,.) (0,), ( , )TnnnnZZ ZZNIAAXAZN 若 则可证明：42生成多维正态随机数的方法（二）生成多维正态随机数的具体步

15、骤：1. TAAA 找到一个矩阵 使其满足 122. (0,1) ,.nNnz zz由 分布独立抽取 个随机数123. ( ,.)nxAzzz zz计算 ，其中43v用Monte Carlo方法求解Laplace方程v参见书上5.8节 P213P21544马氏链在马氏链在Monte Carlo随机模拟中的应用随机模拟中的应用定义定义 为要模拟服从给定分布的随机变量,用生成一个易于实现的不可约遍历链 作为随机样本,使其平稳分布为 的方法，称为马氏链蒙特卡罗方法马氏链蒙特卡罗方法.0,nXXn蒙特卡罗方法的一个首要步骤是产生服从给定的概率分布函数 的随机变量(或称为随机样本),由概率论知识,熟知下

16、面的结论.45引理引理 生成随机变量U,使其分布满足U0,1,记为UU0,1, F(x)是给定的一个分布函数,记 为F(x)的反函数,则X=F-1(U)分布函数为F(x).)1 , 0(,)(:sup)(1yyxFxyF46).().1 ();1 (,)()()2(;1 , 0),() 1 (:,),()(, 0)(,)(:xZYZYMgYUUUygYZxxMgxMxgx则返回步骤否则直接返回步骤取若抽取量由下列步骤生成随机变满足和常数的分布函数选取一个易于生成是给定的分布函数假设命题命题47).),()()()(/ )(,)(,)(,)(, )()()(:11nZEhyZEhxdxhnZhy

17、nZEhZZnxZxZEhxdxhPnnkknn大数定律知由的近似值作为积分取则对充分大的存在若样本个独立的只需生成服从的概率分布函数为其中计算积分应用举例48米特罗波利斯(Metropolis)等人在1953年最早给出了通过生成一马氏链实现从分布 中采样(生成相关的样本)这一重要基本思想.随后,哈斯汀(Hastings)将其推广到更一般的形式.下面仅叙述状态空间S为至多可数的情形:49.,0,).1 (,);1 (,),(,1 , 0)2(;),()(),()(, 1min),(,) ,().0,() 1 (:0,).(),),(,),(11平稳分布为是不可约遍历马氏链则返回步骤取否则舍去返

18、回步骤则令若抽取并计算抽取由给定由下列步骤生成为参照矩阵称阵不可约遍历概率转移矩的为任选的易于实现条件而为将要生成的概率分布记nXXXXYYXYXUUUYXTXXYTYYXYXTnSXXnXXSjijiTSiinnnnnnnnnnnnnTT50.)(),()(),()(min),()(),()(, 1min),()(),(),()()().,()()(,),)(,(),(,:jiijjiijijPjijTjjiTijiTiijTjjiTijijiTiPiSjiPjPiSjijijiTPX事实上即性有对称转移概率矩阵其一步是不可约遍历马氏链易证证明梗概.0,),(.)()(,的平稳分布是由上式即

19、证得求和上式两边对nXXSiiPjiSjnSjji51Markov chain Monte Carlo (MCMC)问题提出：( ) ( ) ( ) ( ) ( ) nnRxRf xx dxxx设 为 上的密度函数，在很多问题中我们需要计算高维积分 ，其中常常是非常复杂的多维概率分布，我们无法用前面介绍的方法生成 分布的随机数。因此，我们必须探讨一些新的生成方法。52MCMC方法的基本思路MCMC 是一种简单有效的计算方法，在统计物理，Bayes 统计计算，显著性检验，极大似然估计等领域都有着广泛的应用。基本思路：( ) Markov ( ) xx 通过建立一个平稳分布为的链来得到的样本，基于

20、这些样就可以做各种统计推断。53概率转移核的构造（一）MCMC的方法有很多，在此我们只介绍Metropolis-Hastings方法。基本思路： ( , )() (0()1)( ,) ( ,)( ,) ( ,) ( , )1( ,) ( ,) ( ,) x xqx xp x xq x xx xxxp x xq x xx x dxxxp x x 任意选择一个不可约的转移概率以及一个函数，。对任一组合，定义：则易见构成一个概率转移核。54概率转移核的构造（二） ( ) ( | )( ,)1-( ,)ttxXxqxxxx xxx xxx 此方法的实施比较直观：如果链在时刻处于状态，即，则首先由产生一

21、个潜在的转移，然后根据概率接受作为链下一时刻的状态值，而以概率拒绝转移到，从而链在下一时刻仍然处于状态。 ( )( , )( , )xq 我们的目标是使成为马氏链的平稳分布，下面我们就介绍在给定后，如何选择55概率转移核的构造（三）Metropolis-Hastings 算法：( ) ( , )( ,)min1, ( ) ( ,)x q x xx xx q x x 一个常用的选择：( ,) ( ) ( , )( ) ( ,)( ,)( )( , ) ( ) ( , )( ) ( ,)( )q x xx q x xx q x xp x xxq x xx q x xx q x xx 此时有：56 Markov ( ) ( ,)( ) ( , ) (*)( )Markov x p x xx p x xx 定理：由上述过程产生的链是可逆的，即：且是链的平稳分布。:,*( ) ( , ) ( ) ( ,)( ) ( ,)min 1,( ) ( ,) min ( ) ( ,), ( ) ( , )( ) ( ( ) ( , )minproofxxxxx q x xx p x xx q x