第四章R編程實踐

1、第四章 R 編程實踐自此，我們已經(jīng)對R 的功能有個全面的了解。一些數(shù)據(jù)分析都是在R 對話窗口( R Console )中進行的，但對于復(fù)雜的數(shù)據(jù)分析顯然是不行方便的。下面將從統(tǒng)計語言和編程角度來說明，我們將會了解一些在 R編程實踐中采用的簡單的概念。一、 循環(huán)和向量化相比下拉菜單式的程序(R下開發(fā)的統(tǒng)計包R Commander, R的一個優(yōu)勢在于它可以把一系列連續(xù)的操作簡單的程序化。這一點上和所有其他計算機編程語言是一致的，但R 有一些特性使得非專業(yè)人士也可以很簡單的編寫程序。和其他編程語言一樣，R有一些和C語言類似的控制結(jié)構(gòu)。(一)控制結(jié)構(gòu)1. 條件語句：條件語句常用于避免除零或負(fù)數(shù)等數(shù)學(xué)問

2、題。它有兩種形式：1)if ( 條件 ) 表達式 1 else 表達式 22)ifelse( 條件 ,yes,no)例如：if (x=0) sqrt(x) else NAifelse (x=0,sqrt(x),NA)2. 循環(huán)：有三種形式：1)for ( 變量 in 向量 ) 表達式2)while( 條件 ) 表達式3)repeat 表達式例如：假定我們有一個向量x ，對于向量x 中值為 b 的元素，我們把0 賦給另外一個等長向量y 的對應(yīng)元素，否則賦1。我們首先創(chuàng)建一個與向量x 等長的向量y ：y - numeric(length(x)for (i in 1:length(x) if (xi

3、 = b) yi - 0 else yi - 1幾個指令可以放在一個大括弧里面：for (i in 1:length(x) yi - 0x - rnorm(10,-5,0.1)y - rnorm(10,5,2)X - cbind(x,y) # 矩陣X的列名保持x和y for (i in 1:5/2) zi minimum) count-0while (count10) print(count)count=count+1#如果設(shè)置不正確while 循環(huán)可以能導(dǎo)致無限循環(huán)，要小心#computeEstimate()x0-1tol-1e-5repeatx1-rnorm(1,1,2)if (abs(x

4、1-x0)tol)breakelsex0-x1#next,returnfor(i in 1:100)if(i=20)#skipnext#do somthing here3. 分支控制語句switch(statement,list)switch(2, banana , apple , other )num-10mode z z for (i in 1:length(z) zi yx= b yx!= b x y X apply(X,2,mean)lapply() 可以用于一個列表對象：它的語法類似apply 并且返回一個列表對象。 forms lapply(forms,lm)sapply() 是函

5、數(shù) lapply() 一個更為靈活的變種，它可以接受向量或者矩陣作為主要參數(shù)，并且返回形式更為友好的結(jié)果，常常是表格方式。還有sweep() 、 tapply() 函數(shù)等。例如：attach(mtcars)tapply(mpg,cy1,mean)col.mean=apply(mtcars,2,mean)sweep(mtcars,2,col.mean,-)表示從數(shù)據(jù)框mtcars 中按列計算col.mean ，并從 mtcars 中減去。2、 用 R 寫程序一般情況下，一個 R程序以ASCII格式保存，擴展名為.R。如果一個工作要重復(fù)好多次，用 R 程序是一個不錯的選擇。在我們的第一個例子中，我

6、們想對三個不同種屬的鳥繪制一樣的圖，而 且數(shù)據(jù)在三個不同的文件中。我們將一步一步的演示看R用不同的方式去完成這個簡單的問題。首先，我們憑直覺連續(xù)鍵入一系列命令，而且預(yù)先分割圖形設(shè)備。layout(matrix(1:3,3,1) # 分割圖形界面data - read.table(Swal.dat) # 讀入數(shù)據(jù)plot(data$V1,data$V2,type=l)title(swallow) # 增加標(biāo)題data - read.table(Wren.dat)plot(data$V1,data$V2,type=l)title(wren)data - read.table(Dunn.dat)pl

7、ot(data$V1,data$V2,type=l)title(dunnock)字符#用于在程序中添加注釋行：R會自動跳過注釋行。第一個程序的問題是在我們加入另外一個物種數(shù)據(jù)時，它過長。此外，一些命令多次執(zhí)行， 因 此它們可以放在一起，在執(zhí)行的時候僅僅修改一些參數(shù)。這里的策略是把參數(shù)放到一個字符型的向 量中去，然后用下標(biāo)去訪問這些不同的值。layout(matrix(1:3,3,1) # 分割圖形界面species - c(swallow,wren,dunnock)file - c(Swal.dat,Wren.dat,Dunn.dat)for(i in 1:length(species) da

8、ta - read.table(filei) #讀入數(shù)據(jù)plot(data$V1,data$V2,type=l) title(speciesi) # 增加標(biāo)題 注意代碼 read.table() 里面的參數(shù)filei 上面沒有雙引號， 因為這個參數(shù)是字符型。現(xiàn)在的代碼比較緊湊。 它比較容易加入新的物種， 因為設(shè)置物種名字和數(shù)據(jù)文件的向量都程序 的前端。如果擴展名為.dat的數(shù)據(jù)文件在 R的工作目錄下面，程序可以正常運行，否則用戶要設(shè)置工作目錄，或者設(shè)置絕對路徑( 例如： file source(Mybirds.R)和所有以文件作為輸入對象的函數(shù)一樣，如果該文件不在當(dāng)前工作目錄下面，用戶需要提供

9、該文件的絕對路徑。3、 編寫你自己的函數(shù)(一)編寫自己的函數(shù)1. 函數(shù)的基本操作大多數(shù) R 的工作是通過函數(shù)來實現(xiàn)，而且這些函數(shù)的輸入?yún)?shù)都放在一個括弧里面。用戶可以編寫他們自己的函數(shù)，并且這些函數(shù)和R里面的其他函數(shù)有一樣的特性。stdev=function(x)# 計算標(biāo)準(zhǔn)差 sqrt(var(x)stdevz=1:3stdev(z)函數(shù)的返回值可以是任何R對象，甚至可以是另一個函數(shù)?？梢酝ㄟ^顯式地調(diào)用return(),把一個值返回給主調(diào)函數(shù)。如果不用這條語句，默認(rèn)會將最后執(zhí)行的語句的值作為返回值。 例如：oddcount-function(x)k-0#k 的初始值 for(n in x)i

