1、蒙特卡羅方法初探附 C 語言實例程序詳解2007-04-30 17:04蒙特卡羅( Monte Carlo )方法，又稱隨機抽樣或統計試驗方法，屬于計算數學的一個分支，它是在 本世紀四十年代中期為了適應當時原子能事業的發展而發展起來的。傳統的經驗方法由于不能逼近真實的物理過程，很難得到滿意的結果，而蒙特卡羅方法由于能 夠真實地模擬實際物理過程，故解決問題與實際非常符合，可以得到很圓滿的結果。這也是我們采用 該方法的原因。蒙特卡羅方法的基本原理及思想如下： 當所要求解的問題是某種事件出現的概率，或者是某個隨機變量的期望值時，它們可以通過某 種“試驗”的方法， 得到這種事件出現的頻率， 或者這個隨

2、機變數的平均值， 并用它們作為問題的解。 這就是蒙特卡羅方法的基本思想。蒙特卡羅方法通過抓住事物運動的幾何數量和幾何特征，利用數學 方法來加以模擬，即進行一種數字模擬實驗。它是以一個概率模型為基礎，按照這個模型所描繪的過 程，通過模擬實驗的結果，作為問題的近似解。可以把蒙特卡羅解題歸結為三個主要步驟：構造或描 述概率過程；實現從已知概率分布抽樣；建立各種估計量。蒙特卡羅解題三個主要步驟：構造或描述概率過程： 對于本身就具有隨機性質的問題，如粒子輸運問題，主要是正確描述和模擬這個概率過程，對 于本來不是隨機性質的確定性問題，比如計算定積分，就必須事先構造一個人為的概率過程，它的某 些參量正好是所

3、要求問題的解。即要將不具有隨機性質的問題轉化為隨機性質的問題。實現從已知概率分布抽樣： 構造了概率模型以后，由于各種概率模型都可以看作是由各種各樣的概率分布構成的，因此產 生已知概率分布的隨機變量(或隨機向量)，就成為實現蒙特卡羅方法模擬實驗的基本手段，這也是 蒙特卡羅方法被稱為隨機抽樣的原因。最簡單、最基本、最重要的一個概率分布是 (0,1) 上的均勻分 布(或稱矩形分布)。隨機數就是具有這種均勻分布的隨機變量。隨機數序列就是具有這種分布的總 體的一個簡單子樣，也就是一個具有這種分布的相互獨立的隨機變數序列。產生隨機數的問題，就是 從這個分布的抽樣問題。在計算機上，可以用物理方法產生隨機數，

4、但價格昂貴，不能重復，使用不 便。另一種方法是用數學遞推公式產生。這樣產生的序列，與真正的隨機數序列不同，所以稱為偽隨 機數，或偽隨機數序列。不過，經過多種統計檢驗表明，它與真正的隨機數，或隨機數序列具有相近 的性質，因此可把它作為真正的隨機數來使用。由已知分布隨機抽樣有各種方法，與從 (0,1) 上均勻 分布抽樣不同，這些方法都是借助于隨機序列來實現的，也就是說，都是以產生隨機數為前提的。由 此可見，隨機數是我們實現蒙特卡羅模擬的基本工具。建立各種估計量： 一般說來，構造了概率模型并能從中抽樣后，即實現模擬實驗后，我們就要確定一個隨機變量， 作為所要求的問題的解，我們稱它為無偏估計。建立各種

5、估計量，相當于對模擬實驗的結果進行考察 和登記，從中得到問題的解。蒙特卡羅方法與一般計算方法有很大區別，一般計算方法對于解決多維或因素復雜的問題非常 困難，而蒙特卡羅方法對于解決這方面的問題卻比較簡單。其特點如下：1 、直接追蹤粒子，物理思路清晰，易于理解。2 、采用隨機抽樣的方法，較真切的模擬粒子輸運的過程，反映了統計漲落的規律。3 、不受系統多維、多因素等復雜性的限制，是解決復雜系統粒子輸運問題的好方法。4 、MC程序結構清晰簡單。5 、研究人員采用 MC方法編寫程序來解決粒子輸運問題，比較容易得到自己想得到的任意中間 結果，應用靈活性強。6 、MC方法主要弱點是收斂速度較慢和誤差的概率性

6、質，其概率誤差正比于，如果單純以增大抽 樣粒子個數 N 來減小誤差，就要增加很大的計算量。近十年來，蒙特卡羅方法發展很快，從 1983 年到 1988 年期刊論文數量增長了五倍，有幾本好 書是關于電子 ? 光子蒙特卡羅問題的 注 1 ，蒙特卡羅方法的代碼被認為是黑匣子，它已成為計算數 學中不可缺少的組成部分，這主要是因為以下原因：1 、傳統的分析方法受到了問題復雜性的限制。2 、MC方法直觀，對實驗者很有吸引力。3 、計算機變得更快更便宜。4 、量子理論的發展為我們提供了輻射與物質相互作用的截面數據。 (以上內容來自網絡)簡而言之，蒙特卡羅方法就是用頻率來代替概率，當實驗樣本足夠大的時候，就可

