生物信息學資料-r相關statistical-simulation-based-on-r

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

2、連續(xù)隨機模擬和離散隨機模擬。而按變量是否隨時間變化又可分為動態(tài)隨機模擬和靜態(tài)隨機模擬。常用的統(tǒng)計模擬方法主要有以下幾種：1. 蒙特卡羅法2. 系統(tǒng)模擬方法3. 其它方法：包括Bootstrap(自助法)、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：

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

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

5、是否使用隨機產(chǎn)生的n 組隨機數(shù)驗證模型中的關系表達式是否成立；計算n 次模擬實驗中，使得關系表達式成立的次數(shù)k ；否試驗次數(shù)是否達到n次當n 時，以k作為此人能趕上n火車的概率p 的近似估計； 進入演示是計算估計結果k/n成立次數(shù)k=k+1成立次數(shù)不變驗證模型產(chǎn)生隨機數(shù)編寫R程序windows(7, 3)prb = replicate(100, #括號內(nèi)程序重復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)計模擬功能1、R軟件優(yōu)秀的隨機數(shù)模擬功能生產(chǎn)某概率分布的隨機數(shù)是實現(xiàn)統(tǒng)計模擬的前提條件， 而使用R命令可以生成以下常用分布的隨機數(shù)：分

7、布產(chǎn)生隨機數(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)計模擬功能2、優(yōu)良的編程環(huán)境和編程語言R所擁有的好的兼容性、拓展性和強大的內(nèi)置函數(shù)有利于統(tǒng)計模擬的實現(xiàn)。3、高效率的向量運算功能使用R擁有的向量運算功能可以大大減少程序運行的時間

8、， 提高程序運行的效率。 下面以求解Pi的程序為例加以說明三、R軟件的統(tǒng)計模擬功能引入向量運算功能改進后的程序為：未采用R向量運算功能的程序為：mc1-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í)行兩個程序，看看有什么差異 程序1 程序2四、應用R軟件

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

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

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