蒙特卡罗方法与拟蒙特卡罗方法的历史、现状及展望

朱 辉, 刘义保, 游 运

(东华理工大学,江西抚州 344000)

蒙特卡罗方法与拟蒙特卡罗方法的历史、现状及展望

朱 辉, 刘义保, 游 运

(东华理工大学,江西抚州 344000)

随着非线性科学的发展,许多物理、工程技术和数学模型都可以转化为非线性方程,如非线性常微分方程、偏微分方程等。非线性方程的求解已经成为非线性科学领域的一个重要研究课题。Zakharov-Kuznetsov方程(简称 ZK方程)作为非线性方程中重要的一类,是由 Zakharov和 Kuznetsov在 1974年提出的,该方程是 KdV方程在二维空间的典型推广形式之一,因此研究该方程具有广泛的理论意义和实践意义。本文用拓展的双曲函数正切法,借助 Riccati方程的解,结合Mathematical数学软件,得到 Zakharov-Kuznetsov方程新的显示精确解,包括周期解和孤立波解.所给的方法还可以用来求解其它的一大类非线性发展方程。

Zakharov-Kuznetsov方程;Riccati方程;周期解;孤立波解

蒙特卡罗 (Monte Carlo,MC)方法,又称随机抽样或统计试验方法,其前身为Metropolis算法,是由Neumann,Ulam和Metropolis在第二次世界大战末为了解决在核武器实验与研制中遇到的裂变物质中子扩散问题而提出的。“蒙特卡罗”由Metropolis(1949)命名,并正式使用。提出蒙特卡罗方法的初衷是用于物理数值模拟问题,后来随着计算机的快速发展,这一方法很快在函数值极小化、计算几何、组合计数等方面应用 (雷桂媛,2003),于是它作为一种独立的方法被提出来,并发展成为一门新兴的计算科学,属于计算数学的一个分支。如今MC方法已是求解科学、工程和科学技术领域大量应用问题的常用数值方法。

MC方法对随机过程和高维积分等问题求解的有效性是公认的,随着拟随机序列的出现,蒙特卡罗方法也已经发展到拟蒙特卡罗方法。事实上,蒙特卡罗和拟蒙特卡罗 (Quasi-Monte Carlo,QMC)方法可能比较耗时,但它们是目前解决高维问题唯一切实可行的方法。随着现代计算机技术的飞速发展,已可以计算越来越复杂的问题。

1 随机数介绍

随机数的产生是 MC模拟的核心问题之一。介绍MC和 QMC方法前,有必要先解释真随机数、伪随机数和拟随机数的概念。

在连续型随机变量的分布中,最简单而且最基本的分布是[0,1]区间上的均匀分布,也称单位均匀分布。由该分布抽取的简单子样ξ1,ξ2,……称为随机数序列,其中每一个体称为随机数,有时称为标准随机数或真随机数(standard random or true random numbers),独立性和均匀性是其必备的两个特点。真随机数是数学上的抽象,真随机数序列是不可预计的,因而也不可能重复产生两个相同的真随机数序列。真随机数只能用某些随机物理过程来产生,如放射性衰变、电子设备的热噪音、宇宙射线的触发时间等。采用随机物理过程来产生MC计算用的随机数,理论上不存在什么问题,但在实际应用时,要做出速度很快而又准确的随机数物理过程产生器是非常困难的。Frigerio等在 1975年做过相关工作 (马文淦,2005),都柏林圣三一学院计算机科学和统计学院从 1998年开始提供随机数的在线生成服务,它是由大气噪声产生。而如何在计算机上由其硬件产生真随机数也是目前很多学者探讨的问题。

实际使用的随机数通常都是用某些数学公式产生的伪随机数 (pseudo-random numbers)。要把伪随机数当真随机数来使用,必须要通过随机数的一系列的统计检验。一个经典且现在仍然广泛使用的方法是线性同余方法(linear congruentialmethod),如 Lehmer(1951)提出的乘同余方法和 Roten-berg(1960)提出的乘加同余方法等。好的伪随机数发生器具备以下特征:良好的统计分布特性,高效率的伪随机数产生,循环周期长,产生程序可移植性好和伪随机数可以重复产生,其中满足均匀性、独立性等统计性质是最重要的。

