用R語言進行數(shù)據分析:編寫函數(shù)

前一段時間，作為數(shù)據分析師的我給大家分享過一些關于r語言的一些經驗，今天我(數(shù)據分析師)會在才分享一些，正如前面內容所暗示的一樣，R 語言允許用戶 創(chuàng)建自己的函數(shù)（function）對象。R 有一些 內部函數(shù)并且可以用在 其他的表達式中。通過這個過程，R 在程序的功能性， 便利性和優(yōu)美性上得到了擴展。學寫這些有用的函數(shù) 是一個人輕松地創(chuàng)造性地使用 R 的 最主要的方式。

需要強調的是，大多是函數(shù)都作為 R 系統(tǒng)的一部分提供，如mean(),var(),postscript()等等。這些函數(shù)都是用 R 寫的， 因此在本質上和用戶寫的沒有差別。

一個函數(shù)是通過下面的語句形式定義的。

>name<- function(arg_1,arg_2, ...)expression

其中expression是一個 R 表達式(常常是一個成組 表達式)，它利用參數(shù)arg_i計算最終的結果。 該表達式的值就是返回給函數(shù)的最終值。

可以在任何地方以name(expr_1,expr_2, ...)的形式調用函數(shù)。

Simple examples: 簡單的例子

這是一個簡單的例子，它用來計算雙樣本的 t-統(tǒng)計量，并且顯示“所有步驟”。這是一個人為的例子， 當然還有其他更簡單的辦法 得到一樣的結果。

函數(shù)定義如下：

     > twosam <- function(y1, y2) {
         n1  <- length(y1); n2  <- length(y2)
         yb1 <- mean(y1);   yb2 <- mean(y2)
         s1  <- var(y1);    s2  <- var(y2)
         s <- ((n1-1)*s1 + (n2-1)*s2)/(n1+n2-2)
         tst <- (yb1 - yb2)/sqrt(s*(1/n1 + 1/n2))
         tst
       }

通過這個函數(shù)，你可以下面的命令 實現(xiàn)雙樣本 t-檢驗。

     > tstat <- twosam(data$male, data$female); tstat

第二個例子是仿效 Matlab 里面的反斜杠命令。它是用來返回向量 y 正交投影到 X 列空間上面的系數(shù)。 (這常常被稱為回歸系數(shù)的 最小二乘法估計。) 這可以用 函數(shù)qr()來實現(xiàn)；但是直接使用這個函數(shù)有時有點 難處理。下面提供了一個簡單 而又安全的函數(shù)。

給定 n × 1 的向量 y 和一個 n × p 的矩陣 X，因此 X \ y 可以定義如下 (X’X)^{-}X’y, 其中 (X’X)^{-} 這就是 X’X的廣義逆矩陣（generalized inverse）。

     > bslash <- function(X, y) {
       X <- qr(X)
       qr.coef(X, y)
     }

當這個對象創(chuàng)建后，它可能用于這樣的命令中：

     > regcoeff <- bslash(Xmat, yvar)

經典的 R 函數(shù)lsfit()可以很好的實現(xiàn)這個功能和其他一些相關 的事情¹。它以一種有點違反直覺的方式 依次用函數(shù)qr()和qr.coef()去完成這部分計算。 如果一部分代碼常常使用， 我們可以把這部分代碼單獨列出來成為函數(shù) 使用。如果是這樣，我們可能期望矩陣的二元操作 能以一種更為便利的方式進行。

Defining new binary operators: 定義新的二元操作符

假定我們給予函數(shù)bslash()一個不同的名字，且以 下面的形式給出

     %anything%

那么它將以二元操作符的形式在表達式中使用， 而不是函數(shù)的形式。例如我們選擇!作為中間的字符。函數(shù)可以如下定義

     > "%!%" <- function(X, y) { ... }

(注意要使用引號)。該函數(shù)然后就可以用于X %!% y。(反斜杠符不是一個很好的選擇， 因為在某些情況下會引入一些特定的問題。)

