Monte Carlo 方法（課堂PPT）

1、.1Monte Carlo 方法黃世萍.2起源起源 這一方法源于美國在第二次世界大戰(zhàn)進(jìn)研制原子彈的曼哈頓計(jì)劃。Monte Carlo方法創(chuàng)始人主要是這四位：Stanislaw Marcin Ulam, Enrico Fermi, John von Neumann和 Nicholas Metropolis。 .3 Monte Carlo 方法基礎(chǔ)Monte Carlo（MC）方法是在簡單的理論準(zhǔn)則基礎(chǔ)上（如簡單的物質(zhì)與物質(zhì)以及物質(zhì)與環(huán)境相互作用） ， 采用反復(fù)隨機(jī)抽樣的手段 ， 解決復(fù)雜系統(tǒng)的問題 。 該方法采用隨機(jī)抽樣的手法 ，可以模擬對象的概率與統(tǒng)計(jì)的問題。通過設(shè)計(jì)適當(dāng)?shù)母怕誓Ｐ?， 該方法

2、還可以解決確定性問題 ， 如定積分等。 隨著計(jì)算機(jī)的迅速發(fā)展 ， MC 方法已在應(yīng)用物理、原子能、固體物理 、化學(xué)、材料、生物、生態(tài)學(xué)、社會(huì)學(xué)以及經(jīng)濟(jì)學(xué)等領(lǐng)域得到了廣泛的應(yīng)用。.4MC meothds classical Monte Carlo: samples are drawn from a probability distribution, often the classical Boltzmann distribution, to obtain thermodynamic properties or minimum-energy structures;quantum Monte Car

3、lo: random walks are used to compute quantum-mechanical energies and wavefunctions, often to solve electronic structure problems, using Schrdingers equation as a formal starting point; volumetric Monte Carlo: random number generators are used to generate volumes per atom or to perform other types of

4、 geometrical analysis;kinetic Monte Carlo: simulate processes using scaling arguments to establish timescales or by introducing stochastic effects into molecular dynamics.5Monte Carlo方法的基本思想Monte Carlo 方法的基本思想是: 為了求解某個(gè)問題 ， 建立一個(gè)恰當(dāng)?shù)母怕誓Ｐ突螂S機(jī)過程 ， 使得其參量（如事件的概率、隨機(jī)變量的數(shù)學(xué)期望等）等于所求問題的解 ， 然后對模型或過程進(jìn)行反復(fù)多次的隨機(jī)抽樣試驗(yàn) ，

5、 并對結(jié)果進(jìn)行統(tǒng)計(jì)分析 ， 最后計(jì)算所求參量 ， 得到問題的近似解。.6Monte Carlo 方法是隨機(jī)模擬方法 ； 它不僅限于模擬隨機(jī)性問題 ， 還可以解決確定性的數(shù)學(xué)問題。 對隨機(jī)性問題 ， 可以根據(jù)實(shí)際問題的概率法則 ， 直接進(jìn)行隨機(jī)抽樣試驗(yàn) ， 即直接模擬方法。 對于確定性問題采用間接模擬方法 ， 即通過統(tǒng)計(jì)分析隨機(jī)抽樣的結(jié)果獲得確定性問題的解。用 Monte Carlo 方法解決確定性的問題主要是在數(shù)學(xué)領(lǐng)域 ， 如計(jì)算重積分、求逆矩陣、解線性代數(shù)方程組、解積分方程、解偏微分方程邊界問題和計(jì)算微分算子的特征值等。 用 Monte Carlo 方法解決隨機(jī)性問題則在眾多的科學(xué)及應(yīng)用技術(shù)

