統(tǒng)計(jì)模擬介紹

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

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

3、的發(fā)展歷史方法的發(fā)展歷史1 1、BuffonBuffon投針實(shí)驗(yàn)：投針實(shí)驗(yàn)：17681768年，法國(guó)數(shù)學(xué)家年，法國(guó)數(shù)學(xué)家Comte de Buffon利用投針實(shí)驗(yàn)估計(jì)利用投針實(shí)驗(yàn)估計(jì) 的的值值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年，古稀之年的蒲豐在家中請(qǐng)來(lái)好些客人玩投針游戲（針長(zhǎng)是線距之半），他事先沒(méi)有給客人講與有關(guān)的事?？腿藗冸m然不知道主人的用意，但是都參加了游戲。他們共投針2212次，其中704次相交。蒲豐說(shuō)，2212/704=3.142，這就是值。這著實(shí)讓人們驚喜不已。2 2、19301930年，年，Enrico FermiEnrico Fermi利用利用Monte CarloMonte Carlo方法研究中子方法研究中子的擴(kuò)散，并設(shè)計(jì)了一個(gè)的擴(kuò)散，并設(shè)計(jì)了一個(gè)Monte CarloMonte Carlo機(jī)械裝置，機(jī)械裝置，F(xiàn)ermiac,F

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

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

10、當(dāng)所求問(wèn)題的解是某個(gè)事件的概率，或者是某個(gè)隨機(jī)變量的數(shù)學(xué)期望，或者是與概率、數(shù)學(xué)期望有關(guān)的量時(shí)，通過(guò)某種試驗(yàn)的方法，得出該事件發(fā)生的頻率，或者該隨機(jī)變量若干個(gè)具體觀察值的算術(shù)平均值，通過(guò)它得到問(wèn)題的解。這就是蒙特卡羅方法的基本思想。 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關(guān)系式成立產(chǎn)生隨機(jī)數(shù)驗(yàn)證模型成立次數(shù)k=k+1否是計(jì)算估計(jì)結(jié)果k/n成立次數(shù)不變?cè)囼?yàn)次數(shù)是否達(dá)到n次是否編寫(xiě)R程序Monte CarloMonte Carlo模擬的應(yīng)

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

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

13、法來(lái)求解物理問(wèn)方法利用隨機(jī)抽樣的方法來(lái)求解物理問(wèn)題題; ;數(shù)值解法數(shù)值解法: :從一個(gè)物理系統(tǒng)的數(shù)學(xué)模型出發(fā)從一個(gè)物理系統(tǒng)的數(shù)學(xué)模型出發(fā), ,通過(guò)求解一通過(guò)求解一系列的微分方程來(lái)的導(dǎo)出系統(tǒng)的未知狀態(tài)系列的微分方程來(lái)的導(dǎo)出系統(tǒng)的未知狀態(tài); ;Monte CarloMonte Carlo方法并非只能用來(lái)解決包含隨機(jī)的過(guò)程的問(wèn)題方法并非只能用來(lái)解決包含隨機(jī)的過(guò)程的問(wèn)題: :許多利用許多利用Monte CarloMonte Carlo方法進(jìn)行求解的問(wèn)題中并不包含方法進(jìn)行求解的問(wèn)題中并不包含隨機(jī)過(guò)程隨機(jī)過(guò)程 例如例如: :用用Monte CarloMonte Carlo方法計(jì)算定積分方法計(jì)算定積分. .

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

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

16、us, R, Gauss, minitab優(yōu)點(diǎn) 能夠比較逼真地描述具有隨機(jī)性質(zhì)的事物的特點(diǎn)及物理實(shí)驗(yàn)過(guò)程。 受幾何條件限制小。 收斂速度與問(wèn)題的維數(shù)無(wú)關(guān)。 具有同時(shí)計(jì)算多個(gè)方案與多個(gè)未知量的能力。 誤差容易確定。 程序結(jié)構(gòu)簡(jiǎn)單，易于實(shí)現(xiàn)。 缺點(diǎn) 收斂速度慢。 誤差具有概率性。 在粒子輸運(yùn)問(wèn)題中，計(jì)算結(jié)果與系統(tǒng)大小有關(guān)?；疖?chē)離站時(shí)刻火車(chē)離站時(shí)刻13:0013:0513:10概率概率0.70.20.1 一列列車(chē)從一列列車(chē)從A A站開(kāi)往站開(kāi)往B B站，某人每天趕往站，某人每天趕往B B站上車(chē)。他站上車(chē)。他已經(jīng)了解到火車(chē)從已經(jīng)了解到火車(chē)從A A站到站到B B站的運(yùn)行時(shí)間是服從均值為站的運(yùn)行時(shí)間是服從均

