蒙特卡羅方法與擬蒙特卡羅方法的歷史、現狀及展望

朱 輝, 劉義保, 游 運

(東華理工大學,江西撫州 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],上海:上?？茖W技術出版社.

許淑艷.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—),男,碩士生,主要從事計算方法及應用研究。