无论伪随机数用什么方法产生,它的局限性在于这些随机数总是一个有限长的循环集合,而且序列偏差的上确界达到最大值。所以若能产生低偏差的确定性序列 (数学上称为 Low Discrepancy Sequences)是很有用的,产生的序列应该具有这样的性质,即任意长的子序列都能均匀地填充函数空间。人们已经产生了若干种满足这个要求的序列,如 Halton序列、Faure序列、Sobol’序列和 Niederreiter序列等。称这些序列为拟随机数 (Quasi-Random Numbers)序列。伪随机序列是为了模拟随机性,而拟随机序列更致力于均匀性。

从二维单位正方形中抽取的三种随机数,图 1为真随机数 (来源于 RANDOM.ORG),图 2为伪随机数 (Matlab中的 rand函数生成),图 3为拟随机数(Halton序列),显示伪随机数更接近于真随机数,而拟随机数则比它们分布的更加均匀。

图 3 拟随机数Fig.3 Q uasi-random number

2 蒙特卡罗方法

2.1 基本原理

原理。若要计算某个量 I,可以任意选择一个数学期望为 I的随机变量γ,从中抽出子样γ1,

2.4 任意分布的抽样

由MC原理可知,当确定γ后,关键的就是从γ的分布中抽取子样γ1,γ2,…,γN。因此,随机变量抽样是MC方法的最重要问题之一。对于任意非单位均匀分布随机变量γ的抽样,均是使用严格数学方法,借助随机数产生,步骤为先抽取若干个随机数ξ1,ξ2,…,ξm,然后经过某种变换 (或关系或运算 )g(ξ1,ξ2,…,ξm)得到 γ子样的一个个体。主要抽样方法有直接抽样、舍选抽样、复合抽样、加抽样、减抽样、乘抽样、乘加抽样、乘减抽样、对称抽样、积分抽样、偏倚抽样、近似-修正抽样、分层抽样、系统抽样、相关抽样、分裂和轮盘赌抽样、等概率抽样、驿站重要抽样、拉丁超立方抽样等 (Mckay et al.,1979;裴鹿成等,1980;徐钟济,1985;马文淦,2005;许淑艳,2006)。如何减小方差提高效率是抽样中的一个重要问题。

2.5 MC方法应用与发展

蒙特卡罗方法首先在核领域取得成功,但早期由于计算费用昂贵及性价比不高,MC方法只是作为确定论方法的补充。1970年代,组合几何的发展,解决了三维复杂几何和反应堆堆芯结构的描述,使MC方法理论得到重要发展①邓力.2009.输运问题蒙特卡罗模拟现状概述[J].第十届全国蒙特卡罗方法学术交流会资料[C].。21世纪采用连续点截面的Monte Carlo中子输运与燃耗耦合计算程序MCNP,fluka,EGS,Geant4等对全堆的模拟,进一步巩固了MC方法在核领域的地位,并已成为模拟各种粒子输运问题、反应堆临界安全分析、燃料循环、辐射输运计算的首选工具。输运问题中的反应堆分析、辐射输运、核探测领域是MC方法发展的新领域(刘立坡等,2007;刘立坡等,2008;Zhang Tao et al.,2009)①邓力.2009.输运问题蒙特卡罗模拟现状概述[J].第十届全国蒙特卡罗方法学术交流会资料[C].②张涛,刘义保,杨波.2009.平行 gamma射线束深穿透的蒙特卡罗模拟[J].第十届全国蒙特卡罗方法学术交流会资料[C].。随着计算机的快速发展,MC方法在统计物理、生物医学、核医学、量子力学、分子动力学、仿真学、仪器设计与校正、可靠性设计、石油测井、物探、金融、信息和航天、参数估计等领域内也有广泛应用。

3 拟蒙特卡罗方法

