生物信息學(xué)資料-r相關(guān)statistical-simulation-based-on-r

1、基于R軟件的統(tǒng)計(jì)模擬奚潭（南京財(cái)經(jīng)大學(xué)統(tǒng)計(jì)系2006級(jí)）主要內(nèi)容1. 統(tǒng)計(jì)模擬的基本概念2. 趕火車問(wèn)題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í)際問(wèn)題之間的橋梁。它與統(tǒng)計(jì)建模的關(guān)系如下圖。實(shí)際問(wèn)題 統(tǒng)計(jì)、邏輯 模型 計(jì)算機(jī)模擬（程序、算法） 實(shí)際解 統(tǒng)計(jì)、計(jì)算機(jī)解 一、統(tǒng)計(jì)模擬的基本概念（二）統(tǒng)計(jì)模擬方法一般地，統(tǒng)計(jì)模擬分類如下：若按狀態(tài)變量的變化性質(zhì)分為

2、連續(xù)隨機(jī)模擬和離散隨機(jī)模擬。而按變量是否隨時(shí)間變化又可分為動(dòng)態(tài)隨機(jī)模擬和靜態(tài)隨機(jī)模擬。常用的統(tǒng)計(jì)模擬方法主要有以下幾種：1. 蒙特卡羅法2. 系統(tǒng)模擬方法3. 其它方法：包括Bootstrap(自助法)、MCMC（馬氏鏈蒙特卡羅法）等。一、統(tǒng)計(jì)模擬的基本概念（三）統(tǒng)計(jì)模擬的一般步驟二、趕火車問(wèn)題一列列車從A站開(kāi)往B站，某人每天趕往B站上車。他已經(jīng)了解到火車從A站到B站的運(yùn)行時(shí)間是服從均值為30min，標(biāo)準(zhǔn)差為2min的正態(tài)隨量?；疖嚧蠹s下午13：00離開(kāi)A站，此人大約13:30到達(dá)B站?；疖囯x開(kāi)A站的時(shí)刻及概率如表1所示，此人到達(dá)B站的時(shí)刻及概率如表2所示。問(wèn)此人能趕上火車的概率有多大？表1：

3、火車離開(kāi)A站的時(shí)刻及概率表2：某人到達(dá)B站的時(shí)刻及概人到站時(shí)刻率13:2813:3013:3213:34概率0.30.40.20.1火車離站時(shí)刻13:0013:0513:10概率0.70.20.1二、趕火車問(wèn)題問(wèn)題的分析這個(gè)問(wèn)題用概率論的方法求解十分困難， 它涉及此人到達(dá)時(shí)刻、火車離開(kāi)站的時(shí)刻、火車運(yùn)行時(shí)間幾個(gè)隨量，而且火車運(yùn)行時(shí)間是服從正態(tài)分布的隨量，沒(méi)有有效的解析方法來(lái)進(jìn)行概率計(jì)算。在這種情況下可以用計(jì)算機(jī)模擬的方法來(lái)解決。二、趕火車問(wèn)題進(jìn)行計(jì)算機(jī)統(tǒng)計(jì)模擬的基礎(chǔ)是抽象現(xiàn)實(shí)系統(tǒng)的數(shù)學(xué)模型為了便于建模，對(duì)模型中使用的變量作出如下假定：T1 ：火車從A站出發(fā)的時(shí)刻；：火車從A站到B站的運(yùn)行時(shí)間；

4、T2：某人到達(dá)B站的時(shí)刻；T3m ：隨s：隨量 T 服2從正態(tài)分布的均值；量 T 服2從正態(tài)分布的標(biāo)準(zhǔn)差；二、趕火車問(wèn)題為了分析簡(jiǎn)化，假定13時(shí)為時(shí)刻t=0，則變量T1、 T3的分布律為：此人能及時(shí)趕上火車的充分必要條件為：T1，+ T2 T3所以此人能趕上火車的概率模型為：pT1。+ T2 T3T3 / min28303234P(t)0.30.40.20.1T1 / min0510P(t)0.70.20.1二、趕火車問(wèn)題R軟件求解的總算法：借助區(qū)間(0,1)分布產(chǎn)生的隨機(jī)數(shù)，對(duì)變量、T概T率分布進(jìn)行統(tǒng)計(jì)模擬；13根據(jù)變量T、1T、2概T3率分布及模擬程序、關(guān)系式成立命令產(chǎn)生n 個(gè)隨機(jī)分布數(shù)；

5、是否使用隨機(jī)產(chǎn)生的n 組隨機(jī)數(shù)驗(yàn)證模型中的關(guān)系表達(dá)式是否成立；計(jì)算n 次模擬實(shí)驗(yàn)中，使得關(guān)系表達(dá)式成立的次數(shù)k ；否試驗(yàn)次數(shù)是否達(dá)到n次當(dāng)n 時(shí)，以k作為此人能趕上n火車的概率p 的近似估計(jì)； 進(jìn)入演示是計(jì)算估計(jì)結(jié)果k/n成立次數(shù)k=k+1成立次數(shù)不變驗(yàn)證模型產(chǎn)生隨機(jī)數(shù)編寫(xiě)R程序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.4, 0.2, 0.1)plo