矩陣的乘法操作符%*%和外積操作符%o%同樣是這種方式定義的 二元操作符。

Named arguments and defaults: 參數(shù)命名和默認值

和 Generating regular sequences 提示的一樣，如果調用的函數(shù)的參數(shù) 以“name=object”的方式給出， 它們可以用任何順序。但是，參數(shù)賦值序列可能以 未命名的，位置特異性的方式給出，同時也有可能 在這些位置特異性的參數(shù)后加上命名參數(shù)賦值。

因此，如果有下面方式定義的函數(shù)fun1

     > fun1 <- function(data, data.frame, graph, limit) { [function body omitted] }

那么函數(shù)將會被好幾種方式調用，如

     > ans <- fun1(d, df, TRUE, 20)
     > ans <- fun1(d, df, graph=TRUE, limit=20)
     > ans <- fun1(data=d, limit=20, graph=TRUE, data.frame=df)

上面所有的方式是等價的。

許多時候，參數(shù)會被設定一些默認值。 如果默認值適合你要做的事情，你可以省略這些參數(shù)。 例如，函數(shù)fun1用下面的方式 定義

     > fun1 <- function(data, data.frame, graph=TRUE, limit=20) { ... }

它可以被如下命令調用

     > ans <- fun1(d, df)

這和前面的三種情況等價。

     > ans <- fun1(d, df, limit=10)

這就改變了一個默認值。

特別說明一下，默認值可以是任何表達式，甚至是函數(shù)本身 所帶有的其他參數(shù)；它們沒有要求 是常數(shù)。我們的例子采用常數(shù)只是使問題簡單容易說明。

The ellipsis argument (…): 省略符號的參數(shù)(…)

還有一種常常出現(xiàn)的情況就是要求一個函數(shù)的參數(shù)設置 可以傳遞給另外一個函數(shù)。例如圖形函數(shù)如果調用了函數(shù)par()和其他如plot()類的函數(shù)，par()函數(shù)的圖形設置將會傳遞給圖形輸出的設備控制。 (See The par() function, 后面的章節(jié)會給出函數(shù)par()更為詳細的內容。)這個可以通過給函數(shù) 增加一個額外的參數(shù)來實現(xiàn)。這個參數(shù)字面上就是 …，它可以被傳遞。 一個概述性的例子可以如下所示。

     fun1 <- function(data, data.frame, graph=TRUE, limit=20, ...) { [省略一些語句] if (graph)
         par(pch="*", ...) [省略其他語句] }

Assignment within functions: 函數(shù)內部的賦值

注意任何在函數(shù)內部的普通賦值都是局部的 暫時的，當退出函數(shù)時都會丟失。因此 函數(shù)中的賦值語句X <- qr(X)不會影響 調用該函數(shù)的程序賦值情況。

如果要徹底理解 R 賦值管理的原則， 讀者需要熟悉解析 框架（Evalution frame）的概念。這屬于高級內容， 不在這里深入討論。

如果想在一個函數(shù)里面全局賦值或者永久賦值，那么可以采用 “強賦值”（superassignment）操作符<<-或者采用函數(shù)assign()。用help命令可以得到更細節(jié)的內容。 S-Plus 用戶需要注意<<-在 R 里面有著不同的語義（semantics）。

More advanced examples: 高級的例子

區(qū)組設計中的效率因子

下面是一個有點枯燥但較為完整的例子。它是用來 計算一個區(qū)組設計中的效率因子。

