基于r軟件的統(tǒng)計模擬

1、基于R軟件的統(tǒng)計模擬,奚 潭（南京財經(jīng)大學(xué)統(tǒng)計系2006級）,主要內(nèi)容,1.統(tǒng)計模擬的基本概念,2.趕火車問題,3.R軟件的統(tǒng)計模擬功能,4.應(yīng)用R軟件模擬驗證大數(shù)定律,5.應(yīng)用R軟件模擬驗證中心極限定理,一、統(tǒng)計模擬的基本概念,（一）統(tǒng)計模擬的定義,統(tǒng)計模擬即是計算機統(tǒng)計模擬，它實質(zhì)上是計算機建模，而這里的計算機模型就是計算機方法、統(tǒng)計模型(如程序、流程圖、算法等)，它是架于計算機理論和實際問題之間的橋梁。它與統(tǒng)計建模的關(guān)系如下圖

2、。,實際問題,統(tǒng)計、邏輯 模型,計算機模擬（程序、算法）,統(tǒng)計、計算機解,實際解,,,,,一、統(tǒng)計模擬的基本概念,（二）統(tǒng)計模擬方法,一般地，統(tǒng)計模擬分類如下： 若按狀態(tài)變量的變化性質(zhì)分為連續(xù)隨機模擬和離散隨機模擬。 而按變量是否隨時間變化又可分為動態(tài)隨機模擬和靜態(tài)隨機模擬。 常用的統(tǒng)計模擬方法主要有以下幾種： 1.蒙特卡羅法 2.系統(tǒng)模擬方法 3.其它方法：包括Bootstr

3、ap(自助法)、MCMC（馬氏鏈蒙特卡羅法）等。,一、統(tǒng)計模擬的基本概念,（三）統(tǒng)計模擬的一般步驟,二、趕火車問題,一列列車從A站開往B站，某人每天趕往B站上車。他已經(jīng)了解到火車從A站到B站的運行時間是服從均值為30min，標準差為2min的正態(tài)隨機變量?；疖嚧蠹s下午13：00離開A站，此人大約13:30到達B站。火車離開A站的時刻及概率如表1所示，此人到達B站的時刻及概率如表2所示。問此人能趕上火車的概率有多大？,表1：火車離開A站的

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

5、服從正態(tài)分布的標準差；,二、趕火車問題,?進行計算機統(tǒng)計模擬的基礎(chǔ)是抽象現(xiàn)實系統(tǒng)的數(shù)學(xué)模型,?為了便于建模，對模型中使用的變量作出如下假定：,,,,,?此人能及時趕上火車的充分必要條件為： ，所以此人能趕上火車的概率模型為： 。,二、趕火車問題,,,?為了分析簡化，假定13時為時刻t=0，則變量 、 的分布律為：,二、趕火車問題,?R軟件求解的總算法：,,產(chǎn)生隨機數(shù),驗證模型,成立次

6、數(shù)k=k+1,,,,,,否,是,計算估計結(jié)果k/n,,成立次數(shù)不變,,,,,,,,,是,否,編寫R程序,,①借助區(qū)間(0,1)分布產(chǎn)生的隨機數(shù)，對變量 、 概率分布進行統(tǒng)計模擬；,②根據(jù)變量 、 、 概率分布及模擬程序、命令產(chǎn)生n 個隨機分布數(shù)；,③使用隨機產(chǎn)生的n 組隨機數(shù)驗證模型中的關(guān)系表達式是否成立；,④計算n 次模擬實驗中，使得關(guān)系表達式成立的次數(shù)k ；,⑤當 時，以 作為此人能趕上火

7、車的概率p 的近似估計；,?進入演示,windows(7, 3)prb = replicate(100, { #括號內(nèi)程序重復(fù)100次    x = sample(c(0, 5, 10), 1, prob = c(0.7, 0.2, 0.1))    y = sample(c(28, 30, 32, 34), 1, prob = c(0.3, 0.

8、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,

9、2)    points(x, 1, pch = 15)    i = 0    while (i = y)            points(y, 1, pch = 19)        Sys.sleep(0.1)    }  &

10、#160; points(y, 1, pch = 19)    title(ifelse(x + r y})mean(prb),?進入模擬,三、R軟件的統(tǒng)計模擬功能,1、R軟件優(yōu)秀的隨機數(shù)模擬功能,生產(chǎn)某概率分布的隨機數(shù)是實現(xiàn)統(tǒng)計模擬的前提條件，而使用R命令可以生成以下常用分布的隨機數(shù)：,三、R軟件的統(tǒng)計模擬功能,2、優(yōu)良的編程環(huán)境和編程語言,R所擁有的好的兼容性、拓展性和強大的內(nèi)置函數(shù)有利于統(tǒng)計模擬的實現(xiàn)。