与Monte Carlo方法相似,但理论基础不同的方法 ——“拟蒙特卡罗方法或准蒙特卡罗方法”近年来也获得迅速发展。我国数学家华罗庚、王元提出的“华-王”方法即是其中的一例。这种方法的基本思想是用确定性的超均匀分布序列 (数学上称为Low Discrepancy Sequences)代替Monte Carlo方法中的随机数序列。对某些问题该方法的实际速度一般比MC方法快,并能计算精确度。

3.1 基本原理

点列的均匀性可用偏差来度量,用的比较多的是所谓的星偏差 (star discrepancy):设{X1,X2,…,XN}为 s维的单位立方体 [0,1)s中的点列,μ(V)为[0,1)s中区域 V的勒贝格测度,M(X1,X2,…,XN)为点列{X1,X2,…,XN}含在 V中的点数,则此点列的星偏差定义为

其中 V(f)为 Hardy and Krause意义下 f的有界变差,为 点 列 {X1,X2,…,XN}的 星 偏 差。Koksma-Hlawka不等式给出了QMC积分的误差限,同时说明所用点列偏差越小,结果越精确。随机数序列的偏差期望值可限制在 (log(logN))N-1/2,如果用低偏差序列代替随机数序列可以提高精度。而最小偏差可达 o((logN)s-1N-1)(Wang et al.,1994),当维数 s较大时QMC积分比MC积分的收敛速度 o(N-1/2)更快。QMC方法的优点是排除了 MC方法的随机性,给出一致分布得“最均匀”的确定点列,因而得出的精确度是确定的,而不再是概率的。

3.3 低偏差序列

应用QMC方法的前提是要产生“超均匀分布”的低偏差序列,统计学家们提出了许多种不同的低偏差序列,有 Halton(1960)序列,Sobol’(1967)序列,Faure(1982)序列和 Niederreiter(1987)的 (t,s)序列等。下面介绍 Halton序列的构造:

Halton序列是 Van der Corput序列的推广。Halton序列是通过将一系列整数表示成某个基(base)的数位 (digit)形式,然后将这些数位按反序排列再在前面加小数点而得到的值。将 s维Halton序列表示成 X1,X2,…,其中每一个随机数是一个 s维向量,即 Xi=(xi1,xi2,…,xis)。生成 Halton序列的步骤如下:首先选择 s个基 b1,b2,…bs,比如选择前 s个素数。然后对某个整数m,将m表示成以bj为基的数位,再将这些数位按反序排列再在前面加小数点得到新的值,便是序列中某个随机数向量的第 j个元素。即对某个整数 m,有以下步骤:

(3)置 m=m+1,然后重复 1,2两步。

举个例子,比如 s=2,即维数是 2,设m=7,选择 2和 3作为基。有 7=1112,7=213,于是便有第一个随机数向量 = (0.1112,0.123),即(0.875 000,0.555 555 6)将 m加 1,如此循环,便得到整个序列。图 3是以 2和 3为基,整数 m=1开始所得二维 Halton序列。

3.4 QMC方法的应用、缺陷与改进

拟蒙特卡罗方法已成功应用的领域有 QMC积分、计算机仿真实验、全局最优化等 (AndreaReese,2009)。但是QMC积分方法的误差公式只给出了上限且难以计算,因此常常不能用来估计误差。从而QMC积分只对有限的几类函数比MC积分有效,这也是基于积分更侧重样本点的均匀性而非随机性。对于用拟随机数来进行模拟研究,QMC方法却不是很奏效,现在用来模拟的领域主要是在金融中的高维问题 (Tan et al.,2000),这是一方面由于模拟时无偏估计是至关重要的,同时要对误差进行估计,而另一方面是因为随机性的要求导致的。

有一种改进方法是 Randomized Quasi-Monte Carlo(RQMC)methods,是将低偏差序列与MC方法结合,其思想是使超均匀分布序列随机化,使其服从均匀分布,又保留了它的高度均匀性。Cranley等(1976)第一次提出将拟随机数随机化,很多文献探讨了随机化 (Joe,1990;Faure,1992;Owen,1997b;Owen,1998a; KenSengTan et al.,2000;Braaten et al.,1979;Owen et al.,1997;Owen,1997a)及部分随机化 (KenSengTan et al.,2000)并应用于实际模拟,取得较好效果。

4 总结与展望

