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

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

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

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

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

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

可以在任何地方以name(expr_1,expr_2, ...)的形式調(diào)用函數(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)

經(jīng)典的 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 提示的一樣，如果調(diào)用的函數(shù)的參數(shù) 以“name=object”的方式給出， 它們可以用任何順序。但是，參數(shù)賦值序列可能以 未命名的，位置特異性的方式給出，同時也有可能 在這些位置特異性的參數(shù)后加上命名參數(shù)賦值。

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

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

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

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

它可以被如下命令調(diào)用

     > ans <- fun1(d, df)

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

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

這就改變了一個默認值。

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

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

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

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

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

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

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

如果想在一個函數(shù)里面全局賦值或者永久賦值，那么可以采用 “強賦值”（superassignment）操作符<<-或者采用函數(shù)assign()。用help命令可以得到更細節(jié)的內(nèi)容。 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ù)據(jù)分析師

去除打印數(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ù)內(nèi)部調(diào)用自己。 但是需要注意的是，這些函數(shù)或者變量， 不會被更高層次的解析框架（evaluation frame）中的被調(diào)用函數(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: 作用域

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

在函數(shù)內(nèi)部的變量可以分為三類： 形式參數(shù)，局部變量和自由變量。 形式參數(shù)是出現(xiàn)在函數(shù)的參數(shù)列表中的變量。 它們的值由實際的函數(shù)參數(shù) 綁定形式參數(shù)的過程決定的。 局部變量由函數(shù)內(nèi)部的表達式的值決定的。 既不是形式參數(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ù)據(jù)分析師培訓 在 R 里面解析和在 S-Plus 里面解析不同點在于 S-Plus 搜索 全局變量n而 R 在cube調(diào)用時 首先尋找環(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ù)，然后返回一個包含它們的 的列表。當調(diào)用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 在該目錄下面 被調(diào)用，這個文件就會被載入。這個文件允許用戶 定制它們的工作空間，允許在不同的工作目錄下 設置不同的起始命令。如果在起始目錄中沒有 .Rprofile， R 會在用戶主目錄下面搜索 .Rprofile 文件并且調(diào)用它 (如果它 存在的話)。

在這兩個文件或者 .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ù)據(jù)分析師：下面的一個例子使這個問題清晰。類機制為用戶提供了為特定問題設計和編寫 泛型函數(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ù)據(jù)分析師培訓 

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