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

基于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)系如下圖。實際問題統(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.其它方法：包括Bootstrap(自助法)、MCMC（馬氏鏈蒙特卡羅法）等。一、統(tǒng)計模擬的基本概念（三）統(tǒng)計模擬的一般步驟二、趕火車問題火車離站時刻13:0013:0513:10概率0.70.20.1一列列車從A站開往B站，某人每天趕往B站上車。他已經(jīng)了解到火車從A站到B站的運行時間是服從均值為30min，標(biāo)準(zhǔn)差為2min的正態(tài)隨機變量?；疖嚧蠹s下午13：00離開A站，此人大約13:30到達B站。火車離開A站的時刻及概率如表1所示，此人到達B站的時刻及概率如表2所示。問此人能趕上火車的概率有多大？表1：火車離開A站的時刻及概率

表2：某人到達B站的時刻及概率人到站時刻13:2813:3013:3213:34概率0.30.40.20.1二、趕火車問題——問題的分析——

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

進行計算機統(tǒng)計模擬的基礎(chǔ)是抽象現(xiàn)實系統(tǒng)的數(shù)學(xué)模型

為了便于建模，對模型中使用的變量作出如下假定：

此人能及時趕上火車的充分必要條件為：，所以此人能趕上火車的概率模型為：。二、趕火車問題

為了分析簡化，假定13時為時刻t=0，則變量、的分布律為：05100.70.20.1283032340.30.40.20.1二、趕火車問題

R軟件求解的總算法：關(guān)系式成立產(chǎn)生隨機數(shù)驗證模型成立次數(shù)k=k+1否是計算估計結(jié)果k/n成立次數(shù)不變試驗次數(shù)是否達到n次是否編寫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；⑤當(dāng)時，以作為此人能趕上火車的概率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.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<=r){

i=i+1

segments(x,1,x+i,1)

if(x+i>=y)

points(y,1,pch=19)

Sys.sleep(0.1)

}

points(y,1,pch=19)

title(ifelse(x+r<=y,