蒙特卡罗方法是对问题进行随机采样,所以它不是“确定论”的解法,是与确定性算法相对应的方法。过去主要是应用于核领域,但随着计算机的快速发展以及并行算法的实现,MC方法的应用领域也在迅速拓展 (Jung et al.,2001;Martin,2001;Frenkel et al.,2002;DavidBoas et al.,2002;WangXiaoqun et al.,2003;Bruno,2005;Aline et al.,2008;Belen et al.,2009;Chen et al.,2009;Chin et al.,2009;Zhang et al.,2009)。而作为加速MC收敛速度的 QMC方法,在更关注样本空间均匀性的应用领域,更加显现其优越性,已为越来越多的计算科学家和应用工作者研究和使用。

华罗庚,王元著.1963.数值积分及其应用[M].科学出版社.

雷桂媛.2003.关于蒙特卡罗及拟蒙特卡罗方法的若干问题研究[D].浙江大学博士学位论文.

刘立坡,刘义保.2007.特征 X射线激发效率的蒙特卡罗模拟[J].核技术.30(08):672-674.

刘立坡,刘义保,王娟.2008.积累因子影响因素的蒙特卡罗模拟[J].辐射防护.28(02):51-54.

马文淦.2005.计算物理学[M].北京:科学出版社.

茆诗松,王静龙,濮晓龙.1998.高等数理统计[M].北京:高等教育出版社,施普林格出版社.

裴鹿成,张孝泽.1980.蒙特卡罗方法及其在粒子输运问题中的应用[M].北京:科学出版社.

魏宗舒.2005.概率论与数理统计教程[M].北京:高等教育出版社.徐钟济.1985.蒙特卡罗方法[M],上海:上海科学技术出版社.

许淑艳.2006.蒙特卡罗方法在实验核物理中的应用[M].北京:原子能出版社.

Aline B,Koen L.2008.Monte Carlo localization for mobile wireless sensor networks[J].Ad Hoc Networks,6(5),718-733.

Andrea R.2009.Random number generators in genetic algorithms for unconstrained and constrained optimization[J],Nonlinear Analysis 71:e679-e692.

Belen J,Rafae lMiro,JuanM.Campayob,SergioDez,GumersindoVerdu.2009.Radiotherapy treatment planning based on Monte Carlo techniques[J].Nuclear Instruments andMethods in PhysicsResearch A.1-6.

Braaten E,Weller G.1979.An improved low discrepancy sequence for multidimensional quasi-Monte Carlo integration[J].Journal of Computational Physics,33:249-258.

Bruno B.2005.Associating domain-dependent knowledge and Monte Carlo approacheswithin a Go program[J].Infor mation Sciences,175(4):247-257.

Caflisch R E.1998.Monte carlo and quasi-monte carlo methods[J].Acta,Numerica,7:1-49.

ChinM PW,Spyrou N M.2009.A detailedMonte Carlo accounting of radiation transport in the brain during BNCT[J].Applied Radiation and Isotopes 67:164-167.

Cranley R,Patterson,TN L.1976..Randomization of number theoretic methods formultiple Integration[J].S.I.A.M Journal of NumericalAnalysis 23,904-914.

David B,Culver J,Stott J,Dunn A.2002.Three dimensional Monte Carlo code for photon migration through complex heterogeneousmedia including the adult human head[J].Optics Express,10(3):159-170.

Faure H.1992.Good permutation for extreme discrepancy[J].Journal ofNumber Theory 42,47-56.

FrenkelD,MulderB M.2002.The hard ellipsoid-of-revolution fluid I.Monte Carlo simulations[J].Molecular Physics,100(1):201-217.

Halton J H.1960.On the efficiencyof certain quasi-random sequencesof points in evaluation rnulti-dimensional integrals[J].Numer.Math.,2:84-90.

Halton J H,Smith,Algorith G B.1964.Radical-inverse quasi-random point sequence[J].Comm.ACM,7:701-702.

Henderson S G,NelsonB L.(Eds.),Chapter12 Quasi-Random Number Techniques Handbook in OR&MS,13.http://www.random.org/