17、值為30min30min，標(biāo)準(zhǔn)差為，標(biāo)準(zhǔn)差為2min2min的正態(tài)隨機(jī)變量?；疖?chē)大約下午的正態(tài)隨機(jī)變量?；疖?chē)大約下午1313：0000離開(kāi)離開(kāi)A A站，此人大約站，此人大約13:3013:30到達(dá)到達(dá)B B站?；疖?chē)離開(kāi)站。火車(chē)離開(kāi)A A站的時(shí)刻站的時(shí)刻及概率如表及概率如表1 1所示，此人到達(dá)所示，此人到達(dá)B B站的時(shí)刻及概率如表站的時(shí)刻及概率如表2 2所示。所示。問(wèn)此人能趕上火車(chē)的概率有多大？問(wèn)此人能趕上火車(chē)的概率有多大？表1：火車(chē)離開(kāi)A站的時(shí)刻及概率 表2：某人到達(dá)B站的時(shí)刻及概率 人到站時(shí)刻人到站時(shí)刻13:2813:3013:3213:34概率概率0.30.40.20.1這個(gè)問(wèn)題用概率論的

18、方法求解十分困難，它涉及此人到達(dá)時(shí)刻、火車(chē)離開(kāi)站的時(shí)刻、火車(chē)運(yùn)行時(shí)間幾個(gè)隨機(jī)變量，而且火車(chē)運(yùn)行時(shí)間是服從正態(tài)分布的隨機(jī)變量，沒(méi)有有效的解析方法來(lái)進(jìn)行概率計(jì)算。在這種情況下可以用計(jì)算機(jī)模擬的方法來(lái)解決。：火車(chē)從：火車(chē)從A A站出發(fā)的時(shí)刻；站出發(fā)的時(shí)刻；：火車(chē)從：火車(chē)從A A站到站到B B站的運(yùn)行時(shí)間；站的運(yùn)行時(shí)間；：某人到達(dá)：某人到達(dá)B B站的時(shí)刻；站的時(shí)刻；：隨機(jī)變量：隨機(jī)變量 服從正態(tài)分布的均值；服從正態(tài)分布的均值；：隨機(jī)變量：隨機(jī)變量 服從正態(tài)分布的標(biāo)準(zhǔn)差服從正態(tài)分布的標(biāo)準(zhǔn)差；1T2T3T2T2T此人能及時(shí)趕上火車(chē)的充分必要條件為：此人能及時(shí)趕上火車(chē)的充分必要條件為： ，所以此人能趕上火車(chē)

19、的概率模型為：所以此人能趕上火車(chē)的概率模型為： 。123TTT123p TTT為了分析簡(jiǎn)化，假定為了分析簡(jiǎn)化，假定1313時(shí)為時(shí)刻時(shí)為時(shí)刻t=0，則變量，則變量 、 的分布律為：的分布律為：1T3T05100.70.20.1283032340.30.40.20.11/minT( )P t3/minT( )P t關(guān)系式成立產(chǎn)生隨機(jī)數(shù)驗(yàn)證模型成立次數(shù)k=k+1否是計(jì)算估計(jì)結(jié)果k/n成立次數(shù)不變?cè)囼?yàn)次數(shù)是否達(dá)到n次是否編寫(xiě)R程序借助區(qū)間借助區(qū)間(0,1)(0,1)分布產(chǎn)生的隨機(jī)數(shù)，對(duì)分布產(chǎn)生的隨機(jī)數(shù)，對(duì)變量變量 、 概率分布進(jìn)行統(tǒng)計(jì)模擬；概率分布進(jìn)行統(tǒng)計(jì)模擬；1T3T根據(jù)變量根據(jù)變量 、 、 概率分

