統(tǒng)計模擬介紹

1、專業(yè)必修課教學目的：掌握統(tǒng)計模擬在實際中的應用36理論課+18實驗課(PBL)考試形式：閉卷考試(60%)+小論文(40%)1.統(tǒng)計模擬，Ross著，王兆軍，陳廣雷，鄒長亮譯 2007年7月由人民郵電出版社出版2.統(tǒng)計建模與R軟件，薛毅,陳廣萍編著，2007年4月清華大學出版社Monte CarloMonte Carlo方法：方法：亦稱統(tǒng)計模擬方法，亦稱統(tǒng)計模擬方法，statistical simulation methodstatistical simulation method 利用隨機數(shù)進行數(shù)值模擬的方法利用隨機數(shù)進行數(shù)值模擬的方法Monte CarloMonte Carlo名字的由來：

2、名字的由來： 是由是由MetropolisMetropolis在二次世界大戰(zhàn)期間提出的：在二次世界大戰(zhàn)期間提出的：ManhattanManhattan計劃，研究與原子計劃，研究與原子彈有關的中子輸運過程；彈有關的中子輸運過程； Monte CarloMonte Carlo是摩納哥（是摩納哥（monaco)monaco)的首都，該城以賭博聞名的首都，該城以賭博聞名Nicholas Metropolis (1915-1999)Monte-Carlo, MonacoMonte CarloMonte Carlo方法簡史方法簡史簡單地介紹一下簡單地介紹一下Monte CarloMonte Carlo方法

3、的發(fā)展歷史方法的發(fā)展歷史1 1、BuffonBuffon投針實驗：投針實驗：17681768年，法國數(shù)學家年，法國數(shù)學家Comte de Buffon利用投針實驗估計利用投針實驗估計 的的值值dLp2dLProblem of Buffons needle:Problem of Buffons needle:If a needle of lengthIf a needle of length l l is dropped at random on the is dropped at random on the middle of a horizontal surface ruled with p

4、arallel middle of a horizontal surface ruled with parallel lines a distance lines a distance d d l l apart, what is the probability apart, what is the probability that the needle will cross one of the lines?that the needle will cross one of the lines?SolutionSolution:The positioning of the needle Th

5、e positioning of the needle relative to nearby lines can be relative to nearby lines can be described with a random vector described with a random vector which has components:which has components:),0),0dAThe random vector is uniformly distributed on the region The random vector is uniformly distribu

6、ted on the region 0,0,d d) )0,0, ). Accordingly, it has probability density function 1/d). Accordingly, it has probability density function 1/d . .The probability that the needle will cross one of the lines is given The probability that the needle will cross one of the lines is given by the integral

7、by the integraldldAdpld20sin01 1777年，古稀之年的蒲豐在家中請來好些客人玩投針游戲（針長是線距之半），他事先沒有給客人講與有關的事?？腿藗冸m然不知道主人的用意，但是都參加了游戲。他們共投針2212次，其中704次相交。蒲豐說，2212/704=3.142，這就是值。這著實讓人們驚喜不已。2 2、19301930年，年，Enrico FermiEnrico Fermi利用利用Monte CarloMonte Carlo方法研究中子方法研究中子的擴散，并設計了一個的擴散，并設計了一個Monte CarloMonte Carlo機械裝置，機械裝置，F(xiàn)ermiac,F

8、ermiac,用于計算核反應堆的臨界狀態(tài)用于計算核反應堆的臨界狀態(tài)3 3、Von NeumannVon Neumann是是Monte CarloMonte Carlo方法的正式奠基者方法的正式奠基者, ,他與他與Stanislaw UlamStanislaw Ulam合作建立了概率密度函數(shù)、反累積分布函數(shù)合作建立了概率密度函數(shù)、反累積分布函數(shù)的數(shù)學基礎，以及偽隨機數(shù)產(chǎn)生器。在這些工作中，的數(shù)學基礎，以及偽隨機數(shù)產(chǎn)生器。在這些工作中， Stanislaw UlamStanislaw Ulam意識到了數(shù)字計算機的重要性意識到了數(shù)字計算機的重要性合作起源于合作起源于ManhattanManhatta

9、n工程：利用工程：利用ENIAC(Electronic ENIAC(Electronic Numerical Integrator and Computer)Numerical Integrator and Computer)計算產(chǎn)額計算產(chǎn)額一些人進行了實驗，其結果列于下表 ：實驗者年份投計次數(shù)的實驗值沃爾弗(Wolf)185050003.1596斯密思(Smith)185532043.1553?？怂?Fox)189411203.1419拉查里尼(Lazzarini)190134083.1415929蒙特卡羅方法又稱計算機隨機模擬方法。它是以概率統(tǒng)計理論為基礎的一種方法。 由蒲豐試驗可以看出，