6、領(lǐng)域得到廣泛的應(yīng)用 ， 如中子在介質(zhì)中的擴(kuò)散問題、庫存問題、隨機(jī)服務(wù)系統(tǒng)中的排隊(duì)問題、動(dòng)物的生態(tài)競爭、傳染病的蔓延等。.7簡單的例子對積分進(jìn)行變換 ， 構(gòu)造新的被積函數(shù) g（x） ， 使得該函數(shù)滿足下列條件 ：g（x）是連續(xù)隨機(jī)變量 的概率密度函數(shù)。 定積分是概率積分 ， 其積分值等于概率 Pr （a b） ， 即這個(gè)步驟就是將一個(gè)積分轉(zhuǎn)化為一個(gè)概率模型的過程 ； 然后 ， 反復(fù)多次的隨機(jī)抽樣試驗(yàn) ，以抽樣結(jié)果的統(tǒng)計(jì)平均作為索求概率的近似值 ， 從而求得該積分。.8具體試驗(yàn)步驟如下 ：（1） 產(chǎn)生服從給定分布函數(shù) g（x）的隨機(jī)變量值 xi（2）檢查 xi是否落入積分區(qū)域（a x b） ， 如

7、果滿足條件 ， 則記錄一次。反復(fù)進(jìn)行上述試驗(yàn)。 假設(shè)在 N 次試驗(yàn)后 ， xi落入積分區(qū)域的總次數(shù)為 m ， 那么 ， 積分值近似表示為對于隨機(jī)性問題 ， 可直接將實(shí)際的隨機(jī)問題抽象為概率數(shù)學(xué)模型 ， 然后與求解確定性問題一樣進(jìn)行抽樣試驗(yàn)和統(tǒng)計(jì)計(jì)算。.9Monte Carlo 方法解決實(shí)際問題的過程中 ， 主要有以下幾個(gè)內(nèi)容 建立簡單而又便于實(shí)現(xiàn)的概率統(tǒng)計(jì)模型 ， 使所求的解正是該模型的某 一事件的概率或數(shù)學(xué)期望 ， 或該模型能夠直接描述實(shí)際的隨機(jī)過程。 根據(jù)概率統(tǒng)計(jì)模型的特點(diǎn)和計(jì)算的需求 ， 改進(jìn)模型 ， 以便減小方差和減低費(fèi)用 ， 提高計(jì)算效率。 建立隨機(jī)變量的抽樣方法 ， 包括偽隨機(jī)數(shù)和

8、服從特定分布的隨機(jī)變量的產(chǎn)生方法。 給出統(tǒng)計(jì)估計(jì)值及其方差或標(biāo)準(zhǔn)誤差。.10One-Dimensional IntegralsnMethodical approachesrectangle rule, trapezoid rule, Simpsons rule nQuadrature formulan uniformly separated pointsSum areas of shapes approximating shape of curve( )baIf x dxEvaluating the general integral( )f xx11()()nniiiibaIxf xf xn

9、.11假設(shè)對于下列積分 ( )baIf x dx首先，按照均勻分布，隨即采樣區(qū)域a,b，記采樣點(diǎn)為Xi，每個(gè)采樣點(diǎn)的值p(x)=1/(b-a)。利用Monte Carlo算法，上述積分可以模擬為： 當(dāng)N趨近于無窮大的時(shí)候，我們認(rèn)為兩者的值是一致的，所以可以利用上述離散的方式模擬這個(gè)積分 .12Monte Carlo IntegrationnStochastic approachnSame quadrature formula, different selection of points1()niibaIf xnn points selected from uniform distributio

10、n p(x)xp(x).13我們在對于a,b進(jìn)行采樣的時(shí)候，完全沒有必要進(jìn)行均勻采樣，這樣做只是簡單而已，我們可以根據(jù)一種概率分布來進(jìn)行采樣，從而使得某些區(qū)域的采樣密度更大。加入采樣點(diǎn)的概率分布函數(shù)為P(x)，那么我們可以利用如下公式計(jì)算定積分的值： 證明：.143 Monte Carlo方法的收斂性和基本特點(diǎn)設(shè)所求的量 x是隨機(jī)變量的數(shù)學(xué)期望 E（x） ， 那么 ， Monte Carlo 方法通常使用隨機(jī)變量的簡單子樣， ， ，N 的算術(shù)平均值 ， 即作為所求量 x 的近似值。 由柯爾莫哥羅夫（Kolmogorov）大數(shù)定理可知 ，即當(dāng) N 充分大時(shí) ， 有成立的概率等于 ， 亦即可以用

