用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)数据分析师培训 

读者可以参考 完整描述这一机制的正式文档。