6、t(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 = 0while (i = y)points(y, 1, pch = 19) Sys.sleep(0.1)points(y, 1, pch = 19)title(ifelse(x + r y)mean(prb)三、R軟件的統(tǒng)計(jì)模擬功能1、R軟件優(yōu)秀的隨機(jī)數(shù)模擬功能生產(chǎn)某概率分布的隨機(jī)數(shù)是實(shí)現(xiàn)統(tǒng)計(jì)模擬的前提條件， 而使用R命令可以生成以下常用分布的隨機(jī)數(shù)：分

7、布產(chǎn)生隨機(jī)數(shù)序列命令參數(shù)設(shè)置binomialrbinom()n, size, probchi-squaredrchisq()n, df, ncpexponentialexp()n, rateFF()n, df1, df2, ncpnormalnorm()n,mean,sdPoissonpois()n, lambdaStudents tt()n, df, ncpunifomunif()n, min, max三、R軟件的統(tǒng)計(jì)模擬功能2、優(yōu)良的編程環(huán)境和編程語(yǔ)言R所擁有的好的兼容性、拓展性和強(qiáng)大的內(nèi)置函數(shù)有利于統(tǒng)計(jì)模擬的實(shí)現(xiàn)。3、高效率的向量運(yùn)算功能使用R擁有的向量運(yùn)算功能可以大大減少程序運(yùn)行的時(shí)間

8、， 提高程序運(yùn)行的效率。 下面以求解Pi的程序?yàn)槔右哉f(shuō)明三、R軟件的統(tǒng)計(jì)模擬功能引入向量運(yùn)算功能改進(jìn)后的程序?yàn)椋何床捎肦向量運(yùn)算功能的程序?yàn)椋簃c1-function(n) set.seed(1234579) k-0;x-runif(n);y-runif(n);for(i in 1:n) if(xi2+yi21) 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(xx2+y2 下面用R軟件分別執(zhí)行兩個(gè)程序，看看有什么差異 程序1 程序2四、應(yīng)用R軟件

9、模擬驗(yàn)證大數(shù)定律1、驗(yàn)證的大數(shù)定律有：（1）伯努利大數(shù)定理發(fā)生A的次數(shù)。在每A設(shè)是nA次n獨(dú)立重復(fù)試驗(yàn)中是次試驗(yàn)中發(fā)生的概率，則對(duì)于任意正數(shù)0,有elim P| nA - p | e = 1nn（2）辛設(shè)隨期望理：量X1 , X 2 , 相,互X n獨(dú),立，服從同一分布，且具有數(shù)學(xué)e，E( X，k )則= 對(duì)m(于k =任1,意2, 正)數(shù)，有l(wèi)im P 1X- m e = 1nk =1k nn四、應(yīng)用R軟件模擬驗(yàn)證大數(shù)定律2、在R軟件實(shí)現(xiàn)的算法思想：由大數(shù)定律可知，當(dāng) n ，樣本的均值趨向與理論分布的期望，因此利用樣本容量n逐漸增大這一趨勢(shì)來(lái)模擬 n 這一趨勢(shì)，在這種趨勢(shì)下，樣本的均值與理論

10、分布期望的誤差e 應(yīng)該呈現(xiàn)出越來(lái)越小的趨勢(shì)，同時(shí)，根據(jù)上述思想，分別對(duì)五種常用分布下的大數(shù)定律進(jìn)行驗(yàn)證。四、應(yīng)用R軟件模擬驗(yàn)證大數(shù)定律大數(shù)定律模擬算法設(shè)置循環(huán)的跳躍步長(zhǎng)ste、ps的n第一次抽樣的n2樣本容量初始值 和上n限1 值 ；利用函數(shù) seq( from = n1,to = n2 ,by = steps)產(chǎn)生由各模擬樣本空間大小組成的m 維序列；選擇隨機(jī)數(shù)Xi的分布類型，本文中的相關(guān)程序僅選擇了常用的隨機(jī)分布：正態(tài)分布、指數(shù)分布、均勻分布、泊松分布、二項(xiàng)分布、兩點(diǎn)分布；利用R軟件產(chǎn)生n個(gè)服從同一分布的隨機(jī)數(shù)X i (i = 1,2,.；., n)試驗(yàn)次數(shù)是否達(dá)到m次nn1n X（k計(jì)算

11、或 A ）的值；nk =1否若循環(huán)次數(shù) i之0(k和= 1, 2的,標(biāo))準(zhǔn)化變量：kknnnn Xk- E( Xk ) Xk- nmYn = i=1k =1 i=1nsnD( X )ki=1的分布函數(shù)對(duì)于任意滿足： n Xk- nm1x2lim F (x) = lim P i=1 x =e-t/2dt = F(x)nns2p-nn五、應(yīng)用R軟件模擬驗(yàn)證中心極限定理（2）De Moivre-Laplace(棣莫弗-拉普拉斯)中心極限定理設(shè)相互獨(dú)立的隨對(duì)于任意實(shí)數(shù) x ，有量hn (n =服1,從2, 參)數(shù)為 p 的兩點(diǎn)分布，則n Xi- np1x x = 2/2dt = f(x)e-tlim P i=12pnp(1- p)-n五、應(yīng)用R軟件模擬驗(yàn)證中心極限定理中心極限定理模擬算法選擇隨量X的i 分布類型，主要分布類型有正態(tài)分布、指數(shù)分布、均勻分布、泊松分布、二項(xiàng)分布和兩點(diǎn)分布；設(shè)置模擬試驗(yàn)總次數(shù)m 及每次模擬試驗(yàn)中隨量的個(gè)數(shù)n 的值;利用R軟件模擬產(chǎn)生n個(gè)服從同一分布的隨機(jī)數(shù)xi , (i = 1, 2,； , n)使用產(chǎn)生的n 個(gè)隨機(jī)數(shù)計(jì)算標(biāo)準(zhǔn)化隨量n值x - nmky j =