11、N 作為所求量 x 的估計(jì)值。.15根據(jù)中心極限定理 ， 如果隨機(jī)變量 的標(biāo)準(zhǔn)差 不為零 ， 那么 Monte Carlo 方法的誤差為式中 ， 為正態(tài)差 ， 是與置信水平有關(guān)的常量。 Monte Carlo 方法的收斂速度的階為 o（N -1/2） ， 誤差是由隨機(jī)變量的標(biāo)準(zhǔn)差 s和抽樣次數(shù) N 決定的。提高精度一位數(shù) ， 抽樣次數(shù)要增加 倍 ； 減小隨機(jī)變量的標(biāo)準(zhǔn)差 ， 可以減小誤差 。.16Monte Carlo 方法具有以下四個(gè)重要特征 ： 由于 Monte Carlo 方法是通過大量簡單的重復(fù)抽樣來實(shí)現(xiàn)的 ， 因此 ， 方法和程序的結(jié)構(gòu)十分簡單 。 收斂速度比較慢 ， 因此 ， 較適

12、用于求解精度要求不高的問題。 收斂速度與問題的維數(shù)無關(guān) ， 因此 ， 較適用于求解多維問題。 問題的求解過程取決于所構(gòu)造的概率模型 ， 而受問題條件限制的影響較小 ， 因此 ， 對各種問題的適應(yīng)性很強(qiáng)。.17隨機(jī)數(shù)的產(chǎn)生1 隨機(jī)數(shù)與偽隨機(jī)數(shù)Monte Carlo 方法的核心是隨機(jī)抽樣。 在該過程中往往需要各種各樣分布的隨機(jī)變量其中最簡單、最基本的是在 ，區(qū)間上均勻分布的隨機(jī)變量。 在該隨機(jī)變量總體中抽取的子樣 ， ， ，N 稱為隨機(jī)數(shù)序列 ， 其中每個(gè)個(gè)體稱為隨機(jī)數(shù)。用數(shù)學(xué)的方法產(chǎn)生隨機(jī)數(shù)是目前廣泛使用的方法。 該方法的基本思想是利用一種遞推公式 ：對于給定的初始值 ， 逐個(gè)地產(chǎn)生 ， ， .

13、18數(shù)學(xué)方法產(chǎn)生的隨機(jī)數(shù)存在兩個(gè)問題 ： 整個(gè)隨機(jī)數(shù)序列是完全由遞推函數(shù)形式和初始值唯一確定的 ， 嚴(yán)格地說不滿足隨機(jī)數(shù)的相互獨(dú)立的要求。 存在周期現(xiàn)象?；谶@兩個(gè)原因 ， 將用數(shù)學(xué)方法所產(chǎn)生的隨機(jī)數(shù)稱為偽隨機(jī)數(shù)。 偽隨機(jī)數(shù)的優(yōu)點(diǎn)是適用于計(jì)算機(jī) ， 產(chǎn)生速度快 ， 費(fèi)用低廉。 目前 ， 多數(shù)計(jì)算機(jī)均附帶有“隨機(jī)數(shù)發(fā)生器” 。.19選擇遞推函數(shù)必須注意以下幾點(diǎn) ： 隨機(jī)性好 ； 在計(jì)算機(jī)上容易實(shí)現(xiàn) ； 省時(shí) ； 偽隨機(jī)數(shù)的周期長。.202 偽隨機(jī)數(shù)的產(chǎn)生方法最基本的偽隨機(jī)數(shù)是均勻分布的偽隨機(jī)數(shù)。該方法是首先給一個(gè) r位的數(shù) ， 取其中間的 r位數(shù)碼作為第一個(gè)偽隨機(jī)數(shù) ， 然后將這個(gè)數(shù)平方 ， 構(gòu)

14、成一個(gè)新的 r位的數(shù) ， 再取中間的 r位數(shù)作為第二個(gè)偽隨機(jī)數(shù)。 如此循環(huán)可得到一個(gè)偽隨機(jī)數(shù)序列。 該方法的遞推公式為x表示對 x 取整 ， 運(yùn)算 B（Mod M）表示 B被 M 整除后的余數(shù)。 數(shù)列i 是分布在 ，上的。 該方法由于效率較低 ， 有時(shí)周期較短 ， 甚至?xí)霈F(xiàn)零 。.21同余法該方法也是由選定的初始值出發(fā) ， 通過遞推產(chǎn)生偽隨機(jī)數(shù)序列。由于該遞推公式可寫成數(shù)論中的同余式 ， 故稱同余法。 該方法的遞推公式為a ， c ， m分別稱作倍數(shù)（multiplier） 、增值（Increment ）和模（modulus） ， 均為正整數(shù) 。 x 稱為種子或初值 ， 也為正整數(shù)。 該方法