10、f(n%2=1)k=k+1#%是模 return(k)這個函數(shù)返回參數(shù)x 中奇數(shù)的個數(shù)。如果函數(shù)的輸出包括多個變量，可以把他們集中在一個列表中，以列表為參數(shù)調(diào)用函數(shù)，讓這個函數(shù)返回列表。創(chuàng)建函數(shù)之后可以通過函數(shù)formals() 和 body() 來查看其參數(shù)。formals(oddcount)body(oddcount)可以使用 page() 一個函數(shù)的代碼，確定它的功能和用法。page(abline)函數(shù)內(nèi)變量的賦值只在函數(shù)內(nèi)有效。編寫自己的函數(shù)可以讓你有效的，靈活的與合理的使用R。我們再次使用前面讀數(shù)據(jù)并且畫圖的例子。如果我們想在其他情況下進行這樣的操作，寫一個函數(shù)是一個不錯的想法：my

11、fun - function(S,F) data - read.table(F)plot(data$V1,data$V2,type=l)title(S) 執(zhí)行時，這個函數(shù)必須載入內(nèi)存。當(dāng)然，這有多種方式實現(xiàn)。和所有其他命令一樣，函數(shù)的各行可以直接通過鍵盤鍵入， 或者從一個文本編輯器里面拷貝粘貼。 如果函數(shù)保存在一個文本文件中，可以命令 source() 載入。如果用戶期望一些函數(shù)在 R 啟動時就被載入，可以把它們保存在工作目錄下面的文件.RData中。另外一種方式是，配置文件.Rpro le或Rpro - le( 詳見?Startup fordetails) 。最后，還可以創(chuàng)建一個包，但是這里