10、當所求問題的解是某個事件的概率，或者是某個隨機變量的數(shù)學期望，或者是與概率、數(shù)學期望有關的量時，通過某種試驗的方法，得出該事件發(fā)生的頻率，或者該隨機變量若干個具體觀察值的算術平均值，通過它得到問題的解。這就是蒙特卡羅方法的基本思想。 l=1;d=2;m=0;n=10000for k=1:n;x=unifrnd(0,d);y=unifrnd(0,pi);if x1*sin(y)m=m+1elseendendp=m/npi_m=1/p關系式成立產(chǎn)生隨機數(shù)驗證模型成立次數(shù)k=k+1否是計算估計結果k/n成立次數(shù)不變試驗次數(shù)是否達到n次是否編寫R程序Monte CarloMonte Carlo模擬的應

11、用：模擬的應用：自然現(xiàn)象的模擬：自然現(xiàn)象的模擬：宇宙射線在地球大氣中的傳輸過程；宇宙射線在地球大氣中的傳輸過程；高能物理實驗中的核相互作用過程；高能物理實驗中的核相互作用過程；實驗探測器的模擬實驗探測器的模擬數(shù)值分析：數(shù)值分析：利用利用Monte CarloMonte Carlo方法求積分方法求積分金融工程：金融工程：股票期權的模擬定價股票期權的模擬定價物理 化學 生物 環(huán)境工程 醫(yī)學 金融 交通教育 心里 衛(wèi)生 數(shù)學語言 軍事 歷史 經(jīng)濟天文。Monte CarloMonte Carlo模擬在物理研究中的作用模擬在物理研究中的作用Monte CarloMonte Carlo模擬的步驟：模擬的

12、步驟：1.1. 根據(jù)欲研究的物理系統(tǒng)的性質，建立能夠描述該系統(tǒng)特性根據(jù)欲研究的物理系統(tǒng)的性質，建立能夠描述該系統(tǒng)特性的理論模型，導出該模型的某些特征量的概率密度函數(shù)；的理論模型，導出該模型的某些特征量的概率密度函數(shù)；2.2. 從概率密度函數(shù)出發(fā)進行隨機抽樣，得到特征量的一些模從概率密度函數(shù)出發(fā)進行隨機抽樣，得到特征量的一些模擬結果；擬結果；3.3. 對模擬結果進行分析總結，預言物理系統(tǒng)的某些特性。對模擬結果進行分析總結，預言物理系統(tǒng)的某些特性。注意以下兩點：注意以下兩點：Monte CarloMonte Carlo方法與數(shù)值解法的不同方法與數(shù)值解法的不同: :Monte CarloMonte

13、Carlo方法利用隨機抽樣的方法來求解物理問方法利用隨機抽樣的方法來求解物理問題題; ;數(shù)值解法數(shù)值解法: :從一個物理系統(tǒng)的數(shù)學模型出發(fā)從一個物理系統(tǒng)的數(shù)學模型出發(fā), ,通過求解一通過求解一系列的微分方程來的導出系統(tǒng)的未知狀態(tài)系列的微分方程來的導出系統(tǒng)的未知狀態(tài); ;Monte CarloMonte Carlo方法并非只能用來解決包含隨機的過程的問題方法并非只能用來解決包含隨機的過程的問題: :許多利用許多利用Monte CarloMonte Carlo方法進行求解的問題中并不包含方法進行求解的問題中并不包含隨機過程隨機過程 例如例如: :用用Monte CarloMonte Carlo方法

14、計算定積分方法計算定積分. . 對這樣的問題可將其轉換成相關的隨機過程對這樣的問題可將其轉換成相關的隨機過程, , 然后用然后用Monte CarloMonte Carlo方法進行求解方法進行求解Monte CarloMonte Carlo算法的主要組成部分算法的主要組成部分概率密度函數(shù)概率密度函數(shù)(pdf)(pdf) 必須給出描述一個物理系統(tǒng)的一組概必須給出描述一個物理系統(tǒng)的一組概率密度函數(shù)率密度函數(shù); ;隨機數(shù)產(chǎn)生器隨機數(shù)產(chǎn)生器能夠產(chǎn)生在區(qū)間能夠產(chǎn)生在區(qū)間0,10,1上均勻分布的隨機數(shù)上均勻分布的隨機數(shù)抽樣規(guī)則抽樣規(guī)則如何從在區(qū)間如何從在區(qū)間0,10,1上均勻分布的隨機數(shù)出發(fā)上均勻分布的隨

