蒙特卡罗算法

1MonteCarloSimulationMethods(蒙特卡罗模拟方法)主要内容:一.M-C方法概述.二.随机数的生成.三.模拟训练.四.实验题目.成信院数学与信息科学系李胜坤2蒙特卡洛(MonteCarlo)方法，或称计算机随机模拟方法，是一种基于“随机数”的计算方法。这一方法源于美国在第二次世界大战中研制原子弹的“曼哈顿计划”。该计划的主持人之一、数学家冯·诺伊曼用驰名世界的赌城—摩纳哥的MonteCarlo—来命名这种方法，为它蒙上了一层神秘色彩。一.M-C方法概述3基本思想很早以前就被人们所发现和利用。17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”。19世纪人们用投针试验的方法来决定π。高速计算机的出现，使得用数学方法在计算机上大量模拟这样的试验成为可能。4从Buffon(蒲丰)投针问题谈起Buffon投针问题:平面上画很多平行线，间距为a．向此平面投掷长为l(la)的针,此针与任一平行线相交的概率p。22[0,/2][0,]sin,{:sin}.llaXAX随机投针可以理解成针的中心点与最近的平行线的距离X是均匀地分布在区间上的r.v.，针与平行线的夹角是均匀地分布在区间上的r.v.，且X与相互独立，于是针与平行线相交的充要条件为即相交52sin0022(sin)2lllpPXdxdaa于是有：2lap若我们独立重复地作n次投针试验，记()nA为A发生的次数。()nfA为A在n次中出现的频率。假如我们取()nfA作为()pPA的估计，即ˆ()npfA。然后取2ˆ()nlafA作为的估计。根据大数定律，当n时，..ˆ().asnpfAp从而有2ˆ()PnlafA。这样可以用随机试验的方法求得的估计。历史上有如下的试验结果。6试验者时间(年)针长投针次数相交次数π的估计值Wolf18500.80500025323.15956Smith18550.60320412183.15665Fox18840.7510304893.15951Lazzarini19250.83340818083.141592927数值积分问题niiixfnxfExxf110)(1)]([d)(选取(0,1)中随机数序列x1，x2，x3，……xn。则误差约，它并不能和一些高级的数值积分算法比拟，但对多维情况，MC方法却很有吸引力。n18niiiizyxfnzyxzyxf1101010),,(1ddd),,(我们可产生一系列随机数可简单取3个随机数构成一个随机点，即,.......,,,,,654321),.......,,(),,,(654321相应地，niibaxfnxxfab1)(1d)(1niiidcbayxfnyxyxfabcd1),(1dd),())((1一般地，A)inpointsrandomnoverof(average)ofmeasure(d)(fAxfA9MonteCarlo数值积分的优点与一般的数值积分方法比较，MonteCarlo方法具有以下优点：1.MonteCarlo一般的数值方法很难推广到高维积分的情形，而方法很容易推广到高维情形2/1/22.()()dOnOn一般的数值积分方法的收敛阶为，而由中心极限定理可以保证MonteCarlo方法的收敛阶为。此收敛阶与维数无关，且在高维时明显优于一般的数值方法。10M-C的基本思路1.针对实际问题建立一个简单且便于实现的概率统计模型，使所求的量（或解）恰好是该模型某个指标的概率分布或者数字特征。2.对模型中的随机变量建立抽样方法，在计算机上进行模拟测试，抽取足够多的随机数，对有关事件进行统计3.对模拟试验结果加以分析，给出所求解的估计及其精度(方差)的估计4.必要时，还应改进模型以降低估计方差和减少试验费用，提高模拟计算的效率11回顾几种连续型分布1.均匀分布U(a,b)其概率密度函数为其他,0,1)(bxaabxf有),(),c(,}{badabdcdxcP其中12均匀性特点:均匀分布随机变量X落在(a,b)内任意子区间的概率只与子区间的长度有关,而与子区间的位置无关.可假设有这种特性的随机变量服从均匀分布.均匀分布随机变量X的取值具有”均匀性”132.正态分布正态分布随机变量X的概率密度函数是Rxxxf],)(21exp[21)(2正态分布由两个参数和唯一确定.其中是X的均值(数学期望):=E(X),它确定了概率曲线的中心位置,而是X的标准差:,它确定了概率曲线的”宽窄”程度.)(XD14在许多实际问题中,有一类随机变量可以表示成为许多相互独立的随机变量之和,而其中每个随机变量对总和只起微小的影响,这类随机变量往往服从或近似服从正态分布.在实际应用中,如果我们分析到一个随机变量受到较多独立的微小因素的叠加影响,就可以用正态分布来模拟这个变量.如:工厂产品的测量尺寸,农作物的收获量,某地区成年人的身高,体重等可看成服从正态分布的随机变量.153.指数分布指数分布随机变量X的概率密度为0,00,)(xxexfx指数分布常用来描述寿命问题.16二.随机数的生成1.蒙特卡罗模拟的关键是生成优良的随机数。2.在计算机实现中,我们是通过确定性的算法生成随机数,所以这样生成的序列在本质上不是随机的,只是很好的模仿了随机数的性质(如可以通过统计检验)。我们通常称之为伪随机数(pseudo-randomnumbers)。3.在模拟中，我们需要产生各种概率分布的随机数，而大多数概率分布的随机数产生均基于均匀分布U(0,1)的随机数。17U(0,1)随机数的生成乘同余法：1101mod,,/iiiiixaxmuaxxmxm其中均为整数,可以任意选取。称为种子,a是乘因子,m是模数0x18一个简单的例子1i+116mod11,u/116,11iiixxxam（）01,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如果令得到序列:如果令得到序列：19一个简单的例子(续)上面的例子中，第一个随机数生成器的周期长度是10，而后两个生成器的周期长度只有它的一半。我们自然希望生成器的周期越长越好，这样我们得到的分布就更接近于真实的均匀分布。0(max在给定的情况下，生成器的周期与和初值种子)选择有关。20线性同余生成器(混合同余法)(LinearCongruentialGenerator)111()mod/iiiixaxcmuxm一般形式：02.1mmaxm线性同余器可以达到的最长周期为，我们可以通过适当的选择和，使无论选取怎样的初值都可以达到最大周期(一般选取为质数)1.c是非负整数.通过适当选取参数c可以改善随机数的统计性质(独立性,均匀性).21常用的线性同余生成器ModulusmMultiplieraReference2^31-1=214748364716807Lewis,Goodman,andMiller39373L’Ecuyer742938285FishmanandMoore950706376FishmanandMoore1226874159FishmanandMoore214748339940692L’Ecuyer214748356340014L’Ecuyer22复杂一些的生成器Multiplerecursivegenerator1122102(......)mo,.......)d/iiikikikkixaxaxaxmuxxxxm需(要选取种子12--1-1(,......)1ijiiikkkxmxxxmm每个有种选择，所以向量可以取个不同的值，所以这样的随机数生成器的最大周期可以达到，大大提高了简单同余生成器的周期。23算法实现许多程序语言中都自带生成随机数的方法，如c中的random()函数，Matlab中的rand()函数等。但这些生成器生成的随机数效果很不一样，比如c中的函数生成的随机数性质就比较差，如果用c，最好自己再编一个程序。Matlab中的rand()函数，经过了很多优化。可以产生性质很好的随机数，可以直接利用。24从U(0,1)到其它概率分布的随机数1.离散型随机数的模拟2.连续型随机数的模拟3.正态随机数的模拟251.离散型随机数的模拟设随机变量的分布律为令将作为区间(0,1)的分点.若随机变量,有X),2,1(}{ipxXPii),2,1(,)(,0)0(1npnPPnii)1,0(~UU),2,1(,)1()()}()1({npnPnPnPUnPPn)}({nP26令则有据此,可得产生的随机数的具体过程为:每产生一个(0,1)区间上均匀分布随机数,若则令取值.nx}{)}()1({nxXnPUnPXnnpxXP}{X)()1(nPUnPU27例1:离散型随机变量X有如下分布律:X012P(x)0.30.30.4设是(0,1)上均匀分布的随机数,令则是具有X分布律的随机数.iiiiUUUx6.0,26.03.0,13.00,0NUUU,,,21Nxxx,,,21282.连续型随机数的模拟a.逆变换方法(常用)(InverseTransformMethod)b.舍取方法(Acceptance-RejectionMethod)定理:设随机变量Y的分布函数F(y)是连续函数,而U是在(0,1)上均匀分布的随机变量,令,则Y与X有相同的分布.)(1UFX2911()(0,1)()FxUuFu由定理，要产生来自的随机数，只要先产生来自随机数，然后计算即可。具体步骤如下：(1)(0,1)U上生成均匀分布的随机数。-1(2),()()XXUFFx计算则为来自分布的随机数.11()()(())(())()FUPXxPFUxPUFxFx由的定义和均匀分布的分布函明数可证：得：30-11~(,)0()1,()()01(0,1)(-)(,)XUabxaxaFxaxbbaxbFyabayyUUabaUUab例：设，则其分布函数为，生成随机数，则是来自的随机数。31/12~exp()()1,0()log(1)log(1)()1xXXFxexFyyXUUUU例：设服从指数分布，则的分布函数为：通过计算得，则：服从指数分布其中服从均匀分布又因为－和有着同样的分布，所以也可以取：log()XU32舍取方法舍取方法(Acceptance-Rejection)方法最早由VonNeumann提出，现在已经广泛应用于各种随机数的生成。基本思路：通过一个容易生成的概率分布g和一个取舍准则生成另一个与g相近的概率分布f。33具体步骤：()()()()1,fxgxfxcgxcx假设和均为集合上的概率密度函数。且满足：。易于从g(x)中生成样本X.用如下步骤生成有Y:1.();2.(0,1),U3.()/()gxXUXUfXcgX生成的样本生成U~且与独立；如果，则取Y=X,返回步骤1，否则舍去X,返回步骤1.34下面我们验证由上述步骤生成的随机数Y确实具有密度函数f(x)()(|()/())(,()/())(()/()()()1/()()()()()()AAAPYAPXAUfXcgXPXAUfXcgXPUfXcgXfxgxdxccgxfxPYAcgxdxfxdxcgxx对于任何的，我们有：）分母＝即Y的概率密度函数为f().35(()/()1/()()(0,1)PUfXcgXcfxcgxU由