7、以得到比較 精確的結果。下面就利用蒙特卡羅方法作一個簡單的任務， 從算法和編程中來理解蒙特卡羅方法的基本原理。 本人從接觸到蒙特卡羅方法到撰寫這篇日志剛好一天的時間，其中有很多誤解和不當之處，敬請各路 高手不吝賜教。題目：利用蒙特卡羅方法模擬計算圓周率PI 。分析：我們根據計算圓的面積公式來求解圓周率 PI ，當圓的半徑為 1 時，圓的面積剛好等于圓 周率 PI 。根據蒙特卡羅方法的思想，我們在以坐標原點為圓心作一個半徑為1 的單位圓，再作一個正方形與此圓相切(其實是否相切不重要，此處只是為了計算方便)。在這個正方形內任取M個點，并比較判斷是否落在圓內，將落在圓內的點計作 N,那么M與N的比值

8、就可以看成正方形和圓的面積之比。 即圓周率 PI=4*N/M 。下面是本題目的C語言源程序： /*/#include"stdio.h" #include"time.h" #include"stdlib.h" #include"graphics.h" main()int gd=DETECT,gm=0;long m=0,n=0,i; double xi,yi,y,y1; initgraph(&gd,&gm,""); /*設置圖形顯示模式 */setbkcolor(BLUE);/* 畫

9、坐標、正方形和圓 */line(200,50,200,400);line(200,50,205,60); line(200,50,195,60); line(50,200,400,200);line(400,200,390,195); line(400,200,390,205); rectangle(100,100,300,300);circle(200,200,100);/*/* 用時間作為產生隨機數的種子 */srand(unsigned)time(NULL);for(i=0;i<1000000;i+)/*設置產生隨機數的個數，這里為一百萬次 */* 產生-1 到 1之間的隨機數 x

10、i 和 yi*/xi=(rand()%(1000-0+1)+0)/500.0-1.0;yi=(rand()%(1000-0+1)+0)/500.0-1.0;putpixel(int)(xi*100+200),(int)(yi*100+200),(int)RED);/*在正方形區域內畫點*/*printf("xi=%f,yi=%fn",xi,yi);*/if(yi*yi>=(-(1-xi*xi)&&yi*yi<=1-xi*xi)/*判斷是否在圓內并計數*/m+,n+;else m+;y=4.0*n/m;/*計算 PI 值 */*outtextxg(

11、100,100,"PI=");*/printf("ntPI=%fn",y);getch();closegraph();/* */程序十分簡單，即產生隨機數、畫點、判別區域、統計，最后根據統計結果模擬計算PI值。我調試的結果如下：循環1000 次:3.2323.0603.0723.0803.22410000次:3.1143.1343.1373.1623.13610萬次：3.1453.1503.1513.1473.146100萬次：3.1433.1463.1453.1473.1471000 萬次：3.146 (時間長，只測了一次，呵呵)結果分析：1 、循環的

12、次數越多，結果越精確。2 、C語言rand()函數產生的并不是真正懂得隨機數，而是偽隨機數，會影響結果的精度。有興 趣的朋友不妨使用其它能產生真正意義的語言或算法試試。3 、算法過于簡單，計算機在處理數據時會產生很多的數據丟失，影響結果精度。此文作為學習交流之用，程序在TC2.0環境下調試通過。時間倉促，語焉不詳，望各路高手不吝賜教。#include<stdio.h>#include<stdlib.h>#include<math.h>#define L 40double funiseed (double iseed) 同余法產生隨機數 double mult

13、=1277,modulo,rmodulo;modulo=rmodulo=pow(2,17);iseed=iseed*mult; if(iseedvmodulo);else iseed=iseed-modulo;iseed=iseed/rmodulo;return iseed;main()int i,j,laL+1L+1,ipL+1,imL+1,count=0;int ici,ien,mcs,mcsmax,n0,mcstep;float w9,m,a;double jkt,iseed;/ 初始化數組for(i=0;i<=L;i+) for(j=0;j<=L;j+) laij = -1

14、;a=1./(L*L);/*for(i=1;i<=L;i+) ipi=i+1;imi=i-1;ipL=1;im1=L;/*printf("輸入樣本步數：”)； scanf("%d",&mcstep);printf("輸入 nO 的值：");scanf("%d",&n0);printf("輸入 jkt 的值：”)；scanf("%lf",&jkt);此方法用在2維Ising模型中/蒙特卡羅部分m=-(L*L);for(mcs=1;mcs<=mcstep+nO;mcs+) for(i=1;i<=L;i+) |for(j=1;jv=L;j+) ici=laij;ien=la ipi j+la imi j+lai ipj +lai imj;|ien=ici*ien;_|/ printf("put in the seed:"