蒙特卡罗算法+经典算法

1MonteCarloSimulationMethods(蒙特卡罗模拟方法)主要内容:1.各种随机数的生成方法.2.MCMC方法.2从Buffon投针问题谈起Buffon投针问题:平面上画很多平行线，间距为a．向此平面投掷长为l(la)的针,求此针与任一平行线相交的概率p。22[0,/2][0,]sin,{:sin}.llaXAX随机投针可以理解成针的中心点与最近的平行线的距离X是均匀地分布在区间上的r.v.，针与平行线的夹角是均匀地分布在区间上的r.v.，且X与相互独立，于是针与平行线相交的充要条件为即相交3Buffon投针问题2sin0022(sin)2lllpPXdxdaa于是有：2lap若我们独立重复地作n次投针试验，记()nA为A发生的次数。()nfA为A在n次中出现的频率。假如我们取()nfA作为()pPA的估计，即ˆ()npfA。然后取2ˆ()nlafA作为的估计。根据大数定律，当n时，..ˆ().asnpfAp从而有2ˆ()PnlafA。这样可以用随机试验的方法求得的估计。历史上有如下的试验结果。4试验者时间(年)针长投针次数相交次数π的估计值Wolf18500.80500025323.15956Smith18550.60320412183.15665Fox18840.7510304893.15951Lazzarini19250.83340818083.141592925数值积分问题1()()()~[0,1]}...1()nniinfxdxEfXfxXUkniidfUnwithprobabilityasn10kk计算积分我们可以将此积分看成的数学期望。其中(均匀分布)。于是可以将上式积分看作是f(X)的数学期望.若{U,1为U~U[0,1].则可以取作为的估计，由大数定律，可以保证收敛性，即：这表明可以用随机模拟的方法计算积分。6MonteCarlo数值积分的优点与一般的数值积分方法比较，MonteCarlo方法具有以下优点：1.MonteCarlo一般的数值方法很难推广到高维积分的情形，而方法很容易推广到高维情形2/1/22.()()dOnOn一般的数值积分方法的收敛阶为，而由中心极限定理可以保证MonteCarlo方法的收敛阶为。此收敛阶与维数无关，且在高维时明显优于一般的数值方法。7随机模拟计算的基本思路1.针对实际问题建立一个简单且便于实现的概率统计模型，使所求的量（或解）恰好是该模型某个指标的概率分布或者数字特征。2.对模型中的随机变量建立抽样方法，在计算机上进行模拟测试，抽取足够多的随机数，对有关事件进行统计3.对模拟试验结果加以分析，给出所求解的估计及其精度(方差)的估计4.必要时，还应改进模型以降低估计方差和减少试验费用，提高模拟计算的效率8随机数的生成1.蒙特卡罗模拟的关键是生成优良的随机数。2.在计算机实现中,我们是通过确定性的算法生成随机数,所以这样生成的序列在本质上不是随机的,只是很好的模仿了随机数的性质(如可以通过统计检验)。我们通常称之为伪随机数(pseudo-randomnumbers)。3.在模拟中，我们需要产生各种概率分布的随机数，而大多数概率分布的随机数产生均基于均匀分布U(0,1)的随机数。9U(0,1)随机数的生成一个简单的随机数生成器：1101mod,,/iiiiixaxmuaxxmxm其中均为整数,可以任意选取。111(),()iiiixfxugx随机数生成器的一般形式：10一个简单的例子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如果令得到序列:如果令得到序列：11一个简单的例子(续)上面的例子中，第一个随机数生成器的周期长度是10，而后两个生成器的周期长度只有它的一半。我们自然希望生成器的周期越长越好，这样我们得到的分布就更接近于真实的均匀分布。0(max在给定的情况下，生成器的周期与和初值种子)选择有关。12线性同余生成器(LinearCongruentialGenerator)111()mod/iiiixaxcmuxm一般形式：0cc1.可以证明的选择对生成的随机数的均匀性影响不大，所以为了提高计算速度，一般都令。02.1mmaxm线性同余器可以达到的最长周期为，我们可以通过适当的选择和，使无论选取怎样的初值都可以达到最大周期(一般选取为质数)13常用的线性同余生成器ModulusmMultiplieraReference2^31-1=214748364716807Lewis,Goodman,andMiller39373L’Ecuyer742938285FishmanandMoore950706376FishmanandMoore1226874159FishmanandMoore214748339940692L’Ecuyer214748356340014L’Ecuyer14复杂一些的生成器（一）1.CombiningGenerators:,1,,1,111,11111111111mod,/1,2....(1)mod(1)/0(1)/0()jijjijjijijJjijijiiiijJxaxmuxmjJxxmuxmxmmxmm考虑个简单线性同余生成器：其中为所有中最大的一个15复杂一些的生成器（二）2.Multiplerecursivegenerator1122102(......)mo,.......)d/iiikikikkixaxaxaxmuxxxxm需(要选取种子12--1-1(,......)1ijiiikkkxmxxxmm每个有种选择，所以向量可以取个不同的值，所以这样的随机数生成器的最大周期可以达到，大大提高了简单同余生成器的周期。16算法实现许多程序语言中都自带生成随机数的方法，如c中的random()函数，Matlab中的rand()函数等。但这些生成器生成的随机数效果很不一样，比如c中的函数生成的随机数性质就比较差，如果用c，最好自己再编一个程序。Matlab中的rand()函数，经过了很多优化。可以产生性质很好的随机数，可以直接利用。17由rand()函数生成的U[0,1]随机数18由rand函数生成的2维随机点19从U(0,1)到其它概率分布的随机数U(0,1)的均匀分布的随机数，是生成其他概率分布随机数的基础，下面我们主要介绍两种将U(0,1)随机数转换为其他分布的随机数的方法。1.逆变换方法(InverseTransformMethod)2.舍取方法(Acceptance-RejectionMethod)20InverseTransformMethod11()()inf{:()},01()()XFxFyxFxyyFyFx设随机变量的分布函数为，定义：即是的反函数。11(0,1)()()UXFUFx定理：设随机变量服从上的均匀分布，则的分布函数为。11()()(())(())()FUPXxPFUxPUFxFx由的定义和均匀分布的分布函明数可证：得：21InverseTransformMethod11()(0,1)()FxUuFu由定理，要产生来自的随机数，只要先产生来自随机数，然后计算即可。具体步骤如下：(1)(0,1)U上生成均匀分布的随机数。-1(2),()()XXUFFx计算则为来自分布的随机数.22几个具体例子（一）-11~(,)0()1,()()01(0,1)(-)(,)XUabxaxaFxaxbbaxbFyabayyUUabaUUab例：设，则其分布函数为，生成随机数，则是来自的随机数。23几个具体例子（二）/12~exp()()1,0()log(1)log(1)()1xXXFxexFyyXUUUU例：设服从指数分布，则的分布函数为：通过计算得，则：服从指数分布其中服从均匀分布又因为－和有着同样的分布，所以也可以取：log()XU24几个具体例子（三）12013()~(),......()0,()niiiijjiiXFxcccPXcpqqpFcq例离散分布随机数的生成设为取值的离散随机变量且，令，则有，我们按照如下步骤生成随机数：(0,1);U(1)生成上均匀分布随机数1,1;kkkqUqkn(2)寻找满足()KXcXFxkk-1kk(3)令，显然P(X=c)=P(qUq)=p则的分布函数是。25标准正态分布随机数的生成正态分布是概率统计中最重要的分布，在此我们着重讨论如何生成标准正态分布随机数。引理：121/21121/221212,(0,1)(2ln)cos(2)(2ln)sin(2)UUUZUUZUUZZ设是独立同分布的变量，令：则证明：与独由随机立，且均服从向量的变换公标准式易正态分布。得，略。26Box-Muller算法121.,UU生成两个独立的均匀分布随机数1/21121/221212(2ln)cos(2)(2ln)sin(2)ZUUZUUZZ2.令则，为独立的正态随机数27逆变换方法（一）我们无法通过具体的数学表达式计算正态分布函数的逆函数，我们必须通过数值的方法逼近正态函数下面我们介绍Beasley-Springer-Moro方法。111))0.51uuu－－－由于标准正态分布的对称性：（1-（所以我们只需计算时的值即可。28逆变换方法（二）32103200.50.92(0.5)()1(0.5)nnnnnnuauubu-1当采用有理函数逼近：8100.921Chebyshev()[log(log(1))]nnnuucu当采用逼近：29逆变换方法（三）,,nnnabc在上面的逼近中均为常数。值由下表给出：a0=2.50662823884b0=-8.47351093090a1=-18.61500062529b1=23.08336743743a2=41.39119773534b2=-21.06224101826a3=-25.44106049637b3=3.13082909833c0=0.3374754822726147c5=0.0003951896511919c1=0.9761690190917186c6=0.0000321767881768c2=0.1607979714918209c7=0.0000002888167364c3=0.0276438810333863c8=0.0000003960315187c4=0.0038405729373609在matlab中可以直接通过norminv()函数直接计算标准正态分布函数的逆。30Matlab生成的正态随机数31Acceptance-RejectionMethod（一）Acceptance-Rejection方法最早由VonNeumann提出，现在已经广泛应用于各种随机数的生成。基本思路：通过一个容易生成的概率分布g和一个取舍准则生成另一个与g相近的概率分布f。32Acceptance-RejectionMethod（二）具体步骤：()()()()1,fxgxfxcgxcx假设和均为集合上的概率密度函数。且满足：。易于从g(x)中生成样本X.用如下步骤生成有Y:1.();2.(0,1),U3.()/()gxXUXUfXcgX生成的样本生成U~且与独立；如果，则取Y=X,返回步骤1，否则舍去X,返回步骤1.33Acceptance-RejectionMeth