终于轮到介绍 R 的精髓了，那就是包，package。欲知包，必先知函数，function。
我们已经遇到过很多函数了。R 中大多数工作都是函数完成的。经验告诉我们，调用函数的格式是：函数名(变量1 = 某个值，变量2 = 某个值，...)。我们以前用过的，都是 R 基础安装包里有限的一些预先设置好的函数。比如：
x <- 1:5
sd(x) # sd 是函数名, x 是自变量。
## [1] 1.581
sd()用来计算标准差。当我们发出这条指令时，R 在幕后到底忙活个啥呢？输入函数名就能看到了：
sd
## function (x, na.rm = FALSE)
## {
## if (is.matrix(x)) {
## msg <- "sd(<matrix>) is deprecated.\n Use apply(*, 2, sd) instead."
## warning(paste(msg, collapse = ""), call. = FALSE, domain = NA)
## apply(x, 2, sd, na.rm = na.rm)
## }
## else if (is.vector(x))
## sqrt(var(x, na.rm = na.rm))
## else if (is.data.frame(x)) {
## msg <- "sd(<data.frame>) is deprecated.\n Use sapply(*, sd) instead."
## warning(paste(msg, collapse = ""), call. = FALSE, domain = NA)
## sapply(x, sd, na.rm = na.rm)
## }
## else sqrt(var(as.vector(x), na.rm = na.rm))
## }
## <bytecode: 0x0576f9bc>
## <environment: namespace:stats>
甚至 x <- 1:5 这一句，其实背后运行的也是个函数，等同于：
assign("x", 1:5)
内置函数虽然功能强大，但是毕竟有限，要是能随心所欲自己定义新函数就好了。这一点 R 贴心地考虑到了。比如说，当年有一回全班考得都很惨，老师心软了，说为了提高及格率，把卷面分数开方乘十作为新分数吧。为了以后经常用来提高及格率，我们可以专门定义这样一个函数，完全仿照内置函数的格式：
newscore <- function(x) # newscore 是自定义的函数名，它有一个自变量 x。函数的返回值是 sqrt(x) * 10。
{
sqrt(x) * 10
}
以后当考了40分的时候，可以这样调用你的新函数：
newscore(x = 40)
## [1] 63.25
开胃小菜：更妙的操作方式，想一次可以用很久喔！有人說，學電腦的人，動腦筋就是為了偷懶。– 语出大家來學VIM
上面这个例子中，自变量 x 只是用来在函数内部传递信息用的，不会影响函数之外的对象。看看这个例子：
x <- 36
y <- 81
newscore(x = y) # 函数内部的 x 把 81 的值传递进来，而不是36。
## [1] 90
可以有多个自变量：
news <- function(x, n) {
sqrt(x) * 10 + n
}
news(x = 36, n = 10)
## [1] 70
news(36, 10) # 懒人为了省事儿，按自变量的默认顺序写就行了。
## [1] 70
news(n = 10, x = 36) # 如果打乱顺序，就必须指定谁是谁。
## [1] 70
每次调用自定义函数 newscore 的时候，必须提供所有自变量的取值，否则就会报错：
newscore()
## Error: 'x' is missing
为避免这个问题，需要给个默认值：
newscore <- function(x = 36) # x 默认是 36。
{
sqrt(x) * 10
}
newscore()
## [1] 60
给函数定义时，我们用了花括号{}，这意味着可以把一组操作都放进去，哪怕这一组操作有千万行，以后只用一行就可以调用一遍了！比如第 05 篇提到过的指数增长，可以定义一个函数exponentialGrowth：
exponentialGrowth <- function(N0, r = 0.01, tmax = 10) # 三个自变量：初始值，增长率，时间。
{
N <- N0
for (t in 1 : (tmax - 1)) {
N[t + 1] <- N[t] + r * N[t]
}
N # 这是最后一行，作为函数的返回值。
}
exponentialGrowth(66.8)
## [1] 66.80 67.47 68.14 68.82 69.51 70.21 70.91 71.62 72.33 73.06
plot(exponentialGrowth(66.8, 0.01, 100))
练习10.1 自定义一个名为 kaifang 的函数，用来开平方。
练习10.2 自定义一个名为 cv 的函数，用来计算变异系数，即 标准差除以平均值的商。
上面我们自定义了几个函数。世界上很多角落都有想把 36 分变成 60 分的苦命同学，为了让他们也能方便地调用我们自定义的函数，我们可以把它们打包上传到服务器上，这样别人下载了就可以直接用。这就是包。包就是一组预设函数的集合，有时候也包含一些数据。一个包里可能只有一个函数，也可能有成千上万个。R 能走到今天，是一个聚沙成塔、集腋成裘的过程，其中的沙和腋，正是众多热心人花心血写成并奉献出来的扩展包。每个人献出一滴水，终于创造出如今的汪洋大海任你畅游。到底有多少个扩展包呢？用这条命令：
length(unique(rownames(available.packages())))
扩展包是 R 的生命力所在。找到一个合适的扩展包，能起到事半功倍的效果。甚至可以说，会用扩展包，比本文前9篇介绍的所有内容都重要！很多人用 R 就是奔着扩展包来的。
还是举个例子吧。
北京有个著名的广场，常年根据预测的日出日落时间来确定升降国旗仪式的时间，这个是怎么预测的？参看这里。日出日落时间的计算是可以当作新闻来报导的，这让 dapeng 觉得挺神秘。直到有一天，dapeng 需要把某个观测点半年的气温数据（每半小时一条）分为白天和黑夜两组，那么就要判断当地每天的日出和日落时间，不得不设法揭开这个神秘面纱了。上网一搜，哦买告的，计算过程不是一般的复杂啊（见这里），需要有三角函数知识、立体几何知识、天文学知识等等等等，最要命的是还得有足够的耐心。dapeng 花了大概一天的工夫，硬着头皮算出了个数，却跟实际对不上号。
后来，惊喜地发现了 maptools 这个扩展包，安装这个包之后，用其中的 sunriset() 函数一条指令轻松搞定。前天是一个重要的日子，我们计算一下当天的升降国旗时间。
install.packages("maptools") # 第一次使用某个扩展包时要先安装。
require(maptools) # 调用扩展包，让 R 识别其中的函数。
## Loading required package: maptools
## Warning: package 'maptools' was built under R version 2.15.3
## Loading required package: foreign
## Loading required package: sp
## Warning: package 'sp' was built under R version 2.15.3
## Loading required package: grid
## Loading required package: lattice
## Checking rgeos availability: FALSE Note: when rgeos is not available,
## polygon geometry computations in maptools depend on gpclib, which has a
## restricted licence. It is disabled by default; to enable gpclib, type
## gpclibPermit()
position <- c(116.39, 39.91) # 旗杆的经纬度。
mydate <- "2013-03-24" # 要计算的日期。
sunriset(matrix(position, nrow = 1), as.POSIXct(mydate, tz = "Asia/Shanghai"),
direction = c("sunrise"), POSIXct.out = TRUE)$time # 日出时间。
## [1] "2013-03-24 06:11:53 CST"
sunriset(matrix(position, nrow = 1), as.POSIXct(mydate, tz = "Asia/Shanghai"),
direction = c("sunset"), POSIXct.out = TRUE)$time # 日落时间。
## [1] "2013-03-24 18:30:15 CST"
跟官方公布的一点也不差。一个完整的扩展包包括了帮助信息，所以我们的法宝问号仍然管用。自己试试 ‘?sunriset’。若要了解整个扩展包中所有的函数，可以 google 搜索 ‘cran maptools’，也可以在本地计算机 R 的安装路径下面 library 文件夹中找到。
让我们更进一步，自定义一个函数，计算该广场任意一段时期的升降旗时刻：
flag <- function(date.start = "2013-03-24", date.length = 7) # 函数名为flag，默认是计算从前天起一周的升降器时刻。
{
mydate <- seq(as.POSIXct(date.start, tz="Asia/Shanghai"), by = 3600 * 24, length.out = date.length)
data.frame(sunrise = sunriset(matrix(c(116.39, 39.91), nrow = 1), as.POSIXct(mydate, tz="Asia/Shanghai"), direction=c("sunrise"), POSIXct.out = TRUE)$time,
sunset = sunriset(matrix(c(116.39, 39.91), nrow = 1), as.POSIXct(mydate, tz="Asia/Shanghai"), direction=c("sunset"), POSIXct.out = TRUE)$time)
}
flag("2013-10-01") # 好了，以后调用这个函数就能很方便计算。
## sunrise sunset
## 1 2013-10-01 06:10:24 2013-10-01 17:57:17
## 2 2013-10-02 06:11:22 2013-10-02 17:55:40
## 3 2013-10-03 06:12:21 2013-10-03 17:54:03
## 4 2013-10-04 06:13:20 2013-10-04 17:52:27
## 5 2013-10-05 06:14:20 2013-10-05 17:50:51
## 6 2013-10-06 06:15:19 2013-10-06 17:49:16
## 7 2013-10-07 06:16:19 2013-10-07 17:47:41
如果你懒得自己算，只想知道某地的日出日落时刻，请给 dapeng 留言。这里是已经计算的北京等地今后三年的日出日落时刻表。
练习10.3 利用google earth 查出你所在地点的经纬度，然后利用 maptools 扩展包，计算你所在地点 2013 -2113 年 100 年的日出日落时间，然后通知记者来报导“退休职工某某某计算出某地日出日落时间表”。
再举个例子。谢益辉开发了一个制作动画的扩展包，很有趣：
install.packages("animation")
require(animation)
demo("fireworks") # 会用网页浏览器打开一个动画。
citation("animation") # 看看作者。
看着那些琳琅满目的扩展包，dapeng 有种逛苹果商场看app的感觉。网上看到说，“陈红是陈凯歌的第四任妻子，从此之后，陈凯歌再也没有过绯闻。记者采访陈红，拴住大导演的心的秘诀是什么？陈红简单的回答，做百变女人。”扩展包让 R 成了魅力十足的百变女人，教人如何不爱她？
有用的信息：
| – | – |
|---|---|
| 自定义函数 | function() |
| 扩展包 | install.packages(), require(), citation() |
( 连载中，待续 )
yangliufr 网友提了个好建议，希望给代码加上高亮。dapeng 试了一下，下面是个例子。但是 dapeng 懒得一一去改了，多包涵，凑合着看吧。
flag <- function(date.start = "2013-03-24", date.length = 7) # 函数名为flag，默认是计算从前天起一周的升降器时刻。
{
mydate <- seq(as.POSIXct(date.start, tz="Asia/Shanghai"), by = 3600 * 24, length.out = date.length)
data.frame(sunrise = sunriset(matrix(c(116.39, 39.91), nrow = 1), as.POSIXct(mydate, tz="Asia/Shanghai"), direction=c("sunrise"), POSIXct.out = TRUE)$time, sunset = sunriset(matrix(c(116.39, 39.91), nrow = 1), as.POSIXct(mydate, tz="Asia/Shanghai"), direction=c("sunset"), POSIXct.out = TRUE)$time)
}
flag("2013-10-01") # 好了，以后调用这个函数就能很方便计算。