用R語言進行數據分析:編寫函數

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

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

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

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

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

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

Simple examples: 簡單的例子

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

函數定義如下：

     > 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
       }

通過這個函數，你可以下面的命令 實現雙樣本 t-檢驗。

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

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

給定 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)
     }

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

     > regcoeff <- bslash(Xmat, yvar)

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

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

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

     %anything%

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

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

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

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

Named arguments and defaults: 參數命名和默認值

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

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

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

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

     > 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)

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

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

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

它可以被如下命令調用

     > ans <- fun1(d, df)

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

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

這就改變了一個默認值。

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

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

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

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

Assignment within functions: 函數內部的賦值

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

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

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

More advanced examples: 高級的例子

區組設計中的效率因子

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

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

     > 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)
     }

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

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

去除打印數組中的名字

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

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

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

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

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

     > no.dimnames(X)

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

遞歸式的數值積分

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

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

這個例子同樣可以部分的作為一個 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 一些重要的 差異。數據分析師培訓

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

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

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

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

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

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

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

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

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

     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: 定制環境

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

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

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

因此，這些文件的執行順序是 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 }

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

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

Object orientation: 面向對象

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

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

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

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

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

     > methods(class="data.frame")

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

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

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