Joe S.1990.Randomization of lattice rules for numerical multiple integration[J].Journal of Computational and Applied Mathematics 31,299-304.

Jung H,Salam G P.2001.Hadronic final state predictions from CCFM:the hadron-levelMonte Carlo generator Cascade[J].Theoretical physics,19:351-360.

LehmerD H.1951.Mathematicalmethods in largescal computing units.Proc.Sym.on largescal digital calcul[J].machinery,Harvard Univ.Press:141.

Martin H.2001.Speeding up the hybridMonte Carlo algorithm for dynamical fermions[J].PhysicsLettersB,519(1-2):177-182.

MckayM D,Beckman R J,ConoverW J.A.1979.comparision of three methods for selecting valuesof input variables in the analysisof output from a computer code[J].Technometrics,21,2.

MetropolisN,Ulam S.1949.The monte carlo method[J].Arn.stat.Ass.,44:335-341.

Niederreiter H.1978.Quasi-monte carlo methods and pscudo-random numbers[J].Bull.Amer.Math.Soc.,84(6):957-1041.

Niederreiter H.1992.Random Number Generation and Quasi-Monte CarloMethods[J].S IAM CBMSNSF Regional Conference Series in AppliedMathematics,63.S IAM,Philadelphia.

Owen A B.1997a.Monte Carlo variance of scrambled net quadrature[J].S.I.A.M Journal ofNumericalAnalysis 34(5),1884-1910.

Owen A B.1997b.Scrambled net variance for integrals of smooth functions[J].Annals of Statistics 25(4),1541-1562.

Owen A B,TavellaD.1997.Scrambled nets for value at risk calculations[J]. In:Grayling,S. (Ed.),VaR:Understanding and Applying Value-at-Risk.R ISK,London.

Owen A B.1998a.Latin supercube sampling for very high-dimensional simulations[J].ACM Transactions onModeling and Computer Simulation 8(1):71-102.

RotenbergA.1960.A new pseudo-random number generator[J].Assoc Comp Mach,7.

Chen S W,Xiaowei Liu,Chunxiang Zhang,Qiang Tang.2009.The Monte Carlo simulation of the absorbed dose in quartz[J].Radiation Measurements,44:626-628.

Sobol IM.1998.On quasi-Monte Carlo integrations.Mathematics and Computers in Simulation[J]:47.

SobolL M.1967.The distribution of points in a cube and the approximate evaluation of integrals[J].USSR Comp.Math.and Math.Phys.,7:86-112.

Tan K S,Phelim P,Boyle.2000.Applications of randomized low discrepancy sequences to the valuation of complex securities[J],Journal of Economic Dynamics&Control,24:1747-1782.

Wang Y,Fang K-T.1994.Number theoretic methods intatistics[M].chapman&Hall,New York.

Wang X-Q,Fang K-T.2003.The effective dimension and quasi-Monte Carlo integration[J].Journal of Complexity,19(2):101-124.

Zhang T,Liu Y-B,Wu H-X,and Gu J-H.2009.Monte Carlo simulation of the exposure factor.China PhysicsB,18(06):2217-2222.

Monte CarloMethod and Quasi-Monte CarloMethod

ZHU Hui, L IU Yi-bao, YOU Yun
(East China Institute of Technology,Fuzhou,JX 344000,China)

Originated in the field of nuclear technology,Monte Carlo method is widely used in other related fields.As a non-deterministic numerical method,it can solve some unsolvable problems by deterministic algorithm.Development of computers,implementation of parallel algorithm,and application of quasi-Monte Carlo method compensateMonte carlo method’s deficiencies.W ith the introduction of random numbers,in this paper,Monte Carlo method and quasi-Monte Carlo method are introduced,including their basic principles,convergence,development,application prospects and problems.

pseudo-random numbers;quasi-random numbers;monte carlo;quasi-monte carlo;star descrepancy;halton sequence

O242.2

:A

:1674-3504(2010)04-357-06

10.3969/j.issn.1674-3504.2010.04.009

2010-01-11

江西省自然科学基金 (2009GZ W0001)

朱 辉 (1980—),男,硕士生,主要从事计算方法及应用研究。