區(qū)組設計（block design）需要考慮兩個因子blocks(b個水平) 和varieties(v個水平)。如果R 和 K 分別是 v × v 和 b × b 重復（replications）及 區(qū)組大小（block size）矩陣而 N 則是 b × v 的關聯(lián)矩陣，那么 那么效率因子就是這個矩陣的特征值。 E = I_v – R^{-1/2}N’K^{-1}NR^{-1/2} = I_v – A’A, where A = K^{-1/2}NR^{-1/2}. 寫這個函數(shù)的一種方式就是

     > bdeff <- function(blocks, varieties) {
         blocks <- as.factor(blocks)             # minor safety move b <- length(levels(blocks))
         varieties <- as.factor(varieties)       # minor safety move v <- length(levels(varieties))
         K <- as.vector(table(blocks))           # 去掉 dim 屬性 R <- as.vector(table(varieties))       
 # 去掉 dim 屬性 N <- table(blocks, varieties)
         A <- 1/sqrt(K) * N * rep(1/sqrt(R), rep(b, v))
         sv <- svd(A)
         list(eff=1 - sv$d^2, blockcv=sv$u, varietycv=sv$v)
     }

這種情況下，奇異值分解 比求解特征值效率高。

函數(shù)的結果是一個列表。它不僅以第一個分量的形式給出了 效率因子，還給出了區(qū)組和規(guī)范對照信息， 因為有些時候這些會給出額外有用的 定量信息。cda數(shù)據分析師

去除打印數(shù)組中的名字

為了顯示一個大的數(shù)組或者矩陣，常常需要 需要以一個完整的塊的形式顯示，同時去掉數(shù)組名和編號。 去掉dimnames屬性是不能達到這個要求的，因為 R 環(huán)境會把空的字符串賦給dimnames屬性。 為了打印一個矩陣X

     > tempdimnames(temp) <- list(rep("", nrow(X)), rep("", ncol(X)))
     > temp; rm(temp)

這個可以非常便利地通過下面的函數(shù)no.dimnames()實現(xiàn)。它是利用一種“卷繞”（wrap around） 的方式實現(xiàn)的。這個例子還說明一些非常高效有用的用戶函數(shù) 也可以是非常簡潔的。

     no.dimnames <- function(a) {
       ## 為了更緊湊的打印，可以去除數(shù)組中的維度名字 d <- list()
       l <- 0
       for(i in dim(a)) {
         d[[l <- l + 1]] <- rep("", i)
       }
       dimnames(a) <- d
       a
     }

通過這個函數(shù)，數(shù)組可以用一種緊湊的方式 顯示

     > no.dimnames(X)

這對大的整數(shù)數(shù)非常的有用，因為這些數(shù)組 表現(xiàn)出來的式樣（pattern）可能比它們的值更為重要。

遞歸式的數(shù)值積分

函數(shù)可以是遞歸的，可以在函數(shù)內部調用自己。 但是需要注意的是，這些函數(shù)或者變量， 不會被更高層次的解析框架（evaluation frame）中的被調用函數(shù) 所繼承。如果它們在搜索路徑中，這種情況就會出現(xiàn)。

下面的例子顯示了一個最簡單的一維數(shù)值積分方法。 被積函數(shù)在所積范圍的兩端和中點的值會被計算。 如果單面板梯形法則（one-panel trapezium rule）的結果和 雙面板的非常相似，那么就以后者作為返回值。 否則同樣的過程會遞歸用于各個面板。 這是一個自適應的積分過程。它會集中各個 被積函數(shù)接近線性的區(qū)域函數(shù)計算值。 但是這種方法開銷有點過大¹。 相對其他積分算法，它的優(yōu)勢體現(xiàn)在被積函數(shù)既平滑又很難求值時。

這個例子同樣可以部分的作為一個 R 編程的難題給出。

     area <- function(f, a, b, eps = 1.0e-06, lim = 10) {
       fun1 <- function(f, a, b, fa, fb, a0, eps, lim, fun) {
         ## function `fun1' is only visible inside `area' d <- (a + b)/2
         h <- (b - a)/4
         fd <- f(d)
         a1 <- h * (fa + fd)
         a2 <- h * (fd + fb)
         if(abs(a0 - a1 - a2) < eps || lim == 0)
           return(a1 + a2)
         else {
           return(fun(f, a, d, fa, fd, a1, eps, lim - 1, fun) +
                  fun(f, d, b, fd, fb, a2, eps, lim - 1, fun))
         }
       }
       fa <- f(a)
       fb <- f(b)
       a0 <- ((fa + fb) * (b - a))/2
       fun1(f, a, b, fa, fb, a0, eps, lim, fun1)
     }