11、,3、高效率的向量運算功能,使用R擁有的向量運算功能可以大大減少程序運行的時間，提高程序運行的效率。,? 下面以求解Pi的程序為例加以說明,未采用R向量運算功能的程序為：,mc1<-function(n){ set.seed(1234579) k<-0; x<-runif(n); y<-runif(n); for(i in 1:n){ if(x[i]^2+y[i]^2&l

12、t;1) k<-k+1; } data.frame(Pi=4*k/n)},,引入向量運算功能改進后的程序為：,mc1<-function(n){ set.seed(1234579) k<-0; x<-runif(n); y<-runif(n); k <- length(x[x^2+y^2 < 1]) data.frame(Pi=4*k

13、/n)},,--> 下面用R軟件分別執(zhí)行兩個程序，看看有什么差異 程序1 ...... 程序2,三、R軟件的統(tǒng)計模擬功能,四、應(yīng)用R軟件模擬驗證大數(shù)定律,1、驗證的大數(shù)定律有：,（1）伯努利大數(shù)定理—— 設(shè) 是 次獨立重復(fù)試驗中事件 發(fā)生的次數(shù)。 是事件 在每次試驗中發(fā)生的概率，則對于任意正數(shù) >0,有,,,,,,,（2）辛欽定理： 設(shè)隨機變量 相互獨立，服從同一

14、分布，且具有數(shù)學(xué)期望 ， ，則對于任意正數(shù) ，有,,,四、應(yīng)用R軟件模擬驗證大數(shù)定律,2、在R軟件實現(xiàn)的算法思想：,由大數(shù)定律可知，當 ，樣本的均值趨向與理論分布的期望，因此利用樣本容量 逐漸增大這一趨勢來模擬 這一趨勢，在這種趨勢下，樣本的均值與理論分布期望的誤差 應(yīng)該呈現(xiàn)出越來越小的趨勢，同時，根據(jù)上述思想，分別對五種常用分布下的大數(shù)定律進行驗

15、證。,,,,四、應(yīng)用R軟件模擬驗證大數(shù)定律,?大數(shù)定律模擬算法,,設(shè)置參數(shù)值,產(chǎn)生m維序列,,,,,繪圖,,,是,否,編寫R程序,,,選擇分布類型,,產(chǎn)生隨機數(shù),,計算樣本均值y,①設(shè)置循環(huán)的跳躍步長 、 的第一次抽樣的樣本容量初始值 和上限值 ；,②利用函數(shù) 產(chǎn)生由各模擬樣本空間大小組成的m 維序列；,③選擇隨機數(shù) 的

16、分布類型，本文中的相關(guān)程序僅選擇了常用的隨機分布：正態(tài)分布、指數(shù)分布、均勻分布、泊松分布、二項分布、兩點分布；,,④利用R軟件產(chǎn)生n個服從同一分布的隨機數(shù) ；,,⑤計算 （或 ）的值；,,,⑥若循環(huán)次數(shù) i<m ，則回轉(zhuǎn)④，否則轉(zhuǎn)⑦;,⑦以x軸代表樣本容量n ，y 軸代表每次抽樣所得的樣本均值，描繪出整個試驗的過程。,? 進入演示……,五、應(yīng)用R軟件模擬驗證中心極限定理,1、驗證的中

17、心極限定理有,（1）獨立同分布的中心極限定理： 設(shè)隨機變量 相互獨立，服從同一分布，且具有數(shù)學(xué)期望和方差： ，則隨機變量之和 的標準化變量：,,的分布函數(shù)對于任意滿足：,,（2）De Moivre-Laplace(棣莫弗-拉普拉斯)中心極限定理,,,設(shè)相互獨立的隨機變量 服從參數(shù)為 p 的兩點分布，則對于任意實數(shù) x ，有,

18、,,五、應(yīng)用R軟件模擬驗證中心極限定理,,選擇分布類型,確定參數(shù)m和n,,,,,統(tǒng)計檢驗和描述性分析,,,是,否,編寫R程序,,,產(chǎn)生隨機數(shù),,計算標準化隨機變量,,設(shè)置參數(shù)j和step,?中心極限定理模擬算法,①選擇隨機變量 的分布類型，主要分布類型有正態(tài)分布、指數(shù)分布、均勻分布、泊松分布、二項分布和兩點分布；,②設(shè)置模擬試驗總次數(shù)m 及每次模擬試驗中隨機變量的個數(shù)n 的值;,③利用R軟件模擬產(chǎn)生n個服從同一分布的隨機數(shù)