15、所產(chǎn)生偽隨機(jī)數(shù)的質(zhì)量 ， 如周期的長度、獨(dú)立性和均勻性都與式中三個(gè)參數(shù)有關(guān)。.22.233 偽隨機(jī)數(shù)的統(tǒng)計(jì)檢驗(yàn)偽隨機(jī)數(shù)的特性好壞將直接影響到 Monte Carlo 的計(jì)算結(jié)果 ， 因此要對所產(chǎn)生的偽隨機(jī)數(shù)序列進(jìn)行隨機(jī)性檢驗(yàn)。 隨機(jī)性檢驗(yàn)主要包括均勻性檢驗(yàn) 、獨(dú)立性檢驗(yàn)、組合規(guī)律性檢驗(yàn)和無連貫性檢驗(yàn)。 檢驗(yàn)是偽隨機(jī)數(shù)檢驗(yàn)最常用的方法。1. 均勻性就是偽隨機(jī)數(shù)列的 N 個(gè)數(shù)是否均勻分布在 ，區(qū)間上。 若將 ，區(qū)間分成 k個(gè)相等的子區(qū)間（一般 k ， ， ） ， 若所得偽隨機(jī)數(shù)在 ，區(qū)間上是均勻分布的 ， 則虛假設(shè)H 應(yīng)為“每個(gè)偽隨機(jī)數(shù)屬于第 i組的概率為.24頻率檢驗(yàn) 檢驗(yàn)每組觀測頻數(shù) ni與理

16、論頻數(shù)mi N /k之間相差的顯著性3. 獨(dú)立性按先后順序排列的 N 個(gè)偽隨機(jī)數(shù)中 ， 每個(gè)數(shù)的出現(xiàn)是否與其前后各個(gè)數(shù)獨(dú)立無關(guān)。 對于兩組偽隨機(jī)數(shù)來說 ， 獨(dú)立性就是指它們不相關(guān) 。4. 組合規(guī)律性將 N 個(gè)偽隨機(jī)數(shù)按一定的規(guī)律組合起來 ， 則各種組合的出現(xiàn)具有一定的概率。5. 無連貫性將一次出現(xiàn)的 N 個(gè)偽隨機(jī)數(shù) ， 按其大小分為兩類或 k類 ， 則各類數(shù)的出現(xiàn)是否沒有連貫現(xiàn)象。.25確定性問題的 Monte Carlo 方法求解Monte Carlo 方法所能解決的問題可以分為兩大類確定性問題隨機(jī)性問題確定性問題主要包括： 求解線性和非線性方程組、逆矩陣、橢圓型差分方程的邊值、積分方程以及

17、多重積分等。用 Monte Carlo 方法求解確定性問題的基本思想是 ： 對于給定的確定性問題 ， 設(shè)計(jì)一個(gè)概率統(tǒng)計(jì)模型（或隨機(jī)過程模型） ， 然后采用一定的抽樣方法 ， 按照所設(shè)計(jì)的概率統(tǒng)計(jì)模型進(jìn)行抽樣 ， 最后把這個(gè)模型產(chǎn)生的一個(gè)數(shù)字特征作為該確定性問題的近似解。.26蒲豐試驗(yàn) 年法國著名學(xué)者蒲豐（Buffon）就發(fā)表了采用隨機(jī)抽樣法計(jì)算圓周率的論文。 這就是蒲豐隨機(jī)投針試驗(yàn) ， 即著名的蒲豐問題。蒲豐概率模型是 : 在平面上畫相距均為a的平行線束 ， 向該平面上隨機(jī)投置一枚長為 l 的針。 為了避免針同時(shí)與兩條平行線相交的復(fù)雜情況 ， 設(shè)定 l a。 如圖所示 ， M 為針的中點(diǎn) ，