Scope: 作用域

這一部分的內容相對本文檔其他部分的內容更偏向一些技術性的問題。 但是它會澄清 S-Plus 和 R 一些重要的 差異。數(shù)據分析師培訓

在函數(shù)內部的變量可以分為三類： 形式參數(shù)，局部變量和自由變量。 形式參數(shù)是出現(xiàn)在函數(shù)的參數(shù)列表中的變量。 它們的值由實際的函數(shù)參數(shù) 綁定形式參數(shù)的過程決定的。 局部變量由函數(shù)內部的表達式的值決定的。 既不是形式參數(shù)又不是局部變量的變量是 自由變量。自由變量 如果被賦值將會變成局部變量?？紤]下面的 函數(shù)定義過程。

     f <- function(x) {
       y <- 2*x
       print(x)
       print(y)
       print(z)
     }

在這個函數(shù)中，x是形式參數(shù)，y是局部變量 ，z是自由變量。

在 R 里面，可以利用函數(shù)被創(chuàng)建的環(huán)境中某個變量的第一次出現(xiàn) 解析一個自由變量的綁定。這稱為 詞法作用域（lexical scope）。我們可以定義一個函數(shù)cube。

     cube <- function(n) {
       sq <- function() n*n
       n*sq()
     }

函數(shù)sq中的變量n不是函數(shù)的參數(shù)。 因此它是自由變量。一些作用域的原則可以用來 確定和它相關的值。在靜態(tài)作用域 (S-Plus)，這個值指的是一個和全局變量n相關的值。在詞法作用域(R)，它指的是函數(shù)cube的參數(shù)。因為當sq定義的時候， 它會動態(tài)綁定參數(shù)n。數(shù)據分析師培訓 在 R 里面解析和在 S-Plus 里面解析不同點在于 S-Plus 搜索 全局變量n而 R 在cube調用時 首先尋找環(huán)境創(chuàng)建的變量n。

     ## 首先用 S 解析 S> cube(2)
     Error in sq(): Object "n" not found
     Dumped
     S> ncube(2)
     [1] 18
     ## 同樣的函數(shù)在 R 中解析 R> cube(2)
     [1] 8

詞匯作用域會給予函數(shù)可變狀態(tài)（mutable state）。 下面的例子演示 R 如何模仿一個銀行的 帳戶。真正的銀行帳戶必須同時有跟蹤收支平衡或者總額的變量， 提供提款業(yè)務，取款業(yè)務和顯示 當前余額的函數(shù)。我們可以在account里面 創(chuàng)建三個函數(shù)，然后返回一個包含它們的 的列表。當調用account時，它讀入一個數(shù)值參數(shù)total，并且返回一個包含三個函數(shù)的列表。 因為這些函數(shù)時定義在一個有變量total的環(huán)境中，它們可以訪問它的值。

<<-是一個特別的賦值操作符，它用來更改和total相關的值。這個操作符會回溯到 一個含有標識符total的密閉環(huán)境中，當它找到這個環(huán)境， 它會用操作符右邊的值替換環(huán)境中這個變量的值。 如果在全局變量或者最高層次的環(huán)境中仍然沒有找到標識符total，那么該變量就會被創(chuàng)建并且在那里被賦值。 大多數(shù)用戶用<<-創(chuàng)建全局變量，并且把操作符右邊 的值賦給它¹。僅僅當<<-用于 一個函數(shù)的輸出是另外一個函數(shù)的輸入時， 這里描述的獨特行為才會出現(xiàn)。

     open.account <- function(total) {
       list(
         deposit = function(amount) {
           if(amount <= 0)
             stop("Deposits must be positive!\n")
           total <<- total + amount
           cat(amount, "deposited.  Your balance is", total, "\n\n")
         },
         withdraw = function(amount) {
           if(amount > total)
             stop("You don't have that much money!\n")
           total <<- total - amount
           cat(amount, "withdrawn.  Your balance is", total, "\n\n")
         },
         balance = function() {
           cat("Your balance is", total, "\n\n")
         }
       )
     }
     
     ross <- open.account(100)
     robert <- open.account(200)
     
     ross$withdraw(30)
     ross$balance()
     robert$balance()
     
     ross$deposit(50)
     ross$balance()
     ross$withdraw(500)

