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

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

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

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

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

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

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

7、車的概率p 的近似估計(jì)；,?進(jìn)入演示,windows(7, 3)prb = replicate(100, { #括號(hào)內(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),?進(jìn)入模擬,三、R軟件的統(tǒng)計(jì)模擬功能,1、R軟件優(yōu)秀的隨機(jī)數(shù)模擬功能,生產(chǎn)某概率分布的隨機(jī)數(shù)是實(shí)現(xiàn)統(tǒng)計(jì)模擬的前提條件，而使用R命令可以生成以下常用分布的隨機(jī)數(shù)：,三、R軟件的統(tǒng)計(jì)模擬功能,2、優(yōu)良的編程環(huán)境和編程語言,R所擁有的好的兼容性、拓展性和強(qiáng)大的內(nèi)置函數(shù)有利于統(tǒng)計(jì)模擬的實(shí)現(xiàn)。

11、,3、高效率的向量運(yùn)算功能,使用R擁有的向量運(yùn)算功能可以大大減少程序運(yùn)行的時(shí)間，提高程序運(yùn)行的效率。,? 下面以求解Pi的程序?yàn)槔右哉f明,未采用R向量運(yùn)算功能的程序?yàn)椋?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)},,引入向量運(yùn)算功能改進(jìn)后的程序?yà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í)行兩個(gè)程序，看看有什么差異 程序1 ...... 程序2,三、R軟件的統(tǒng)計(jì)模擬功能,四、應(yīng)用R軟件模擬驗(yàn)證大數(shù)定律,1、驗(yàn)證的大數(shù)定律有：,（1）伯努利大數(shù)定理—— 設(shè) 是 次獨(dú)立重復(fù)試驗(yàn)中事件 發(fā)生的次數(shù)。 是事件 在每次試驗(yàn)中發(fā)生的概率，則對(duì)于任意正數(shù) >0,有,,,,,,,（2）辛欽定理： 設(shè)隨機(jī)變量 相互獨(dú)立，服從同一

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

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

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

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

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