15、機數(shù)出發(fā), ,隨隨機抽取服從給定的機抽取服從給定的pdfpdf的隨機變量的隨機變量; ;模擬結果記錄模擬結果記錄記錄一些感興趣的量的模擬結果記錄一些感興趣的量的模擬結果誤差估計誤差估計必須確定統(tǒng)計誤差（或方差）隨模擬次數(shù)以及必須確定統(tǒng)計誤差（或方差）隨模擬次數(shù)以及其它一些量的變化；其它一些量的變化；減少方差的技術減少方差的技術利用該技術可減少模擬過程中計算的次數(shù)；利用該技術可減少模擬過程中計算的次數(shù)；并行和矢量化并行和矢量化可以在先進的并行計算機上運行的有效算法可以在先進的并行計算機上運行的有效算法1.具有統(tǒng)計功能的軟件 Excel, Matab, C, Fortran 2. 專業(yè)的統(tǒng)計軟件

16、SPSS, SAS, S-Plus, R, Gauss, minitab優(yōu)點1) 能夠比較逼真地描述具有隨機性質的事物的特點及物理實驗過程。2) 受幾何條件限制小。3) 收斂速度與問題的維數(shù)無關。4) 具有同時計算多個方案與多個未知量的能力。5) 誤差容易確定。6) 程序結構簡單，易于實現(xiàn)。 缺點1) 收斂速度慢。2) 誤差具有概率性。3) 在粒子輸運問題中，計算結果與系統(tǒng)大小有關?；疖囯x站時刻火車離站時刻13:0013:0513:10概率概率0.70.20.1 一列列車從一列列車從A A站開往站開往B B站，某人每天趕往站，某人每天趕往B B站上車。他站上車。他已經(jīng)了解到火車從已經(jīng)了解到火車

17、從A A站到站到B B站的運行時間是服從均值為站的運行時間是服從均值為30min30min，標準差為，標準差為2min2min的正態(tài)隨機變量?；疖嚧蠹s下午的正態(tài)隨機變量。火車大約下午1313：0000離開離開A A站，此人大約站，此人大約13:3013:30到達到達B B站?；疖囯x開站?；疖囯x開A A站的時刻站的時刻及概率如表及概率如表1 1所示，此人到達所示，此人到達B B站的時刻及概率如表站的時刻及概率如表2 2所示。所示。問此人能趕上火車的概率有多大？問此人能趕上火車的概率有多大？表1：火車離開A站的時刻及概率 表2：某人到達B站的時刻及概率 人到站時刻人到站時刻13:2813:3013

18、:3213:34概率概率0.30.40.20.1這個問題用概率論的方法求解十分困難，它涉及此人到達時刻、火車離開站的時刻、火車運行時間幾個隨機變量，而且火車運行時間是服從正態(tài)分布的隨機變量，沒有有效的解析方法來進行概率計算。在這種情況下可以用計算機模擬的方法來解決。：火車從：火車從A A站出發(fā)的時刻；站出發(fā)的時刻；：火車從：火車從A A站到站到B B站的運行時間；站的運行時間；：某人到達：某人到達B B站的時刻；站的時刻；：隨機變量：隨機變量 服從正態(tài)分布的均值；服從正態(tài)分布的均值；：隨機變量：隨機變量 服從正態(tài)分布的標準差服從正態(tài)分布的標準差；1T2T3T2T2T此人能及時趕上火車的充分必要

19、條件為：此人能及時趕上火車的充分必要條件為： ，所以此人能趕上火車的概率模型為：所以此人能趕上火車的概率模型為： 。123TTT123p TTT為了分析簡化，假定為了分析簡化，假定1313時為時刻時為時刻t=0，則變量，則變量 、 的分布律為：的分布律為：1T3T05100.70.20.1283032340.30.40.20.11/minT( )P t3/minT( )P t關系式成立產(chǎn)生隨機數(shù)驗證模型成立次數(shù)k=k+1否是計算估計結果k/n成立次數(shù)不變試驗次數(shù)是否達到n次是否編寫R程序借助區(qū)間借助區(qū)間(0,1)(0,1)分布產(chǎn)生的隨機數(shù)，對分布產(chǎn)生的隨機數(shù)，對變量變量 、 概率分布進行統(tǒng)計模