Customizing the environment: 定制環(huán)境

用戶可以有好幾種辦法定制環(huán)境?？梢孕薷?位置初始化文件，并且每個目錄都有它特有的一個 初始化文件。還有就是利用函數(shù).First和.Last。

位置初始化文件的路徑可以通過 環(huán)境變量 R_PROFILE 設置。如果該變量沒有設置， 默認是R安裝目錄下面的子目錄 etc 中的 Rprofile.site。這個文件包括你每次執(zhí)行 R 時一些自動運行的命令。第二個定制文件是 .Rprofile，它可以放在任何目錄下面。如果 R 在該目錄下面 被調用，這個文件就會被載入。這個文件允許用戶 定制它們的工作空間，允許在不同的工作目錄下 設置不同的起始命令。如果在起始目錄中沒有 .Rprofile， R 會在用戶主目錄下面搜索 .Rprofile 文件并且調用它 (如果它 存在的話)。

在這兩個文件或者 .RData 中任何叫.First()的函數(shù) 都有特定的狀態(tài)的。它會在 R 對話的開始時自動執(zhí)行并且初始化環(huán)境。 下面例子中的定義允許 將提示符改為$，以及設置其他有用的東西。 這些設置同樣會在其他會話中起作用。

因此，這些文件的執(zhí)行順序是 Rprofile.site， .Rprofile，.RData 然后是.First()。后面文件中 定義會屏蔽掉前面文件中的定義。

     > .First <- function() {
       options(prompt="$ ", continue="+\t")  #$ is the prompt options(digits=5, length=999)    
     # custom numbers and printout x11() 
                                # for graphics par(pch = "+")  
                      # plotting character source(file.path(Sys.getenv("HOME"), "R", "mystuff.R"))
                                             # my personal functions library(MASS)      
                   # attach a package }

相似的是，如果定義了函數(shù).Last()，它(常常)會在對話 結束時執(zhí)行。一個例子就是

     > .Last <- function() {
       graphics.off()                        # 一個小的安全措施。 cat(paste(date(),"\nAdios\n"))        # 該吃午飯了？ }

Object orientation: 面向對象

一個對象的類決定了它會如何被一個 泛型函數(shù)處理。相反，一個泛型函數(shù) 通過參數(shù)類的特異參數(shù)來完成特定工作或者事務的。 如果參數(shù)缺乏任何類屬性， 或者在該問題中有一個不能被任何泛型函數(shù)處理的類， 泛型函數(shù)會有一種默認的處理方式。

數(shù)據分析師：下面的一個例子使這個問題清晰。類機制為用戶提供了為特定問題設計和編寫 泛型函數(shù)的便利。 在其他泛型函數(shù)中，plot()用于圖形化顯示 對象，summary()用于各種類型的概述分析， 以及anova()用于比較 統(tǒng)計模型。

能以特定方式處理類的泛型函數(shù)的數(shù)目非常的龐大。 例如，可以在非常時髦的類對象"data.frame"中使用的函數(shù)有

     [     [[<-    any    as.matrix
     [<-   mean    plot   summary

可以用函數(shù)methods()得到當前對某個類對象 可用的泛型函數(shù)列表：

     > methods(class="data.frame")

相反，一個泛型函數(shù)可以處理的類同樣很多。 例如，plot()有默認的方法和變量 處理對象類"data.frame"，"density"，"factor"，等等。一個完整的列表同樣可以通過 函數(shù)methods()得到：

     > methods(plot)數(shù)據分析師培訓 

讀者可以參考 完整描述這一機制的正式文檔。