20、布及模擬概率分布及模擬程序、命令產(chǎn)生程序、命令產(chǎn)生n 個(gè)隨機(jī)分布數(shù)；個(gè)隨機(jī)分布數(shù)；1T2T3T使用隨機(jī)產(chǎn)生的使用隨機(jī)產(chǎn)生的n 組隨機(jī)數(shù)驗(yàn)證模型中組隨機(jī)數(shù)驗(yàn)證模型中的關(guān)系表達(dá)式是否成立；的關(guān)系表達(dá)式是否成立；計(jì)算計(jì)算n 次模擬實(shí)驗(yàn)中，使得關(guān)系表達(dá)次模擬實(shí)驗(yàn)中，使得關(guān)系表達(dá)式成立的次數(shù)式成立的次數(shù)k ；當(dāng)當(dāng) 時(shí)，以時(shí)，以 作為此人能趕作為此人能趕上火車(chē)的概率上火車(chē)的概率p 的近似估計(jì)；的近似估計(jì)；nknwindows(7, 3)prb = replicate(100, x = sample(c(0, 5, 10), 1, prob = c(0.7, 0.2, 0.1) y = sample(c(2

21、8, 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(x + r y )mean(prb)1.1 矩母函數(shù)和生成函數(shù) 定

22、義1.1 設(shè)隨機(jī)變量X的密度函數(shù)為p(x),稱(chēng) 為X的矩母函數(shù)，記為性質(zhì)：(1) (2) 若X和Y相互獨(dú)立，則 Xgt 0nnE Xg X YXYgtgt gtexpEtX定義1.2若為離散型隨機(jī)變量，稱(chēng)為其概率生成函數(shù)，記為性質(zhì)(1)(2)若X和Y獨(dú)立， XE s s 2 1 , 1 1E XE X X YXYsss 1.條件分布其中 為X的邊際分布例子：令X和Y的聯(lián)合密度函數(shù)為試求 ,|Xp x yp y xpx Xpx,01, 01,0 xyxyp x yelse11|43P XY2.條件數(shù)學(xué)期望命題： ,|,yfx y dyE YXxyp y x dyfx y dy |E E YXE

23、Y 3.條件方差條件方差公式222|Var YXEYE YXXE YXE YX |Var YE Var YXVar E YX例1.3 從某大學(xué)任意挑選一個(gè)學(xué)院，然后從此學(xué)院中任意挑選n個(gè)學(xué)生，令X表示這些學(xué)生中來(lái)自武漢市的人數(shù)，令Q 代表該學(xué)院來(lái)自武漢市的人數(shù)所占的比例，因?yàn)閷W(xué)院之間的比例不相同，因此Q也是一個(gè)隨機(jī)變量。若QU(0,1), X|Q=qB(n,q)，求X的方差。隨機(jī)過(guò)程 是一個(gè)隨機(jī)變量集合，狀態(tài)空間S 是隨機(jī)變量 取值集合，集合T 為指標(biāo)集。指標(biāo)集可以是離散的也可以是連續(xù)的。例：天氣變化情況是一個(gè)隨機(jī)過(guò)程 晴，晴，多云，雨，雨， ,tX tTtX命題：設(shè) 為Poisson過(guò)程中事

24、件發(fā)生的間隔時(shí)間序列，則 為獨(dú)立同分布的隨機(jī)變量序列，且共同分布為具有參數(shù) 的指數(shù)分布1234,.XXXX1234,.XXXX 例1.5 顧客依Poisson過(guò)程到達(dá)商店，速率為 人/小時(shí)。已知商店上午9:00開(kāi)門(mén)。試求：到9:30時(shí)僅到一位顧客，而到11:30時(shí)總計(jì)已到5位顧客的概率。.4 2.非齊次Posson過(guò)程一、Markov鏈的定義1Markov鏈,) 1(,)(|) 1(1ninXinXjnXP) 0 (,) 1 (,01iXiX)(|) 1(inXjnXP有限馬氏鏈狀態(tài)空間是有限集I=0,1,2,，k2一步轉(zhuǎn)移概率馬氏鏈在時(shí)刻n處于狀態(tài) i 的條件下，到時(shí)刻n+1轉(zhuǎn)移到狀態(tài) j

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