R 菜鸟入门篇 第10篇 函数和包

终于轮到介绍 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))

plot of chunk unnamed-chunk-9

练习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") # 好了,以后调用这个函数就能很方便计算。

R 菜鸟入门篇 第10篇 函数和包》上有46条评论

  1. Yang

    果然趁热打铁,写的非常好啊。代码如果能像Rstudio里一样高亮,对于菜鸟来说,阅读起来可能更好一点。

    回复
    1. dapeng 文章作者

      谢谢鼓励。你的建议很好,以前我考虑过,只是懒得弄,被你这么一说,这回就试了试,用的是WP-Codebox插件,见正文末尾。
      这个系列是用 Rsutdio 的 markdown 写成的,用 knitr 生成了 html 直接拷贝粘贴到博客里来,就成了这模样,原来 html 里的高亮都不见了,猜想可能是 wordpress 的 CSS 的设置导致的?这个我就外行了。

      回复
      1. Yang

        我用的Pretty R的网页代码生成工具,谷歌搜Pretty R就行。可以选择去掉他们的超链接。然后在编辑器里,点源代码编辑之后拷贝上去就好了。

        回复
        1. 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") # 好了,以后调用这个函数就能很方便计算。

          Created by Pretty R at inside-R.org

        2. dapeng 文章作者

          很不错!以后有集中的代码可以这么发布。只是这个方法对这个系列的帖子来说太麻烦了,懒得手动一段一段去改。

        3. dapeng 文章作者

          是的,如果有必要,我打算用 latex 的 listings 宏包来生成高亮代码的 pdf。

  2. 卡卡

    农民工卡卡计算出驻马店1万年日出日落时间表。

    小编们快来采访我吧!

    回复
  3. causu

    dapeng的博文很生动哈,也比较由浅入深,每篇文章,我估计如果dapeng来讲课的话,恰好可以讲一节课,而且非常的有实操性。。。多谢dapeng 。

    回复
  4. glianyi

    在调用library(tseries)时出现Error in library(tseries) : there is no package called ‘tseries’请问怎么解决 刚接触R谢谢你的博客

    回复
  5. vincenhe

    请问楼主知道有木有月亮出没时间的package啊?小弟搜了一下下发觉没有个感觉

    回复
    1. dapeng 文章作者

      这个,好像真没有……不过有个包可以计算月球的位置,改一改应该可以算出月出月落的时刻。

      回复
  6. 夏至

    谢谢前阵子dapeng把最开始练习用的数据文件又特地上传了。今天抽空又补了几课。 没想到animation 的例子点开就是HAPPY NEW YEAR, 过年前还飘在国外自学软件的孩子看到了好惊喜! 祝dapeng一家新年快乐。

    回复
  7. 路今草

    博主,csv某一列的每个单元格里有多个分类变量,而且是用“;”隔开的。
    如何把每个变量提取出来自成一行,其对应那一行的其他数据跟着重复?菜鸟还不会用HTML格式啊,见谅

    回复
      1. 路今草

        额,这么说吧,csv第一列有三个数据(0.1),(0.2),(0.3)。第二列分别对应(beijing)(shanghai;tianjing),(zhenzhou;beijing;shenzhen)。我的问题之怎么用R软件将其变成第一列(0.1),(0.2),(0.2)(0.3)(0.3)(0.3),对应第二列(beijing)(shanghai)(tianjing)(zhenzhou)(beijing)(shehnzhen)

        回复
        1. 大鹏 文章作者

          你的意思是,第一列是括号分隔的数据?0.2要重复2次,0.3要重复3次?依次跟六个城市对应?我觉得你还是直接把数据贴上来吧。多贴几行。

        2. 路今草

          columnA columnB
          beijing 0.5454
          beijing;shanghai 0.3536

  8. 路今草

    不是括号分割,是用的分号;确实根据第二列里的数据个数,重复第一列的内容

    回复
    1. 大鹏 文章作者

      明白了。这个容易。先用字符串函数得到columnA里分号的数量n,然后让columnB重复n+1次,就可以对应上了,再生成个新列即可。你可以把数据发给我,我写个代码给你,一看就懂了。

      回复
    2. 大鹏 文章作者

      我把你的数据截取了一小段,放在了网上,方便代码读取和讨论,希望你不会介意。下面是代码。

      data <- read.csv(file = 'http://pastebin.com/raw.php?i=vLL11Hr6', stringsAsFactors = FALSE)
      newlist <- strsplit(data$UCSC_RefGene_Group, ";")
      newdata <- data.frame(unlist(newlist))
      names(newdata) <- names(data)[1]
      comma_nr <- lapply(newlist, length)
      data$rep <- unlist(comma_nr)
      newdata$v <- rep(data$X11N, data$rep)
      newdata
      
      回复
      1. 路今草

        非常感谢,数据用没有问题,以后可能还要多向您请教

        回复
      2. oocainiao

        LZ您好,首先感谢您的分享,很受用。针对“路今草”的问题,我写了一段
        data <- read.csv(file.choose())
        data

        m <- data.frame()
        for( i in 1:length(data$columnB)) {
        strA <- toString(data$columnA[i])
        list <- unlist(strsplit(strA,";"))
        if(is.null(list)) {
        m <- cbind(strA,toString(data$columnB[i]))
        } else {
        m <- cbind(list,(data$columnB[i]))
        }
        }
        m

        问题是会出现数据被覆盖的问题,它只显示最后两条,这是为什么

        回复
        1. 大鹏 文章作者

          m <- cbind(list,(data$columnB[i]))
          

          改成

          m <- rbind(m, cbind(list,(data$columnB[i])))
          
  9. 大熊

    sd这个内置函数的完整定义要怎么才能看到啊?

    我在rstudio里面运行sd结果如下。说明文档里面也没有详细的定义

    sd
    function (x, na.rm = FALSE)
    sqrt(var(if (is.vector(x)) x else as.double(x), na.rm = na.rm))

    回复
  10. hnndylf

    很受用,我是你同学的学生,不知还有没有这种简单名了的资料。

    回复

发表评论

电子邮件地址不会被公开。 必填项已用*标注

To create code blocks or other preformatted text, indent by four spaces:

    This will be displayed in a monospaced font. The first four 
    spaces will be stripped off, but all other whitespace
    will be preserved.
    
    Markdown is turned off in code blocks:
     [This is not a link](http://example.com)

To create not a block, but an inline code span, use backticks:

Here is some inline `code`.

For more help see http://daringfireball.net/projects/markdown/syntax