20、擬；概率分布進行統(tǒng)計模擬；1T3T根據(jù)變量根據(jù)變量 、 、 概率分布及模擬概率分布及模擬程序、命令產(chǎn)生程序、命令產(chǎn)生n 個隨機分布數(shù)；個隨機分布數(shù)；1T2T3T使用隨機產(chǎn)生的使用隨機產(chǎn)生的n 組隨機數(shù)驗證模型中組隨機數(shù)驗證模型中的關系表達式是否成立；的關系表達式是否成立；計算計算n 次模擬實驗中，使得關系表達次模擬實驗中，使得關系表達式成立的次數(shù)式成立的次數(shù)k ；當當 時，以時，以 作為此人能趕作為此人能趕上火車的概率上火車的概率p 的近似估計；的近似估計；nknwindows(7, 3)prb = replicate(100, x = sample(c(0, 5, 10), 1, prob

21、= c(0.7, 0.2, 0.1) y = sample(c(28, 30, 32, 34), 1, prob = c(0.3, 0.4, 0.2, 0.1) plot(0:40, rep(1, 41), type = n, xlab = time, ylab = ,axes = FALSE) axis(1, 0:40) r = rnorm(1, 30, 2) points(x, 1, pch = 15) i = 0 while (i = y) points(y, 1, pch = 19) Sys.sleep(0.1) points(y, 1, pch = 19) title(ifelse(

22、x + r y )mean(prb)1.1 矩母函數(shù)和生成函數(shù) 定義1.1 設隨機變量X的密度函數(shù)為p(x),稱 為X的矩母函數(shù)，記為性質：(1) (2) 若X和Y相互獨立，則 Xgt 0nnE Xg X YXYgtgt gtexpEtX定義1.2若為離散型隨機變量，稱為其概率生成函數(shù)，記為性質(1)(2)若X和Y獨立， XE s s 2 1 , 1 1E XE X X YXYsss1.條件分布其中 為X的邊際分布例子：令X和Y的聯(lián)合密度函數(shù)為試求 ,|Xp x yp y xpx Xpx,01, 01,0 xyxyp x yelse11|43P XY2.條件數(shù)學期望命題： ,|,yf x y

23、dyE YXxyp y x dyf x y dy |E E YXE Y 3.條件方差條件方差公式222|Var YXEYE YXXE YXE YX |Var YE Var YXVar E YX例1.3 從某大學任意挑選一個學院，然后從此學院中任意挑選n個學生，令X表示這些學生中來自武漢市的人數(shù)，令Q 代表該學院來自武漢市的人數(shù)所占的比例，因為學院之間的比例不相同，因此Q也是一個隨機變量。若QU(0,1), X|Q=qB(n,q)，求X的方差。隨機過程 是一個隨機變量集合，狀態(tài)空間S 是隨機變量 取值集合，集合T 為指標集。指標集可以是離散的也可以是連續(xù)的。例：天氣變化情況是一個隨機過程 晴，晴

24、，多云，雨，雨， ,tX tTtX命題：設 為Poisson過程中事件發(fā)生的間隔時間序列，則 為獨立同分布的隨機變量序列，且共同分布為具有參數(shù) 的指數(shù)分布1234,.XXXX1234,.XXXX 例1.5 顧客依Poisson過程到達商店，速率為 人/小時。已知商店上午9:00開門。試求：到9:30時僅到一位顧客，而到11:30時總計已到5位顧客的概率。.4 2.非齊次Posson過程一、Markov鏈的定義1Markov鏈,) 1(,)(|) 1(1ninXinXjnXP) 0 (,) 1 (,01iXiX)(|) 1(inXjnXP有限馬氏鏈狀態(tài)空間是有限集I=0,1,2,，k2一步轉移概

25、率馬氏鏈在時刻n處于狀態(tài) i 的條件下，到時刻n+1轉移到狀態(tài) j 的條件概率，即|1iXjXPnn稱為在時刻n的一步轉移概率，記作)(npij注:馬氏鏈由 和條件概率 決定00P Xi11|nnnnP Xi Xi注:由于概率是非負的，且過程從一狀態(tài)出發(fā)，經(jīng)過一步轉移后，必到達狀態(tài)空間中的某個狀態(tài)一步轉移概率滿足3一步轉移矩陣稱為在時刻n的一步轉移矩陣 隨機 矩陣即有有限馬氏鏈狀態(tài)空間I=0，1，2，k)()()()()()(10111001001npnpnpnpnpnpPnn)()()()()()()()()(1011110001001npnpnpnpnpnpnpnpnpPkkkkkk4齊次馬氏鏈即則稱此馬氏鏈為齊次馬氏鏈（即關于時間為齊次）如果馬氏鏈的一步轉移概率)(npij與 n 無關，ijnnpiXjXP|15初始分布設)(00iXPip，Ii，如果對一切Ii都有0)(0ip1)(0ipIi稱)(0ip為馬氏鏈的初始分布注馬氏鏈在初始時刻有可能處于I中任意狀態(tài)，初始分布就是馬氏鏈在初始時刻的概率分布。6絕對分布概率分布)(iXPipnn，Ii，0n稱為馬氏鏈的絕對分布或稱絕對概率13452. 隨例1 （機游動）12 3 4 55),Iii 設一醉漢在， ， ， ，作隨機游動：如果現(xiàn)在位于點（1則下一時刻各以1/3概率向左或向右