18、y為針的中點(diǎn)M 到與之最近的平行線的距離 ， 為針與平行線的夾角。 顯然 ， y a ， 。 該隨機(jī)試驗(yàn)所有可能的集為.27針與平行線相交的充分必要條件是 y lsin ，針與平行線相交事件的集為則針與平行線相交的概率為由于 l和 a均為已知常數(shù) ， 只要通過大量抽樣試驗(yàn)求得該概率 p ， 由上式即可算出圓周率。 設(shè)投針總次數(shù)為 N ， 其中針與平行線相交次數(shù)為 v ， 由貝努里（Bernoulli）定理可知 ， 當(dāng) N充分大時(shí) ， 該事件出現(xiàn)的頻率接近于其概率 ， 即這就是蒲豐試驗(yàn)求圓周率的過程。雖然 ， 該方法要想獲得較高的精確度 ， 需要數(shù)以百萬次的抽樣試驗(yàn) ， 效率很低 ， 但蒲豐試驗(yàn)

19、具有 Monte Carlo 方法解決確定性問題的基本思想。.28A Classic Example The Calculation of .29.30隨機(jī)性問題的 Monte Carlo 模擬該過程是先建立一個(gè)隨機(jī)過程模型 ， 使得該過程的隨機(jī)變量的數(shù)學(xué)期望等于所要求解確定問題的解。這種計(jì)算方法叫做間接模擬方法。 另一方面 ， 采用隨機(jī)試驗(yàn)的方法直接模擬隨機(jī)過程 ， 解決隨機(jī)問題 ， 這就是所謂直接 Monte Carlo 模擬。 如模擬布朗運(yùn)動(dòng)、擴(kuò)散過程、有機(jī)高分子形態(tài)、晶粒生長等隨機(jī)問題的模擬。.311 隨機(jī)行走（random walk）模擬隨機(jī)行走是一種典型的簡單抽樣方法 ， 可用以模

20、擬擴(kuò)散 、溶液中長而柔性的大分子的性質(zhì)等.隨機(jī)行走主要有三種類型:無限制的隨機(jī)行走（Unrestricted random walk ， RW）不退行走（Nonreversal random walk ， NRRW）自回避行走（Self-avoiding walk ， SAW）.32無限制隨機(jī)行走就是指 ， 某一個(gè)質(zhì)點(diǎn)的每一次行走沒有任何限制 ， 既與前一次行走無關(guān) ， 也與以前任何一步所到的位置無關(guān)。 這種模型可以用于模擬質(zhì)點(diǎn)的擴(kuò)散等過程 ， 但是 ， 不能用于模擬高分子的位形。 因?yàn)?， 用隨機(jī)行走方法模擬高分子位形是用隨機(jī)行走的軌跡代表高分子的位形 ， 行走過的位置代表的是構(gòu)成分子的原子

21、或官能團(tuán) ， 因此 ， 無限制隨機(jī)行走忽略了體斥效應(yīng)。不退行走就是禁止在每一步行走后立即倒退 ， 可以解決剛走的一步與上一步重疊的問題。 但不退行走沒有完全解決高分子的體斥效應(yīng)問題。自回避行走就是所有已走過的位置不能再走 ， 這樣就完全解決了體斥效應(yīng)問題。.33二維方格子上的三種隨機(jī)行走的示意圖。 可以用四個(gè)矢量記述從某個(gè)節(jié)點(diǎn)向鄰近節(jié)點(diǎn)的行走方向（設(shè)格子間距為） ：.34N 步無限制隨機(jī)行走的算法如下 ： 取 r （ ， ）（坐標(biāo)原點(diǎn)） ， 并令 k ； 取一個(gè)在 和 之間的隨機(jī)整數(shù) vk ； k k ， rk rk （vk ） ； 若 k N ， 行走結(jié)束 ， 否則回到第 步。對于不退行走

22、， 可選擇的行走方向不再是 ， 而是 。 禁止在每一步行走后立即倒退 ， 即第k步的方向矢量不能與第 k 步的方向矢量相逆。 由式（ ）可以看出方向矢量 （） ， （）互為逆方向 ， （） ， （）互為逆方向。.352 Markov鏈對于簡單抽樣 ， 每一次的抽樣都是獨(dú)立的。如上述的隨機(jī)行走過程 ， 每行走一步都與前一步無關(guān) ， 更與初始位置無關(guān)。Sampling.36Integrate Over a Simple Shape? nStatistical-mechanics integrals typically have significant contributions from mini