12、不想多討論。一旦函數(shù)載入后，我們就可以鍵入一條命令以讀入數(shù)據(jù)和畫出我們想要的圖，如myfun(swallow,Swal.dat) 。因此，現(xiàn)在我們的程序有第三個實現(xiàn)的版本了：layout(matrix(1:3,3,1)myfun(swallow,Swal.dat)myfun(wren,Wrenn.dat)myfun(dunnock,Dunn.dat)我們還可以用 sapply() 實現(xiàn)程序的第四個版本：layout(matrix(1:3,3,1)species - c(swallow,wren,dunnock)file foo x foo()1 1x 不是為了在函數(shù)foo() 里面創(chuàng)建對象，

13、因此 R 將會在封裝環(huán)境中搜索是否有個名為 x 的對象，和打印它的值( 否則一條錯誤信息將會顯示，執(zhí)行中斷) 。如果 x 是我們定義的函數(shù)中一個對象的名字，全局環(huán)境中變量x 值將會被采用。 x foo2 - function() x foo2()1 2 x1 1此時， print() 使用在它所在的環(huán)境中定義的x ，即環(huán)境foo2 中的 x 。前面提及的 封裝 是關(guān)鍵所在。在前面兩個演示函數(shù)中，有兩個環(huán)境：全局環(huán)境和函數(shù)foo 或foo2 的局部環(huán)境。 如果有三個或者更多的嵌套環(huán)境，對象搜索將逐層搜索直到全局環(huán)境。rm(list=ls()ls()ls.str() 可以獲得更多信息。有兩種方式指

14、定一個函數(shù)的參數(shù)： 通過它們的定義時的位置或者名字( 又稱為標(biāo)簽參數(shù)) 。 例如，假定一個函數(shù)有三個參數(shù)：foo - function(arg1,arg2,arg3) foo()在執(zhí)行時可以不用名字argl,，如果相應(yīng)的參數(shù)對象放在相應(yīng)的位置上，如：foo(x,y,z) 。但是，如果使用了參數(shù)的名字，位置信息將會失效，如 foo(arg3 = z,arg2 = y,arg1 = x) 。 R 函數(shù)的另外一個特性是函數(shù)可能采用定義時的默認(rèn)設(shè)置。例如：foo - function(arg1,arg2 = 5,arg3 = FALSE) 命令 foo(x) ， foo(x,5,FALSE) 和 fo

15、o(x,arg3 = FALSE) 將會得到一樣的結(jié)果。使用一個函數(shù)的默認(rèn)設(shè)置非常有用，特別在使用標(biāo)簽參數(shù)的時候，如 foo(x,arg3 = TRUE) 僅僅改變一個默認(rèn)設(shè) 置。環(huán)境層次中，某一層次的代碼對它上級層次中所有變量至少有讀的權(quán)限。但是，通過標(biāo)準(zhǔn)的運 算符 -直接對上級層次的變量進行寫操作是不可行的。#觀察 w 的值是否有變化w-12f-function(y)d-8w-w+1y-y-2cat(w=,w,n)h-function() return(d*(w+y) return(h() f(2) cat(w=,w,n) 如果你希望對一個全局變量，或者更一般情況下，對當(dāng)前層次的上級層次中

16、任意變量進行寫操 作，可以在當(dāng)前層次使用超賦值運算符- ，或者函數(shù)assign() 。(二)函數(shù)的調(diào)用函數(shù)也是對象，也可以給它們賦值，把它們用作其他函數(shù)的參數(shù)，如下所示：#函數(shù)的調(diào)用f1-function(a,b) return(a+b)f2-function(a,b) return(a-b)f-f1f(3,2)f-f2f(3,2)g-function(h,a,b) h(a,b)g(f1,3,2)g(f2,3,2)#在同一幅畫上繪制若干個函數(shù)的圖形g1-function(x) return(5*sin(x)+5)g2-function(x) return(sqrt(xA2+1)g3-funct

17、ion(x) return(2*x-1)plot(c(0,6.5),c(-5,15)# 準(zhǔn)備繪圖，確定 x和y的區(qū)域for(f in c(g1,g2,g3)plot(f,0,6.5,add=T)#在原有的圖形上繪圖(三)函數(shù)編輯R提供了一個edit函數(shù)，讓使用者可以在一個特定的窗口中更改自己的函數(shù)；fix使用預(yù)設(shè)的編輯器。fix(stdev)edit(file=for.r)我們來看另外一個例子。盡管這個例子不是純的統(tǒng)計學(xué)例子，但是它很好地展示了R語言的靈活性。假定我們想研究一個非線性模型的行為：這個模型( Ricker模型)的定義如下：Nt 1 Nt exp r 1 N K這個模型廣泛地用于種

18、群動態(tài)變化( population dynamics )的研究里面，特別是魚類的種群變 化。我們想用一個函數(shù)去模擬這個模型關(guān)于增長率r和初始群體大小 N的變化情況(承載能力K常常設(shè)定為1且以這個值作為默認(rèn)值)的影響；結(jié)果將以種群大小相對時間的圖表示。我們還將設(shè)定一個可選項允許用戶只顯示最后若干步中種群大小(默認(rèn)所有結(jié)果都會被繪制出來)。下面的函數(shù)就是做Ricker模型的數(shù)值模擬的。ricker - function(nzero,r,K=1,time=100,from=0,to=time) N - numeric(time+1)N1 - nzerofor (i in 1:time) Ni+1 -

19、 Ni*exp(r*(1 - Ni/K)Time layout(matrix(1:3,3,1) ricker(0.01,0.1);title(r = 0.1) ricker(0.05,0.2);title(r = 0.2) ricker(0.1,0.3);title(r = 0.3)四、防錯、除錯與調(diào)試1 .防錯與除錯R具備在函數(shù)中加入防錯功能的應(yīng)用函數(shù)：stop()和warning。stop()會傳回錯誤信息并終止執(zhí)行函數(shù)；warning。則在傳回錯誤訊息的同時繼續(xù)執(zhí)行函數(shù)。stdev=function(x) #計算標(biāo)準(zhǔn)差if (!is.vector(x)stop(向量才能計算標(biāo)準(zhǔn)差)7el

20、se sqrt(var(x)z=3z=as.matrix(z)stdev(z)# 利用輸出cat 語句進行確認(rèn)g=function(x,y)xA2-y+1z-0x-yA2+3*g(z,2)cat(x=,x,n)#利用print 語句進行調(diào)試x-yA2+3*g(z,2)cat(x=,x,n)w-28q0)u-1print(the if was done)elsev-10print(the else was done)2. 調(diào)試R 的核心調(diào)試工具由瀏覽器(browser )構(gòu)成，它可以讓你逐行運行代碼， 并在運行過程中檢查?？梢杂糜袃?nèi)置除錯函數(shù)debug() 或 browser()函數(shù)來打開這一瀏

21、覽器。R 的調(diào)試工具是針對單個函數(shù)的 ，如果你認(rèn)為函數(shù) f() 中有漏洞，就可以調(diào)用 debug(f) 來設(shè)置 函數(shù)的調(diào)試狀態(tài)。調(diào)用undebug(f) 將取消函數(shù)的調(diào)試狀態(tài)。debug(stdev)stdev(z)如果你只想對f() 進行一次調(diào)試會話，可以使用 debugonce(f) 函數(shù)，會在f() 第一次執(zhí)行時對其設(shè)置調(diào)試狀態(tài)。你果你認(rèn)為漏洞離函數(shù)的開始部分較遠，可以在函數(shù)的某一行加入 browser() 語句，那么瀏覽器只在程序執(zhí)行到這一行時被打開。#設(shè)置斷點f=function(x)xA3-2*x-1root=function(x,y)if(f(x)*f(y)=0)browser(

22、)x1-x;x20)stop( 區(qū)間端點的 f() 值之正負(fù)號不相反)elseif(f(x)0) x1-y;x2-xrepeatx00)x2-x0elsex1-x0if (abs(x1-x2)為變更為Browserd ( d 表示函數(shù)調(diào)用鏈的深度) 。在命令提示符之后你可以輸入以下命令：n (代表 next ) ：執(zhí)行下一行然后暫停，直接鍵入 Enter 也一樣。c (代表continue ):與n類似，但是如果當(dāng)前在某個循環(huán)中，將會執(zhí)行剩余的循環(huán)退出循環(huán)才會暫停。如果你不在某個循環(huán)中，那么下一次暫停前函數(shù)的剩余部分將被執(zhí)行。任意的R命令：在瀏覽器中，你依然處在R的交互模式，因此可以任意查看變量的取值。where ：會輸出一份棧跟蹤路徑，即顯示到達當(dāng)前