蒙特卡罗方法(II)

蒙特卡罗方法（二）§3随机变量的抽样---离散型直接抽样法111nniiiipp随机变量X：｛x1,x2,x3,···,xN｝例如：x可取3个值x1,x2,x3，出现的几率分别为2/8,5/8,1/8，则随机数小于2/8时实现x1，在区间[2/8,7/8]中时实现x2，大于7/8时实现x3.概率密度f：｛p1,p2,p3,···,pN｝如果从[0,1]区间中均匀抽样得到的随机数ξ满足下式时则物理量x取值为xn。实际需要的大多数随机变量并不是[0,1]区间均匀分布的，而是有各种不同形式分布密度函数的随机变量。因此，对不均匀随机变量抽样的关键问题是如何从均匀分布的伪随机变量样本中，抽取符合所要求的分布密度函数的简单子样。离散型直接抽样法§3随机变量的抽样---离散型直接抽样法!nnpen例如Poisson分布是离散型分布第n个事件发生的条件为100!!iinniieii3MeV光子入射屏蔽铅板的全吸收反射过程反应类型X：光电效应康普顿散射电子对产生反应截面σ：σ1σ2σ3反应几率f：ε1=σ1/σε2=σ2/σε3=σ3/σ归一化：123123100%电子对产生光电效应康普顿效应NNYY1i2ii设连续型变量x在区间[a,b]中取值，将上述的离散情形取连续极限:limlim0nNNbaxN()(')'xaxpxdxξ(a)=0,ξ(b)=1且是单调增加要求变量x，可由上式解析反解出x(ξ)的函数表达式，即求反函数。这对一些简单的几率密度函数解析表达式是很容易做到的。§3随机变量的抽样---连续型直接抽样法求累积函数p(x)已归一化：()1,()1bapxdxpx未归一化的物理量（如截面σ(x)）：()()()baxpxxdx§3随机变量的抽样---连续型直接抽样法如粒子随机运动的自由程分布为指数分布：10()0xexpxotherwise101xtxedte求反函数ln(1)lnx注意（1-ξ）和ξ同样服从[0，1]的均匀分布求累积函数§3随机变量的抽样---变换抽样法()()()()dypxdxgydypxgydx变换抽样法的基本思想是将一个比较复杂的分布p(x)的抽样，变换为已知的简单分布g(y)的抽样找到x↔y之间的对应关系，使得几率密度守恒：101()0ygyotherwise显然，当g(y)取[0,1]均匀分布时，问题即化为：寻找y(x)，使其导数为p(x)，然后在[0,1]区间中对变量y抽样得到均匀分布的随机数，再由x(y)关系得到对应几率密度函数p(x)的随机抽样x。一维：二维：有两个变量x和y的联合分布密度函数为p(x,y)，欲变换至变量u和v，它们的联合分布密度函数为g(u,v)(,)(,)(,)(,)(,)uvpxydxdyguvdudvguvdxdyxy(,)(,)(,)(,)uvpxyguvxy取联合分布密度函数g(u,v)为均匀分布：101,01(,)0ifuvguvotherwise§3随机变量的抽样---变换抽样法则只要寻找变换式x=x(u,v),y=y(u,v)，以使p(x,y)=|∂(u,v)/∂(x,y)|对均匀随机变量(u,v)进行抽样，代入变换式得x和y的抽样。对于Gauss正态几率分布的抽样代换令极坐标系下的角度为2πv，半径为22()/21()2xxpxexxx2/21()2xpxe2lncos2xuv2lnsin2yuv试通过一个两维联合分布的抽样获得该一维分布的抽样。u和v都是[0,1]区间中的均匀分布的随机抽样，则变换关系式为2lnu§3随机变量的抽样---变换抽样法22()/2(,)1()()(,)(,)2xyuvepxpypxyxyJacobian即两维分布为两个独立分布之积。显然抽样x或y都满足正态分布。可见，为了得到满足一个复杂分布的随机抽样，可以用了两个满足简单分布的随机数u和v抽样------Box-Muller法。可得反变换22exp[]2xyu1tan()2yvarcx§3随机变量的抽样---变换抽样法2lncos2xuv2lnsin2yuv§4蒙特卡罗方法应用---数值积分0()banSfxdxSN计算定积分的基本公式S0为方形区域的面积，N是总点数，n是掷入f(x)下面积区域中的点数。例如，用MonteCarlo方法求π值。产生一对在[0,1]区间中均匀分布的随机数作为点的坐标(x,y)值，判断条件x2+y2≤1是否成立，成立则计数n值，当总点数N足够大时：lim4NnN§4蒙特卡罗方法应用---数值积分（1）掷石法•具体FORTRAN程序如下：•programmain•implicitnone•real(kind=8),parameter::m=1.0e8•integer::n,counter=0•real(kind=8)::x(1:m),y(1:m),Pi,r•open(unit=10,file='random_num.txt')•callrandom_seed()•CALLrandom_number(x)•callrandom_number(y)•don=1,m•r=sqrt(x(n)**2+y(n)**2)•if(r=1)then•counter=counter+1•endif•enddo•pi=(counter*1.0)*4.0/(m*1.0)•write(10,*)'ThePiis',pi•endprogrammain附录：程序§4蒙特卡罗方法应用---数值积分0510152025303.023.043.063.083.103.123.143.163.183.20N*10000设想一个物理系统由相互作用的多个粒子组成，如凝聚态物质中的原子、原子中的电子，而每个粒子都可以有几个自由度。如对具有相互作用的m个原子组成的气体，经典配分函数：假设每个坐标取10个点，总共将有103m个网格点，对于m=20，如用万亿次计算机进行计算，需要1048秒（约为宇宙年龄的1029倍）。1exp[()]mijijZdrdrUr1BkT§4蒙特卡罗方法应用---数值积分这个式子只有当m相当小时才可能用数值积分方法计算出来。MC方法真正的威力在于计算多重积分。因此，对于高维的多重积分不应采用固定网格法。若采用简单抽样的MC方法，则高维积分式可写为：121212121,211(,,)1[()](,)nnbbbnnaaanNjjiinijidxdxdxfxxxbafxxxN其中对每个坐标的抽样值是在相应的区间范围内均匀抽取的。§4蒙特卡罗方法应用---数值积分如果积分中的被积函数不光滑，如配分函数式中的势能变化急剧，则用目前的非权重简单抽样MC方法得到的精度可能很差，有一些极大值区域没有被随机抽出。统计力学中Boltzman分布函数在相空间的大部分区域其值都是很小的，真正有贡献的区域范围很窄。简单抽样中如果要提高精度，必须增大抽样量。这种办法是不可取的，因为这样的计算效率很低。引入带权重的抽样法，将既可以保证计算精度又有很高的效率-----重要抽样法。§4蒙特卡罗方法应用---数值积分此外将被积函数变为平坦也可以有效地提高MC方法的精度-----提取法。对于一维积分，设能够构造一个与被积函数f(x)形状相似的函数g(x)，且它的积分值已知运用上面所述的MC方法于一个较为平坦的被积函数f(x)−g(x)()(),()bafxgxgxdxJ()[()()]bbaafxdxfxgxdxJ提取法§4蒙特卡罗方法应用---数值积分在数值计算方法中，奇点必须避开，但因为奇点附近的贡献是相当大的，不可忽略，因此要尽可能的逼近奇点，这在固定网格划分形式下的其它数值方法中是很难做到的，奇异积分指在积分限内被积函数有发散点，尽管在奇点上被积函数是无穷大，但积分值可以是有限的。120121dxx如被积函数在上限处是发散的但用MC方法就比较容易，这是因为，随机选取的x恰为奇点的可能性不大，还可以在包含奇点的附近选择一非常小的排斥区域以避开奇点。不过，简单抽样的效率很低。§4蒙特卡罗方法应用---数值积分重要抽样法设有一个几率分布g(x)与f(x)形状相似：则积分可写成()~1,()1()bafxgxdxgx()()()()bbaafxfxdxgxdxgx§4蒙特卡罗方法应用---数值积分101()0ygyotherwise其中选择y是在[0,1]区间中的均匀分布的，作变量代换x↔y故y(x)即为累积函数，则得重要抽样法的精神是：将随机变量x变换为另一随机变量y，而变换关系y(x)是表征被积函数f(x)曲线变化特征的g(x)的累积函数，对y的简单抽样就是对x带权重g(x)的抽样。g(x)不是归一化的其中对y的抽样是在[0,c]区间。(),(')'xadygxdxygxdx10()()()()()()bbaafxfyfxdxgxdxdygxgy()bacgxdx0()()()bcafyfxdxdygy§4蒙特卡罗方法应用---数值积分重要抽样法下的MonteCarlo积分为：1()1()()NbiaiifxffxdxgNgx例如，计算积分式中的指数项随x增加而迅速减小，g(x)必须是另一有相同行为的函数，用指数函数在x=0附近的渐近展开来近似：它在积分区间内是正的：考虑到y(0)=0，故上式中括号内应取负号原积分写成：340exp[(11)]1yIdyy2(11)xy20'()(1)2(11))24xxxyxdxxxy23()1128482xxxxgx10exp()2xIdx§4蒙特卡罗方法应用---数值积分经典统计力学考虑的是一个多自由度的力学体系，这些自由度一般是粒子的坐标和动量，或者是磁矩即自旋。经典体系是指这些自由度是可对易的，以这些自由度为坐标展开的空间即为相空间。如：N个粒子体系存在6N维相空间3N个位置坐标(q1,…,q3N)3N个动量(p1…p3N)相空间中的一个点代表力学体系的微观状态，相应的6N个坐标组成体系的一个构型。§4蒙特卡罗方法应用---经典粒子体系每个坐标qi和动量pi随时间的演变由正则方程决定,iiiiHHqppq体系的Hamilton量H(q,p)，是所有粒子坐标和动量的函数。（1）随着时间的改变，坐标值(q,p)连续地变化，相空间中的代表点走出一条轨迹。当力学体系从不同的初态出发时，一般情况下在相空间中就沿着不同的轨迹而运动，这些轨迹是不相交的。但当H显含时间t时，经过某一点的轨迹在不同时刻的方向可以不同，因此轨迹可能相交。（2）代表点的轨迹只能在相空间中的一个有限区域内运动，因为有限的体积V限制了q的取值范围，而有限的能量E限制了q和p。一个力学体系的代表点在相空间中随时间变化所走过的轨迹（3）对于能量守恒的保守系统，轨迹限于在相空间中由H(q,p)=E确定的曲面上运动。如果总能处于(E,E+ΔE)的一个区域范围内，则轨迹限制于一个厚为ΔE的曲面壳层里。§4蒙特卡罗方法应用---经典粒子体系统计系综测量仪器的响应很慢，测量结果是多个粒子在一段长时间内作用的平均效应；即使可以作瞬时测量，其瞬时值与平均值差别也是很小的（对于N∼1023的体系，涨落误差是N-1/2）（4）统计力学的中心思想是用几个宏观物理量代替6N个描述微观状态的自由度。而宏观量可以对所有粒子的微观参量求平均值而得到。§4蒙特卡罗方法应用---经典粒子体系对于由微观粒子组成的宏观系统，一般不必知道所有粒子微观力学量（1）没有必要了解那么详细，对实际问题没