23、scule regions of the integration space contributions come only when no spheres overlap e.g., 100 spheres at freezing the fraction is 10-260 nEvaluation of integral is possible only if restricted to region important to integralmust contend with complex shape of regionMC methods highly suited to “impo

24、rtance sampling”()11!()NNNNU rZNUdr U re0Ue.37Importance SamplingnPut more quadrature points in regions where integral receives its greatest contributionsnReturn to 1-dimensional examplenMost contribution from region near x = 1nChoose quadrature pointsnot uniformly, but accordingto distribution (x)l

25、inear form is one possibilitynHow to revise the integral to remove the bias?1203Ix dx2( )3f xx( )2xx.38The Importance-Sampled IntegralnConsider a rectangle-rule quadrature with unevenly spaced abscissasnSpacing between pointsreciprocal of local number of points per unit lengthnImportance-sampled rec

26、tangle ruleSame formula for MC sampling1()niiiIf xx1x2x3xnx1()iibaxnxGreater more points smaller spacing1( )()()niiixbaf xInxchoose x points according to .39在正則系綜中 ， 任意觀察量 A（x）的熱平均為.40重要性抽樣：重要性抽樣：在隨機(jī)過程中，選取一個(gè)隨機(jī)抽樣的分布，使生在隨機(jī)過程中，選取一個(gè)隨機(jī)抽樣的分布，使生成的隨機(jī)數(shù)滿足選取分布形式。成的隨機(jī)數(shù)滿足選取分布形式。根據(jù)一定的分布根據(jù)一定的分布形式進(jìn)行的隨機(jī)抽樣稱為重要性抽樣。形式進(jìn)

27、行的隨機(jī)抽樣稱為重要性抽樣。.41重要抽樣 Monte Carlo 方法的實(shí)質(zhì)是每次抽樣試驗(yàn)不是完全獨(dú)立的 ， 而是與前一次或者與以前的所有抽樣結(jié)果具有一定的概率關(guān)系 ， 如不退隨機(jī)行走和自回避隨機(jī)行走。設(shè)一個(gè)系統(tǒng)的狀態(tài)序列（隨機(jī)變量序列）為 x ，x ， ，xn ， ， 如果對于任何一個(gè)狀態(tài) xn只與前一個(gè)狀態(tài) xn 有關(guān) ， 而與初始狀態(tài)無關(guān) ， 即狀態(tài) n的概率為則稱此序列為 Markov 鏈。.42Markov 鏈?zhǔn)且环N隨機(jī)行走狀態(tài) ， 從狀態(tài) i單步行走到狀態(tài) j 的概率叫做轉(zhuǎn)移概率 ， 或躍遷概率 ， 即設(shè)所有可能的狀態(tài)數(shù)為 N ， 由 pij 構(gòu)成的 N N 矩陣叫做轉(zhuǎn)移矩陣 p

28、 ， 該矩陣的每一行元素的和等于 。 Markov 鏈的重要性質(zhì)是 ， 無論初始狀態(tài)如何 ， 最終狀態(tài)（足夠多的時(shí)間步長次數(shù)）會(huì)遵從某一個(gè)唯一的分布 ， 該分布叫做極限分布 xlim ， 即也就是說 ， 極限狀態(tài)乘轉(zhuǎn)移概率后狀態(tài)不再發(fā)生變化 ， 即系統(tǒng)達(dá)到一個(gè)平衡狀態(tài)。 因此 ，Markov 鏈在平衡態(tài) Monte Carlo 模擬中具有重要的意義。.433 Metropolis Monte Carlo法Monte Carlo 方法主要分為簡單隨機(jī)抽樣方法和重要隨機(jī)抽樣方法。 簡單抽樣就是以平均分布進(jìn)行抽樣 ， 每次抽樣是完全獨(dú)立的。正如前面關(guān)于積分問題中所述 ， 很多問題難以用簡單抽樣方法解決