R/R语法

R的对象数据类型与数据结构

R操作的实体在技术上来说都是对象(object)。当R在运行时，所有变量，数据，函数及结果都以对象的形式存在计算机的活动内存中，并有相应的名字对应。

R的对象

R有五种基本（atomic）的数据类型：

• 字符串 character
• 数值型（实数）numeric
• 整数 integer
• 复合型 complex
• 逻辑型 logistic (True/False)

R最基本的对象是向量

• 向量（vertor）只能储存一种数据类型
• 列表（list）是个特例，用向量表示，可以存储不同数据类型的对象 vector() 函数可以创建一个空向量

v <- vector()

数值

• 数值在R中储存为numeric对象（例如双精度实数）
• 如果想明确的声明一个整数，添加L 后缀

> class(1)
[1] "numeric"
> class(1L)
[1] "integer"

• 特殊数值 Inf 表示无限infinity，Inf也可以用于计算正常计算

> 1/ 0
[1] Inf
> 1 / Inf
[1] 0

• NaN代表未定义的值 (“not a number”)，可以看作一个缺失值

> 0/0
[1] NaN

属性 Attributes

R对象有不同的属性

• 名称，names；行列名，dimnames
• 维度，dimensions (例如矩阵matrix，数组array)
• 类型，class
• 长度，length
• 模式，mode
• 其他自定义的属性attributes/metadata attributes()函数可以查看一个对象的属性。

> x <-matrix(0,4,5)
> x
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    0    0    0    0
[2,]    0    0    0    0    0
[3,]    0    0    0    0    0
[4,]    0    0    0    0    0
> class(x)
[1] "matrix"
> mode(x)
[1] "numeric"
> length(x)
[1] 20

输入

<- 符号进行赋值操作。

> x <- 1
> print(x)
[1] 1
> x
[1] 1
> msg <- "hello"
> msg
[1] "hello"

编程语言的语法决定表达式是否完成。

## 未完成的表达式
> x <- 
# #字符开头表示注释，#字符右边的内容被忽略。

赋值计算Evaluation

输入一个完整的表达式，R语言会计算表达式，并返回结果。

> x <- 5  ## 不打印
> x       ## 自动打印
[1] 5
> print(x)  ## 用print来输出结果
[1] 5

输出 Printing

> x <- 1:20 
> x
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
[16] 16 17 18 19 20

:操作符用来创建整数序列。

R向量vector和列表list

创建向量vector

c()函数可以用来创建向量对象。

> x <- c(0.5, 0.6)       ## numeric
> x <- c(TRUE, FALSE)    ## logical
> x <- c(T, F)           ## logical
> x <- c("a", "b", "c")  ## character
> x <- 9:29              ## integer
> x <- c(1+0i, 2+4i)     ## complex

可以使用vector()函数创建向量

> x <- vector("numeric", length = 10) 
> x
 [1] 0 0 0 0 0 0 0 0 0 0

向量中创建不同数据类型的对象

# 将数值和字符声明到同一向量对象
> y <- c(1.7, "a") 
# 数值强制转换成字符串
> mode(y)
[1] "character"
# 将逻辑值和数值声明到一个向量对象中
> y <- c(TRUE, 2) 
# 逻辑型强制转换为数值型
> mode(y)
[1] "numeric"
# 在一个向量对象中储存字符串和逻辑值
> y <- c("a", TRUE) 
# 逻辑值强制转换为字符串
> mode(y)
[1] "character"

当在同一个向量对象中赋值不同数据类型属，coercion强制转换向量中的每个元素成同一种数据类型。

指定强制转换类型

通过as.*函数，可以将对象转换为指定的数据类型

# 生成一个整数型向量x
> x <- 0:6
# 查看x的数据类型
> class(x)
[1] "integer"
# 转换成数值型
> as.numeric(x)
[1] 0 1 2 3 4 5 6
# 转换成逻辑型
> as.logical(x)
[1] FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
# 转换成字符型
> as.character(x)
[1] "0" "1" "2" "3" "4" "5" "6"

无意义的转换产生NA

# 生成一个字符串型向量x
> x <- c("a", "b", "c")
# 强制转换为数值型产生NA
> as.numeric(x)
[1] NA NA NA
Warning message:
NAs introduced by coercion
# 强制转换为逻辑型产生NA
> as.logical(x)
[1] NA NA NA
# 强制转换为复杂性型向量
> as.complex(x)
[1] 0+0i 1+0i 2+0i 3+0i 4+0i 5+0i 6+0i

矩阵matrix和列表list

矩阵matrix

矩阵是具有二维dimension  (nrow, ncol)维度属性的向量。

> m <- matrix(nrow = 2, ncol = 3) 
> m
     [,1] [,2] [,3]
[1,]   NA   NA   NA
[2,]   NA   NA   NA
> dim(m)
[1] 2 3
> attributes(m)
$dim
[1] 2 3
> dim(m)

矩阵通过扩展列column-wise构建，生成矩阵时，输入从左上角开始，顺列填充。

> m <- matrix(1:6, nrow = 2, ncol = 3) 
> m
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

矩阵可以通过向量添加维度属性来创建。

> m <- 1:10 
> m
[1] 1 2 3 4 5 6 7 8 9 10 
> dim(m) <- c(2, 5)
> m
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    3    5    7    9
[2,]    2    4    6    8   10

cbind和rbind

矩阵可以通过列合并column-binding函数cbind()及行合并row-binding函数rbind()来创建。

> x <- 1:3
> y <- 10:12
> cbind(x, y)
     x  y 
[1,] 1 10 
[2,] 2 11 
[3,] 3 12
> rbind(x, y) 
  [,1] [,2] [,3]
x    1    2    3
y   10   11   12

列表list

列表list是可以包含不同数据类型的特殊向量。列表是R中非常重要的数据类型。

> x <- list(1, "a", TRUE, 1 + 4i) 
> x
[[1]]
[1] 1

[[2]] 
[1] "a"

[[3]]
[1] TRUE

[[4]]
[1] 1+4i

因子Factors

因子用来表示数据分组，可以有序或无序。因子可以视作带标签label的整型向量。

• 因子可通过模型函数lm() 和glm()创建
• 使用标签的因子分组比使用数据更直观，“Male”和“Female”分组比1和2分组更易于理解

> x <- factor(c("yes", "yes", "no", "yes", "no")) 
> x
[1] yes yes no yes no
Levels: no yes
> table(x) 
x
no yes 
 2   3
> unclass(x)
[1] 2 2 1 2 1
attr(,"levels")
[1] "no"  "yes"

水平的先后顺序可以使用factor()函数的levels参数来指定。这在线性模型中很重要，因为第一个水平常被用来作为对照。

> x <- factor(c("yes", "yes", "no", "yes", "no"),
              levels = c("yes", "no"))
> x
[1] yes yes no yes no 
Levels: yes no

缺失值Missing Values

缺失值NA和NaN代表未定义的数学操作。

• is.na()函数用来测试R对象是否为NA
• is.nan()函数用来测试R对象是否为NaN
• NA值也有数据类型，整型NA， 字符串character NA等等
• NaN值属于NA值

> x <- c(1, 2, NA, 10, 3)
> is.na(x)
[1] FALSE FALSE  TRUE FALSE FALSE
> is.nan(x)
[1] FALSE FALSE FALSE FALSE FALSE
> x <- c(1, 2, NaN, NA, 4)
> is.na(x)
[1] FALSE FALSE  TRUE  TRUE FALSE
> is.nan(x)
[1] FALSE FALSE  TRUE FALSE FALSE

数据框dataframe与R对象的命名

数据框dataframe

数据框用来储存表格型的数据

• 数据框是一种特殊类型的多维列表，其中每个列表的长度必须相等
• 一个列表就是数据框的一列，列表的长度就是数据框的行数
• 与矩阵不同，数据框可以在一列中储存不同类型的对象（与list类似）；矩阵中的每个元素都必须是相同的数据类型
• 数据框有一个特殊的属性行名row.names
• 数据框一般通过调用read.table()或read.csv()函数构建
• 数据框可以通过调用data.matrix()函数转换为矩阵

> x <- data.frame(foo = 1:4, bar = c(T, T, F, F)) 
> x
  foo   bar
1   1  TRUE
2   2  TRUE
3   3 FALSE
4   4 FALSE
> nrow(x)
[1] 4
> ncol(x)
[1] 2

R对象的命名

对R对象进行命名，可以提升代码的可读性。

> x <- 1:3
> names(x)
NULL
> names(x) <- c("foo", "bar", "norf") 
> x
foo bar norf 
  1   2    3
> names(x)
[1] "foo"  "bar"  "norf"

对列表list命名。

> x <- list(a = 1, b = 2, c = 3) 
> x
$a
[1] 1

$b 
[1] 2

$c 
[1] 3

矩阵matrix的命名。

> m <- matrix(1:4, nrow = 2, ncol = 2)
> dimnames(m) <- list(c("a", "b"), c("c", "d")) 
> m
  c d 
a 1 3 
b 2 4

R数据类型与结构小结

数据结构是计算机存储、组织数据的方式。数据结构是指相互之间存在一种或多种特定关系的数据元素的集合。

• 基本类：数值numeric（用于直接计算、加减乘除）, 逻辑logical（真假）, 字符串character（连接、转换、提取），整型integer，复合型complex，日期等
• 向量vector, 列表list
• 因子factor
• 缺失值missing values
• 数据框data frames
• 命名names

常见的数据类型及结构

对象	模式	是否允许同一个对象中有多种模式
向量	数值型，字符型，复数型，逻辑型	否
因子	数值型，字符型	否
数组	数值型，字符型，复数型，逻辑型	否
矩阵	数值型，字符型，复数型，逻辑型	否
数据框	数值型，字符型，复数型，逻辑型	是
时间序列	数值型，字符型，复数型，逻辑型	否
列表	数值型，字符型，复数型，逻辑型，函数，表达式，…	是

R文件操作

获取数据的方式

• 利用键盘输入数据
• 读取存储在磁盘上的数据文件
• 通过开放数据库连接（ODBC, Open database Connectivity）访问数据库系统来获取数据(RODBC包)

# 可以使用scan函数一个个输入来创建向量
> patientID <- scan()
1: 1
2: 2
3: 3
4: 4
5: 
Read 4 items
# 也可以直接使用c()函数创建向量，结果是相同的
> patientID <- c(1, 2, 3, 4)
> admdate <- c("10/15/2009", "11/01/2009", "10/21/2009", "10/28/2009")
> age <- c(25, 34, 28, 52)
> diabetes <- c("Type1", "Type2", "Type1", "Type1")
> status <- c("Poor", "Improved", "Excellent", "Poor")
> data <- data.frame(patientID, age, diabetes, status)
> data
  patientID age diabetes    status
1         1  25    Type1      Poor
2         2  34    Type2  Improved
3         3  28    Type1 Excellent
4         4  52    Type1      Poor


data2 <- data.frame(patientID=character(0), age=numeric(0), diabetes=character(), status=character())
# 使用edit和fix函数来添加数据
> data2 <- edit(data2)
> fix(data2)

其他R编辑数据的函数，用法与edit和fix类似：

vi(name = NULL, file = "")
emacs(name = NULL, file = "")
pico(name = NULL, file = "")
xemacs(name = NULL, file = "")
xedit(name = NULL, file = "")

读取数据

R中有很多函数可以读取数据

• read.table, read.csv用来读取制表符分隔的数据
• readLines 读取文本文件的一行
• source 读取R代码文件 (与dump对应)
• dget 读取R代码文件 (与dput对应
• load 读取保存的工作空间
• unserialize 读取一个serialize封装的二进制R对象

写入文件

R中有很多功能类似的函数将数据写入文件

• write.table
• writeLines
• dump
• dput
• save
• serialize

使用read.table函数读取数据

read.table函数是最长用的读取数据的函数。它有几个非常重要的参数：

file 文件的名字或者链接
header 逻辑值表示文件是否有表头
sep 分隔符，指定列分隔符是什么
colClasses 字符串向量表示数据集中每列的类
nrows 数据集的行数
comment.char 字符串表示注释
skip 忽略开始的多少行
stringsAsFactors 逻辑值，设置字符串变量是否用因子结构储存

对于中小型数据集，使用read.table函数的默认参数即可

data<-read.table("foo.txt")

• R会自动忽略以“#”开头的行
• 设置数据集的行数（需要占用多少内存）和每列数据的变量类型可以加快R的运行效率f read.csv函数除了默认分隔符是逗号以外，其他参数和read.table函数相同。

使用read.table函数读取大数据集

对于大型数据集，首先大致估计一下需要使用多少内存来储存数据集，如果数据所需要的内存大于计算机的物理内存，那就需要使用其他方式处理数据。

• 如果数据集中没有注释行将参数comment.char = ""
• 使用colClasses参数可以提高一倍的读取速度，这要求指定数据集中每一列的数据类型。如果每列都是数值，可以直接设置参数colClasses = "numeric"。另一个快速设置每列的类型的方法如下：

initial<-read.table("datatable.txt",nrows=100)
classes<-sapply(initial,class)
tabAll<-read.table("datatable.txt",colClasses=classes)

• 设置nrows参数。不能提速当能够帮助R优化内存的使用，可以使用Linux的wc命令来计算数据集的行数。

了解一些操作系统

使用R处理大数据集时，需要了解一些操作系统的情况。

• 电脑的内存多大
• 有什么应用程序在运行
• 是否有其他用户登陆到该操作系统
• 什么操作系统
• 32还是64位操作系统

计算内存使用

一个1,500,000行和120列的数值型数据框占用内存的计算： 1,500,000 × 120 × 8 bytes/numeric = 1440000000 bytes = 1440000000 / ​ bytes/MB = 1,373.29 MB = 1.34 GB

文本格式

• dump和dput函数的是文本格式，可以直接编辑，在文件出错时，可以恢复。
• dump和dput保存了对象的metadata（R对象的属性信息），这降低了结果的可读性，但读入数据时，不需要在对该数据集进行额外的处理。
• 文本格式文件方便进行版本控制，git等命令只能追踪文本文件的改变。
• 文件格式用处更广泛，如果文本中出现错误，相比二进制文件，更容易被修复。
• 文本格式更符合Unix哲学
• 缺点：占用更多的空间（相比于二进制文件）

dput操作R对象

通过dput函数将R对象写入文件，使用dget函数将文件读入R。

> y <- data.frame(a = 1, b = "a")
> dput(y)
structure(list(a = 1,
               b = structure(1L, .Label = "a",
                             class = "factor")),
          .Names = c("a", "b"), row.names = c(NA, -1L),
          class = "data.frame")
> dput(y, file = "y.R")
> new.y <- dget("y.R")
> new.y
   a  b 
1  1  a

dump函数操作R对象

dump函数可以同时打包多个R对象写入文件，可以使用source函数将文件读回R。

> x <- "foo"
> y <- data.frame(a = 1, b = "a")
> dump(c("x", "y"), file = "data.R") 
> rm(x, y)
> source("data.R")
> y
  a  b 
1 1  a
> x
[1] "foo"

其他数据操作的函数

数据操作可以通过创建一个链接（其他语言称为文件句柄）来实现。

• file 创建一个普通文件的链接
• gzfile 创建一个gzip格式压缩文件的链接
• bzfile 创建一个bzip2格式压缩文件的链接
• url 创建一个读取网页的链接

文件链接

> str(file)
function (description = "", open = "", blocking = TRUE,
          encoding = getOption("encoding"))

• description参数设置读取文件的名称
• open参数设置文件的操作权限
  • “r” 只读
  • “w” 可写并初始化一个新文件
  • “a” 追加
  • “rb”, “wb”, “ab” 读、写、追加二进制格式文件

文件链接可以帮助用户更灵活地操作文件或外部对象，我们不用直接对该链接进行操作，R中的函数可以直接处理。

con <- file("foo.txt", "r")
data <- read.csv(con)
close(con)

等同于

data <- read.csv("foo.txt")

按行读取文件

> con <- gzfile("words.gz") 
> x <- readLines(con, 10) 
> x
 [1] "1080"     "10-point" "10th"     "11-point"
 [5] "12-point" "16-point" "18-point" "1st"
 [9] "2"        "20-point"

writeLines读取字符串向量，将向量中每行对象按行写入文件。

---

readLines在读取网页数据十分有用

> con <- url("http://www.baidu.com", "r")
> x <- readLines(con)
> head(x)
[1] "<!DOCTYPE html>"  "<!--STATUS OK-->" ""                 ""                
[5] ""

提取R对象的子集

R中有很多操作可以提取R对象的子集。

• [ 操作返回跟原R对象类相同的对象；可以同时提取多个对象（例外）
• [[ 用来提取列表list或数据框dataframe对象；该操作只能用来提取单个元素，返回的对象不一定是list或数据框
• $ 通过list和dataframe的命名name来提取元素，语义上类似于 [[

> x <- c("a", "b", "c", "c", "d", "a")
> x[1]
[1] "a"
> x[2]
[1] "b"
> x[1:4]
[1] "a" "b" "c" "c"
> x[x > "a"]
[1] "b" "c" "c" "d"
> u <- x > "a"
> u
[1] FALSE  TRUE  TRUE  TRUE  TRUE FALSE
> x[u]
[1] "b" "c" "c" "d"

提取矩阵的子集

矩阵的子集的提取通常通过 (i,j)来索引。

> x <- matrix(1:6, 2, 3)
> x
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> x[1, 2]
[1] 3
> x[2, 1]
[1] 2

索引可以缺失。

> x[1, ]
[1] 1 3 5
> x[, 2]
[1] 3 4

当提取矩阵的单个元素时，默认返回一个长度为1的向量，而不是一个1 × 1的矩阵。可以通过设置参数drop = FALSE来改变。

> x <- matrix(1:6, 2, 3)
> x[1, 2]
[1] 3
> x[1, 2, drop = FALSE]
     [,1]
[1,]    3

类似的，提取矩阵的一列或一行时默认返回向量，而不是矩阵。

> x <- matrix(1:6, 2, 3)
> x[1, ]
[1] 1 3 5
> x[1, , drop = FALSE]
     [,1] [,2] [,3]
[1,]    1    3    5

提取列表的子集

第一个例子

> x<-list(foo = 1:4, bar = 0.6)
> x[1]
$foo
[1] 1 2 3 4

> x[[1]]
[1] 1 2 3 4
> x$bar
[1] 0.6
> x[["bar"]]
[1] 0.6
> x["bar"]
$bar
[1] 0.6

第二个例子

> x <- list(foo = 1:4, bar = 0.6, baz = "hello")
> x[c(1, 3)]
$foo
[1] 1 2 3 4

$baz
[1] "hello"

[[ 可以使用变量或表达式作为索引；$ 只能按照命名的名称来索引。

> x <- list(foo = 1:4, bar = 0.6, baz = "hello")
> name <- "foo"
## 使用变量索引 ‘foo’
> x[[name]]
[1] 1 2 3 4
## 不存在元素名称为 ‘name’
> x$name
NULL
## ‘foo’存在
> x$foo
[1] 1 2 3 4

提取嵌套列表（多维维列表）的子集

[[ 可以使用整数序列作为索引。

> x <- list(a = list(10, 12, 14), b = c(3.14, 2.81))
> x[[c(1, 3)]]
[1] 14
> x[[1]][[3]]
[1] 14
> x[[c(2, 1)]]
[1] 3.14

局部匹配

[[ 和$允许名称的局部匹配来索引子集。

> x <- list(aardvark = 1:5)
> x$a
[1] 1 2 3 4 5
> x[["a"]]
NULL
> x[["a", exact = FALSE]]
[1] 1 2 3 4 5

去除NA值

在进行数据分析的过程中，经常需要对缺失值(NA)进行处理。

> x <- c(1, 2, NA, 4, NA, 5)
> bad <- is.na(x)
> x[!bad]
[1] 1 2 4 5

去除多个R对象的缺失值，并提取不含缺失值的子集。 第一个例子

> x <- c(1, 2, NA, 4, NA, 5)
> y <- c("a", "b", NA, "d", NA, "f")
> good <- complete.cases(x, y)
> good
[1]  TRUE  TRUE FALSE  TRUE FALSE  TRUE
> x[good]
[1] 1 2 4 5
> y[good]
[1] "a" "b" "d" "f"

第二个列子

> airquality[1:6, ]
  Ozone Solar.R Wind Temp Month Day
1    41     190  7.4   67     5   1
2    36     118  8.0   72     5   2
3    12     149 12.6   74     5   3
4    18     313 11.5   62     5   4
5    NA      NA 14.3   56     5   5
6    28      NA 14.9   66     5   6
> good <- complete.cases(airquality)
> airquality[good, ][1:6, ]
  Ozone Solar.R Wind Temp Month Day
1    41     190  7.4   67     5   1
2    36     118  8.0   72     5   2
3    12     149 12.6   74     5   3
4    18     313 11.5   62     5   4
7    23     299  8.6   65     5   7
8    19      99 13.8   59     5   8

向量化操作

R中的很多操作都是针对向量的，这也使得R代码高效、简洁、可读性强。

> x <- 1:4; y <- 6:9 
> x
[1] 1 2 3 4
> y
[1] 6 7 8 9
> x + y 
[1]  7  9 11 13
> x > 2
[1] FALSE FALSE  TRUE  TRUE
> x >= 2
[1] FALSE  TRUE  TRUE  TRUE
> y == 8
[1] FALSE FALSE  TRUE FALSE
> x * y
[1]  6 14 24 36
> x / y
[1] 0.1666667 0.2857143 0.3750000 0.4444444

矩阵的向量化操作

> x <- matrix(1:4, 2, 2); y <- matrix(rep(10, 4), 2, 2)
> x
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> y
     [,1] [,2]
[1,]   10   10
[2,]   10   10
## element-wise乘法（数组元素依次相乘）
> x * y
     [,1] [,2]
[1,]   10   30
[2,]   20   40
> x / y
     [,1] [,2]
[1,]  0.1  0.3
[2,]  0.2  0.4
## 真正的矩阵乘法
> x %*% y
     [,1] [,2]
[1,]   40   40
[2,]   60   60

R控制结构

R中的控制结构允许用户在不同情况下设置程序的执行流程。常见的结构包括：

• if, else: 判断条件
• for: 执行固定次数的循环
• while: 当 while 的条件为真时，执行循环
• repeat: 执行无限次循环
• break: 退出整个循环
• next: 跳过本次循环
• return: 退出函数 大多数控制结构一般不在R的交互命令对话中使用，主要应用在自定义函数或较长的表达式中。

控制结构：if

if(<condition>) {
        ## do something
} else {
        ## do something else
}
if(<condition1>) {
        ## do something
} else if(<condition2>)  {
        ## do something different
} else {
        ## do something different
}

下面是一些有效的 if/else 结构示例。

if(x > 3) {
        y <- 10
} else {
        y <- 0
}

等同于

y <- if(x > 3) {
        10
} else { 
        0
}

else语句不是必须的。

if(<condition1>) {

}

if(<condition2>) {

}

for

for 循环接收可迭代的变量，一般是来自序列或向量的连续值。for循环经常使用的循环元素来自于列表、向量等R对象。

for(i in 1:10) {
         print(i)
}

该循环将变量 i 作为循环变量，循环依次赋值1, 2, 3, ..., 10，然后退出循环。 下面几个循环的功能相同。

x <- c("a", "b", "c", "d")

for(i in 1:4) {
        print(x[i])
}

for(i in seq_along(x)) {
        print(x[i])
}

for(letter in x) {
        print(letter)
}

for(i in 1:4) print(x[i])

嵌套循环

for 循环可以嵌套使用。

x <- matrix(1:6, 2, 3)

for(i in seq_len(nrow(x))) {
        for(j in seq_len(ncol(x))) {
                print(x[i, j])
        }   
}

慎重使用嵌套循环，嵌套超过2-3层后会难以阅读和理解，也会拖慢程序的运行速度。

while

while循环开始前先测试循环条件。如果值为真，执行循环体。当循环体运行完，再次测试循环条件，如此循环下去。

count <- 0
while(count < 10) {
        print(count)
        count <- count + 1
}

while循环如果没有正确编写，可能会陷入无限循环（死循环），使用时要注意。

---

有时需要同时满足多个判断条件进行循环。

z <- 5

while(z >= 3 && z <= 10) {
        print(z)
        # 随机生成0或1
        coin <- rbinom(1, 1, 0.5) 

        if(coin == 1) { 
                z <- z + 1
        } else {
                z <- z - 1
        } 
}

判断条件是从左到右进行判断。

---

repeat

repeat初始化一个无限循环，有特殊用途，通常不应用在数据统计分析中，退出repeat循环的唯一方法是break。 以下循环是个repeat的无限循环，运行后需要强制终止。

x0 <- 1
tol <- 0

repeat {
        x1 <- rbinom(1, 1, 0.5) 

        if(abs(x1 - x0) < tol) {
                break
        } else {
                x0 <- x1
        } 
}

上面的循环对于初学者不友好，容易陷入死循环。最好设置迭代次数限制（例如使用for循环），然后打印是否得到预期结果。

---

next, return

next用来跳过本次循环。

for(i in 1:100) {
        if(i <= 20) {
                ## 跳过前20次循环
                next 
        }
        ## Do something here
}

return 表示一个函数的结果同时返回函数的运算结果。

---

小结

• R语言中if、while和for 循环和判断结构。
• 避免死循环，即使在理论上是合理的。
• 控制结构语句对R编程十分有用；对于R的交互式命令行操作，apply家族的函数更高效。

函数

通过function()可以自定义函数，函数的储存与其他的R对象相同，函数对象的类为 “function”。

f <- function(<参数>) {
        ## R语句
}

R函数也是R对象，可以向其他R对象一样进行操作：

• 函数可以作为参数传递给其他函数
• 函数可以嵌套，可以在函数中包含函数。
• 函数的返回值是函数体中最后一个表达式。
  • *

函数的参数

函数的 参数 一般会有默认值。

• 形式参数指的是在函数内定义的参数
• formals 函数可以查看一个函数所有的形式参数
• 不是每个R函数都需要使用所有的参数
• 函数的参数可以省略或保持默认值。

参数匹配

R函数参数可以通过position matching或name matching进行匹配。以下调用 sd 函数的操作都是等同的。

> mydata <- rnorm(100)
> sd(mydata)
[1] 0.9418245
> sd(x = mydata)
[1] 0.9418245
> sd(x = mydata, na.rm = FALSE)
[1] 0.9418245
> sd(na.rm = FALSE, x = mydata)
[1] 0.9418245
> sd(na.rm = FALSE, mydata)
[1] 0.9418245

虽然这些操作的结果都相同，但为避免造成误解，一般按照将要函数参数的默认顺序来赋值。

---

可以混合使用position matching和name matching。name matching的优先级高于position matching。

> args(lm)
function (formula, data, subset, weights, na.action, method = "qr", 
    model = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE, 
    contrasts = NULL, offset, ...) 
NULL

下面两种操作等同。

lm(data = mydata, y ~ x, model = FALSE, 1:100)
lm(y ~ x, mydata, 1:100, model = FALSE)

• 命名参数在交互式命令行中使用较多，尤其是参数较多时。
• 命名参数可以帮助记忆函数的参数名称，这在绘图中很有用。

函数参数可以部分匹配，有助于交互式操作。参数匹配的优先级的顺序如下：

1. 检查命名参数的精确匹配Check for exact match for a named argument
2. 检查局部匹配
3. 检查位置匹配

---

定义函数

f <- function(a, b = 1, c = 2, d = NULL) {

}

如果不想对函数的某个参数设置默认值，可以将 NULL 赋值给该参数。

惰性求值（延后计算）

R函数的参数的计算遵循惰性原则，只在需要时进行计算。

f <- function(a, b) {
    a^2
}
f(2)

## [1] 4

该函数实际没有调用参数 b ，因此调用 f(2) 不会报错，根据位置匹配，将2赋值给2 a。

f <- function(a, b) {
    print(a)
    print(b)
}
f(45)

## [1] 45
## Error: argument "b" is missing, with no default

函数首先打印了“45”，然后报错，因为在执行 print(a) 时，参数 b 未调用。执行 print(b) 时程序报错

---

“...” 参数

... 参数用来传递数目不确定的参数给其他函数。

• ... 常被用来扩展一个函数，避免原函数的所有参数列表再次声明

myplot <- function(x, y, type = "l", ...) {
        plot(x, y, type = type, ...)
}

• 一般函数使用 ... 传递额外的参数给方法（Method）

> mean
function (x, ...)
UseMethod("mean")

如果无法提前知道需要传递给函数的参数数目，需要使用“...”参数。

> args(paste)
function (..., sep = " ", collapse = NULL)

> args(cat)
function (..., file = "", sep = " ", fill = FALSE,
    labels = NULL, append = FALSE)

“...” 参数后的参数

“...” 参数后的参数必须通过名称精确匹配，不能局部匹配。

> args(paste)
function (..., sep = " ", collapse = NULL)

> paste("a", "b", sep = ":")
[1] "a:b"

> paste("a", "b", se = ":")
[1] "a b :"

变量的作用域

变量与变量名

对相同的变量名多次赋值时，R如何识别该变量名所对应的值？ 一个例子：

> lm
function (formula, data, subset, weights, na.action, method = "qr", 
    model = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE, 
    contrasts = NULL, offset, ...) 
{......
}
<bytecode: 0x7feae2a23638>
<environment: namespace:stats>
> lm <- function(x) { x * x }
> lm
function(x) { x * x }

当对 lm 赋值后，打印的值发生了变化，不再返回 stats 包中的 lm 所对应的值。

---

当对一个变量赋值时，R会搜索一系列的环境。当R取值时，环境的优先级从高到低如下：

1. 搜索全局环境中该变量名对应的值。
2. 搜索每个扩展包中命名空间中的变量列表。

search 函数可以查看R的变量搜索列表。

> search()

[1] ".GlobalEnv"        "package:stats"     "package:graphics"

[4] "package:grDevices" "package:utils"     "package:datasets"

[7] "package:methods"   "Autoloads"         "package:base"

• 全局变量global environment或用户的工作空间在列表的最前面，base扩展包排在最后。
• 扩展包的前后顺序很重要。
• 用户可以设置R启动时加载的扩展包。
• 使用 library 函数加载扩展包时，该扩展包的命名空间会默认排在搜索列表的第二位，其他变量列表依次后移。
• R中函数与非函数的命名空间是分离的，可以同时存在名称为c的R对象和函数。

变量的作用域

R语言的作用域规则与S语言截然不同。

• 作用域规则决定函数中值如何与自由变量的关联。
• R使用的是静态作用域，一般称为lexical scoping，与之对应的是动态作用域。
• 作用域规则就是R在如何在变量列表中搜寻变量的规则。
• Lexical scoping有助于简化统计计算。

Lexical Scoping

一个示例：

f <- function(x, y) {
        x^2 + y / z
}

该函数有两个形式参数x 和 y，函数体内的 z 就是一个自由变量。编程语言的作用域规则决定值如何赋值给自由变量。自由变量既不是形式参数，也不是本地变量（在函数体内声明的变量）。

---

R的作用域规则：在定义函数的环境中搜索自由变量的值

环境：

• 环境就是一系列对应的变量名符号和值，如 x 是变量名符号，3.14是它对应的值。
• 每个环境都有一个母环境；一个环境可以有多个子环境。
• 没有母环境的环境是空环境。
• 函数 + 环境 = 闭包 或 函数闭包（函数可以使用函数之外定义的自由变量）。

搜索自由变量的值：

• 如果在函数定义的环境中没有找到变量名对应的值，接着搜索母环境的变量。
• 依次搜索环境的母环境，一直到顶层环境；顶层环境一般是全局环境（工作空间）或者扩展包的命名空间。
• 顶层环境后，继续搜索直到碰到空环境。如果到达空环境，还没找到变量名对应的值，R就会报错。

如果在全局环境中定义函数，自由变量的值可以在用户的工作空间中找到，这是大多数情况下的正确做法。然而在R中，允许在函数中定义函数（类似C的编程语言不允许这样做）。在这种情况下，函数的定义环境就是其他函数的函数体。

---

在函数中定义另一个函数：

make.power <- function(n) {
        pow <- function(x) {
                x^n 
        }
        pow 
}

该函数返回函数体内定义的另一个函数作为返回值。

> cube <- make.power(3)
> square <- make.power(2)
> cube(3)
[1] 27
> square(3)
[1] 9

探索函数闭包

函数的环境：

> ls(environment(cube))
[1] "n"   "pow"
> get("n", environment(cube))
[1] 3

> ls(environment(square))
[1] "n"   "pow"
> get("n", environment(square))
[1] 2

静态作用域vs动态作用域

声明两个函数：

y <- 10

f <- function(x) {
        y <- 2
        y^2 + g(x)
}

g <- function(x) { 
        x*y
}

> f(3)
[1] 34

• g 函数中 y 在lexical scoping的规则下，在函数定义的环境中（在本例中为全局环境）搜索 y 对应的值为10。
• 在动态作用域规则下， y 对应的值应该在函数被调用的环境中搜索。
  • R中的调用环境称为parent frame。
  • 根据动态作用规则 y 值应当为2。

当一个函数在全局环境环境中声明并调用时，定义环境和调用环境是相同的，有时表现出动态作用的特征。

> rm(list=ls())
> g <- function(x) { 
a <- 3
x+a+y 
}
> g(2)
Error in g(2) : object "y" not found
> y <- 3
> g(2)
[1] 8

其他支持lexical scoping的编程语言

• Scheme
• Perl
• Python
• Common Lisp

Lexical Scoping的重要性

• R中的所有对象都存储在内存中，所有函数都需要一个指向该函数定义环境的指针（可以在任何地方）。
• S-PLUS一般在全局工作空间中搜索自由变量，因此所有对象可以储存在磁盘上，因为所有函数的定义环境都相同。

应用：优化

• R中的优化函数 optim, nlm, and optimize 要求传递一个参数的向量 (如log-likelihood)
• 目标函数可能依赖于除参数之外的许多东西(如 data)
• 当对程序进行优化，用户更关注特定的参数

最大化标准似然

构建函数：

make.NegLogLik <- function(data, fixed=c(FALSE,FALSE)) {
        params <- fixed
        function(p) {
                params[!fixed] <- p
                mu <- params[1]
                sigma <- params[2]
                a <- -0.5*length(data)*log(2*pi*sigma^2)
                b <- -0.5*sum((data-mu)^2) / (sigma^2)
                -(a + b)
        } 
}

注意：R中的优化函数都是最小化函数，需要使用负对数似然求最大化

> set.seed(1); normals <- rnorm(100, 1, 2)
> nLL <- make.NegLogLik(normals)
> nLL
function(p) {
                params[!fixed] <- p
                mu <- params[1]
                sigma <- params[2]
                a <- -0.5*length(data)*log(2*pi*sigma^2)
                b <- -0.5*sum((data-mu)^2) / (sigma^2)
                -(a + b)
        }
<bytecode: 0x7feaed346eb0>
<environment: 0x7feae35d1bc8>
> ls(environment(nLL))
[1] "data"   "fixed"  "params"

参数估计

> optim(c(mu = 0, sigma = 1), nLL)$par
      mu    sigma
1.218239 1.787343

当σ = 2时

> nLL <- make.NegLogLik(normals, c(FALSE, 2))
> optimize(nLL, c(-1, 3))$minimum
[1] 1.217775

当μ = 1时

> nLL <- make.NegLogLik(normals, c(1, FALSE))
> optimize(nLL, c(1e-6, 10))$minimum
[1] 1.800596

绘制最大似然图

nLL <- make.NegLogLik(normals, c(1, FALSE))
x <- seq(1.7, 1.9, len = 100)
y <- sapply(x, nLL)
plot(x, exp(-(y - min(y))), type = "l")

nLL <- make.NegLogLik(normals, c(FALSE, 2))
x <- seq(0.5, 1.5, len = 100)
y <- sapply(x, nLL)
plot(x, exp(-(y - min(y))), type = "l")

小结

• 目标函数可以包含评估该函数的所有必须数据
• 对于交互式和探索的工作不需要提供太长的参数列表
• 代码可以更简洁

扩展阅读： Robert Gentleman and Ross Ihaka (2000). “Lexical Scope and Statistical Computing,” JCGS, 9, 491–508.

编写R代码

1. 使用文本文件或编辑器 Always use text files / text editor
2. 适当的缩进
3. 限制代码的宽度
4. 限制自定义函数的长度

---

缩进

• 缩进提高代码的可读性
• 限制代码的长度，避免过多的嵌套和太长的函数
• 最少缩进4个空格，8个更完美

日期与时间

R有一套特殊日期和时间的表示方法：

• 日期使用 Date 类来表示
• 时间是通过 POSIXct 或 POSIXlt 类来表示
• R中储存的日期从1970-01-01这天开始的
• R中储存的时间是从1970-01-01的第一秒开始的

日期

日期使用Date类来表示，可以使用 as.Date() 函数从字符串转换成日期。

> x <- as.Date("1970-01-01")
> x
[1] "1970-01-01"
> unclass(x)
[1] 0
> unclass(as.Date("1970-01-02"))
[1] 1

时间

时间使用 POSIXct 或 POSIXlt 类来表示。

• POSIXct 实际是一个十分巨大的整数，当使用数据框来储存时间时十分方便。
• POSIXlt 是一个列表，存储了年月周的哪一天等额外信息

有很多基本函数可以操作日期和时间

• weekdays: 返回一周的第几天
• months: 返回月份名称
• quarters: 返回季度信息(“Q1”, “Q2”, “Q3”, or “Q4”)

通过 as.POSIXlt 或 as.POSIXct 函数将字符串转换为日期。

> x <- Sys.time()
> x
[1] "2019-05-08 21:24:57 CST"
> p <- as.POSIXlt(x)
> names(unclass(p))
 [1] "sec"    "min"    "hour"   "mday"   "mon"    "year"   "wday"   "yday"  
 [9] "isdst"  "zone"   "gmtoff"
> p$sec
[1] 57.742

使用 POSIXct 格式表示日期。

# Sys.time()生成的时间默认为 ‘POSIXct’ 格式
> x 
[1] "2019-05-08 21:24:57 CST"
> unclass(x)
[1] 1557321898
> x$sec
Error in x$sec : $ operator is invalid for atomic vectors
> p <- as.POSIXlt(x)
> p$sec
[1] 57.742

使用 strptime 函数改变日期的格式。

> datestring <- c("May 8, 2019 10:30", "May 9, 2019 09:10")
> x <- strptime(datestring, "%B %d, %Y %H:%M")
> x
[1] "2019-05-08 10:30:00 CST" "2019-05-09 09:10:00 CST"
> class(x)
[1] "POSIXlt" "POSIXt" 
# 查看更多时间格式的帮助文档
> ?strptime

日期和时间的操作

日期和时间也可以用于数学计算（如+, -），也可以做比较 (如 ==, <=)。

> x <- as.Date("2020-05-08")
> y <- strptime("8 May 2019 10:34:21", "%d %b %Y %H:%M:%S") 
> x-y
Error in x - y : non-numeric argument to binary operator
In addition: Warning message:
Incompatible methods ("-.Date", "-.POSIXt") for "-" 
> x <- as.POSIXlt(x) 
> x-y
Time difference of 365.8928 days

也可以计算甚至可以跟踪闰年、闰秒、夏令时和时区。

> x <- as.Date("2019-03-01") 
> y <- as.Date("2019-02-28") 
> x-y
Time difference of 1 days
> x <- as.POSIXct("2019-05-08 01:00:00")
> y <- as.POSIXct("2019-05-08 06:00:00", tz = "GMT") 
> y-x
Time difference of 13 hours

小结

• 日期和时间是R中特殊的类，允许数值和统计计算
• 日期使用 Date 类
• 时间使用 POSIXct 和 POSIXlt 类
• 可以使用 strptime, as.Date, as.POSIXlt,  as.POSIXct 函数将字符串转换为日期/时间类

apply函数家族

命令行的循环

循环对于编程十分重要，但在R交互式命令行中编写循环结构不太方便，R提供了一系列函数来简化循环操作。

• lapply: 对一个列表进行循环，对列表中的，每个元素执行指定的函数操作，返回与原列表等长的结果
• sapply: 与 lapply 相同，简化结果
• apply: 对数组或矩阵执行指定的函数，返回向量、数组或列表
• tapply: 对向量的子集进行操作
• mapply: lapply 和 sapply 的升级版

split 函数辅助 lapply 会产生一些神操作。

---

lapply

lapply 需要三个参数： (1) 列表 X; (2) 函数或函数名FUN; (3) 其他参数通过 ... 参数传递。 如果 X 不是列表，会通过 as.list 函数强制转换为列表。

> lapply
function (X, FUN, ...) 
{
    FUN <- match.fun(FUN)
    if (!is.vector(X) || is.object(X)) 
        X <- as.list(X)
    .Internal(lapply(X, FUN))
}
<bytecode: 0x7fd35e9dc718>
<environment: namespace:base>

实际循环操作是通过内部的C代码完成。

---

lapply 函数或返回一个列表，但会忽略输入数据的类型。

> x <- list(a = 1:5, b = rnorm(10))
> lapply(x, mean)
$a
[1] 3

$b
[1] -0.1843325

> x <- list(a = 1:4, b = rnorm(10), c = rnorm(20, 1), d = rnorm(100, 5))
> lapply(x, mean)
$a
[1] 2.5

$b
[1] -0.1599965

$c
[1] 1.068789

$d
[1] 4.837933

> x <- 1:4
> lapply(x, runif)
[[1]]
[1] 0.9785864

[[2]]
[1] 0.5644791 0.9242951

[[3]]
[1] 0.9274006 0.3646827 0.2551950

[[4]]
[1] 0.3987640 0.5774553 0.7631348 0.7837444

> x <- 1:4
> lapply(x, runif, min = 0, max = 10)
[[1]]
[1] 4.615299

[[2]]
[1] 3.011078 5.549246

[[3]]
[1] 8.5230253 3.4520337 0.1495026

[[4]]
[1] 5.9855820 7.7427082 2.9215029 0.3520867

lapply 也可以使用匿名函数

# 声明包含两个矩阵的列表
> x <- list(a = matrix(1:4, 2, 2), b = matrix(1:6, 3, 2)) 
> x
$a
     [,1] [,2]
[1,]    1    3
[2,]    2    4

$b
     [,1] [,2]
[1,]    1    4
[2,]    2    5
[3,]    3    6


# 提取每个矩阵的第一列的匿名函数
> lapply(x, function(elt) elt[,1])
$a
[1] 1 2

$b
[1] 1 2 3

sapply

sapply 会尽可能的简化 lapply 函数的结果。

• 如果列表中的每个元素的长度为1，结果会返回一个向量。
• 如果列表中的每个元素都是长度大于1的等长向量，结果会返回一个矩阵。
• 如果无法简化，返回列表。

> x <- list(a = 1:4, b = rnorm(10), c = rnorm(20, 1), d = rnorm(100, 5))
> lapply(x, mean)
$a
[1] 2.5

$b
[1] 0.3727818

$c
[1] 1.17464

$d
[1] 5.07108

> sapply(x, mean)
        a         b         c         d 
2.5000000 0.3727818 1.1746396 5.0710795 
> mean(x)
[1] NA
Warning message:
In mean.default(x) : argument is not numeric or logical: returning NA

apply

apply 对数组或矩阵执行指定的函数（通常是匿名函数）。

• 通常对矩阵的一行或一列操作
• 对一组数进行操作，如计算矩阵的行平均值或列平均值。
• 不一定比循环快，但只要一行。

> str(apply)
function (X, MARGIN, FUN, ...)

• X 是一个数组
• MARGIN 是整数型向量，指定对哪一维（如矩阵的行或列）进行操作。
• FUN 是用哪个函数进行操作
• ... 是传递给 FUN 的其他参数

> x <- matrix(rnorm(200), 20, 10)
> apply(x, 2, mean)
 [1]  0.167162527 -0.079293974 -0.062300596 -0.328406829  0.290078933
 [6]  0.480642185  0.009369719 -0.018753792  0.194263160 -0.042819273
> apply(x, 1, sum)
 [1] -0.6314381  1.3082838 -0.9322577  9.2538844 -6.0182467  4.2462860
 [7]  1.3619095 -1.6376867 -2.0452763 -3.1957812 -2.6102388  3.8403223
[13] -1.0428887  3.0235110  1.1093232  3.0340000  0.5430936  1.1590106
[19]  0.4819350  0.9510959

列/行求和与平均值

对于计算矩阵维度的和与平均值，有更快捷的函数。

• rowSums = apply(x, 1, sum)
• rowMeans = apply(x, 1, mean)
• colSums = apply(x, 2, sum)
• colMeans = apply(x, 2, mean)

这些函数速度很快，在不处理大矩阵的时，没有区别。

---

其他使用apply的方法

计算一个矩阵行的分位数。

> x <- matrix(rnorm(200), 20, 10)
> apply(x, 1, quantile, probs = c(0.25, 0.75))
          [,1]      [,2]       [,3]      [,4]       [,5]       [,6]       [,7]
25% 0.09678483 -1.252172 -0.5594222 0.3003727 -0.5394311 -0.7213852 -0.6946071
75% 1.21133570  1.125509  0.4887024 1.2867236  0.4701446  0.3500056  0.1762477
          [,8]       [,9]      [,10]     [,11]      [,12]      [,13]     [,14]
25% -0.2561037 -0.2304411 -0.9268922 -1.079129 -0.7953414 -0.4770530 -0.379590
75%  0.6587894  1.5569296  1.2895396  1.142885 -0.1241847  0.7192662  0.545214
           [,15]      [,16]      [,17]      [,18]      [,19]      [,20]
25% -0.008474735 -0.3288624 -0.5923715 -0.3822122 -0.6172781 -1.2056816
75%  0.998829133  1.0367678  0.8497671  0.6897986  0.3691122  0.1918069

求数组中矩阵的平均值

> a <- array(rnorm(2 * 2 * 10), c(2, 2, 10))
> apply(a, c(1, 2), mean)
            [,1]       [,2]
[1,] 0.006563965 -0.4093383
[2,] 0.057019258  0.1777483
> rowMeans(a, dims = 2)
            [,1]       [,2]
[1,] 0.006563965 -0.4093383
[2,] 0.057019258  0.1777483

mapply

mapply 称为多变量（multivariate）apply可以接收指定函数的一组参数并行计算。

> str(mapply)
function (FUN, ..., MoreArgs = NULL, SIMPLIFY = TRUE, USE.NAMES = TRUE)

• FUN 是用来批处理的函数
• ... 用来批处理的参数
• MoreArgs 是 FUN函数的其他参数
• SIMPLIFY 设置是否简化参数

例如要生成包含四列的列表，分别包含4个1，3个2，2个3，1个4四个向量。 list(rep(1, 4), rep(2, 3), rep(3, 2), rep(4, 1))

可以使用mapply来简化操作

> mapply(rep, 1:4, 4:1)
[[1]]
[1] 1 1 1 1

[[2]]
[1] 2 2 2

[[3]] 
[1] 3 3

[[4]] 
[1] 4

函数的向量化操作

# 声明一个函数noise
> noise <- function(n, mean, sd) {
 rnorm(n, mean, sd)
 }
# 使用单个参数操作
> noise(5, 1, 2)
[1]  0.1629689  2.0816275  2.4775372 -0.9368394  0.9255463
# 使用多个参数操作，结果不是预期的
> noise(1:5, 1:5, 2)
[1] 1.643462 2.794333 2.914787 5.082487 5.377841

向量化操作

> mapply(noise, 1:5, 1:5, 2)
[[1]]
[1] 2.591011

[[2]]
[1] 1.851819 1.861633

[[3]]
[1] 2.7667496 0.4649369 5.8556238

[[4]]
[1] 5.275667 2.682882 4.043813 4.366316

[[5]]
[1] 6.257788 3.618693 4.407559 6.736056 6.720434

等同于

list(noise(1, 1, 2), noise(2, 2, 2),
     noise(3, 3, 2), noise(4, 4, 2),
     noise(5, 5, 2))

tapply

tapply 对向量的子集执行批处理操作。

> str(tapply)
function (X, INDEX, FUN = NULL, ..., simplify = TRUE)

• X 是一个向量
• INDEX 因子或因子的列表（或与因子相关）
• FUN 批处理的函数
• ... 其他传递给 FUN 函数
• simplify 是否简化结果

分组取平均值。

> x <- c(rnorm(10), runif(10), rnorm(10, 1))
> f <- gl(3, 10)
> f
 [1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3
Levels: 1 2 3
> tapply(x, f, mean)
        1         2         3 
0.1045255 0.4867243 0.9131191

不简化分组取平均值的结果

> tapply(x, f, mean, simplify = FALSE)
$`1`
[1] 0.1045255

$`2`
[1] 0.4867243

$`3`
[1] 0.9131191

找到每组的数值范围。

> tapply(x, f, range)
$`1`
[1] -0.8040998  1.0022698

$`2`
[1] 0.04577595 0.95238798

$`3`
[1] -0.4422177  2.3863979

split

split 根据因子向量或因子列表分组。

> str(split)
function (x, f, drop = FALSE, ...)

• x 可以是向量、列表、数据框
• f 因子或因子列表
• drop 是否去除因子水平为空的结果

> split(x, f)
$`1`
 [1]  0.06417511  0.77601085  1.66855356  1.38744423
 [5] -0.90908770  0.39727163 -2.13528805  0.29087121
 [9]  0.82936584  0.53773723

$`2`
 [1] 0.6646064 0.4408925 0.3199122 0.2156969 0.8358507
 [6] 0.1408568 0.4088236 0.2258691 0.9606134 0.7945027

$`3`
 [1]  0.65276220  2.46645556  2.72756544  1.77246304
 [5]  2.94941952  0.11977102 -0.04283368  2.36610370
 [9]  0.44573942  2.31295594

lapply 和 split 配合使用的例子。

> lapply(split(x, f), mean)
$`1`
[1] 0.2907054

$`2`
[1] 0.5007624

$`3`
[1] 1.57704

拆分数据框

> library(datasets)
> head(airquality)
  Ozone Solar.R Wind Temp Month Day
1    41     190  7.4   67     5   1
2    36     118  8.0   72     5   2
3    12     149 12.6   74     5   3
4    18     313 11.5   62     5   4
5    NA      NA 14.3   56     5   5
6    28      NA 14.9   66     5   6

> s <- split(airquality, airquality$Month)
> lapply(s, function(x) colMeans(x[, c("Ozone", "Solar.R", "Wind")]))
$`5`
   Ozone  Solar.R     Wind 
      NA       NA 11.62258 

$`6`
    Ozone   Solar.R      Wind 
       NA 190.16667  10.26667 

$`7`
     Ozone    Solar.R       Wind 
        NA 216.483871   8.941935 

$`8`
   Ozone  Solar.R     Wind 
      NA       NA 8.793548 

$`9`
   Ozone  Solar.R     Wind 
      NA 167.4333  10.1800 

> sapply(s, function(x) colMeans(x[, c("Ozone", "Solar.R", "Wind")])) 
               5         6          7        8        9
Ozone         NA        NA         NA       NA       NA
Solar.R       NA 190.16667 216.483871       NA 167.4333
Wind    11.62258  10.26667   8.941935 8.793548  10.1800
> sapply(s, function(x) colMeans(x[, c("Ozone", "Solar.R", "Wind")], na.rm = TRUE))
                5         6          7          8         9
Ozone    23.61538  29.44444  59.115385  59.961538  31.44828
Solar.R 181.29630 190.16667 216.483871 171.857143 167.43333
Wind     11.62258  10.26667   8.941935   8.793548  10.18000

将数据分割到多个水平

> x <- rnorm(10)
> f1 <- gl(2, 5)
> f2 <- gl(5, 2)
> f1
 [1] 1 1 1 1 1 2 2 2 2 2
Levels: 1 2
> f2
 [1] 1 1 2 2 3 3 4 4 5 5
Levels: 1 2 3 4 5
> interaction(f1, f2)
 [1] 1.1 1.1 1.2 1.2 1.3 2.3 2.4 2.4 2.5 2.5
Levels: 1.1 2.1 1.2 2.2 1.3 2.3 1.4 2.4 1.5 2.5

interaction 函数可以创建没有值的水平。

> str(split(x, list(f1, f2)))
List of 10
 $ 1.1: num [1:2] -1.073 -0.572
 $ 2.1: num(0) 
 $ 1.2: num [1:2] -0.435 -0.976
 $ 2.2: num(0) 
 $ 1.3: num 0.201
 $ 2.3: num 0.691
 $ 1.4: num(0) 
 $ 2.4: num [1:2] 0.498 0.333
 $ 1.5: num(0) 
 $ 2.5: num [1:2] -0.00797 -0.1332

去除无值的分组。

> str(split(x, list(f1, f2), drop = TRUE))
List of 6
 $ 1.1: num [1:2] -1.073 -0.572
 $ 1.2: num [1:2] -0.435 -0.976
 $ 1.3: num 0.201
 $ 2.3: num 0.691
 $ 2.4: num [1:2] 0.498 0.333
 $ 2.5: num [1:2] -0.00797 -0.1332

异常处理

有时R会返回一些信息

• message: message 函数生成的一般提示或诊断信息
• warning: 有些代码书写不恰当，但不会对结果造成太大影响（可忽略）；函数执行期间的提示信息；通过 warning 函数生成警告信息
• error: 程序运行出现错误，运行终止，通过 stop 函数输出
• condition: 选项，用于不能预先知道某些结果时；程序员可以创建自定义条件

Warning

> log(-1)
[1] NaN
Warning message:
In log(-1) : NaNs produced

> printmessage <- function(x) {
        if(x > 0)
                print("x is greater than zero")
        else
                print("x is less than or equal to zero")
        invisible(x)
}

> printmessage(1)
[1] "x is greater than zero"

> printmessage(NA)
Error in if (x > 0) print("x is greater than zero") else print("x is less than or equal to zero") : 
  missing value where TRUE/FALSE needed

> printmessage2 <- function(x) {
        if(is.na(x))
                print("x is a missing value!")
        else if(x > 0)
                print("x is greater than zero")
        else
                print("x is less than or equal to zero")
        invisible(x)
}

> x <- log(-1)
Warning message:
In log(-1) : NaNs produced

> printmessage2(x)
[1] "x is a missing value!"

处理函数的报错

• 注意输入与函数的调用
• 哪些输出，信息或结果是正确的
• 注意当前的输出
• 当前输出与预期结果的差异
• 一开始预想的预期结果是否正确
• 什么情况下会再次出现该问题

R中的Debug调试工具

R自带的调试函数：

• traceback: 当error产生时，调用该函数，输出调用日志
• debug: 将一个函数标记为 “debug” 模式，允许一行一行运行函数就行调试
• browser: 挂起一个函数的运行（无论该函数是否被调用），将该函数加入调试模式
• trace: 将调试代码插入函数的指定位置
• recover: 调整报错的信息，更快地定位函数

这些函数是专门为交互式编程调试设计，将这些调试工具封装为函数。含有一些更为直接的方法，通过调用print/cat函数输出提示信息。

---

traceback

> mean(x)
Error in mean(x) : object 'x' not found
> traceback()
1: mean(x)

> lm(y ~ x)
Error in eval(predvars, data, env) : object 'y' not found
> traceback()
7: eval(predvars, data, env)
6: eval(predvars, data, env)
5: model.frame.default(formula = y ~ x, drop.unused.levels = TRUE)
4: stats::model.frame(formula = y ~ x, drop.unused.levels = TRUE)
3: eval(mf, parent.frame())
2: eval(mf, parent.frame())
1: lm(y ~ x)

调试

> debug(lm)
> lm(y ~ x)
debugging in: lm(y ~ x)
debug: {
    ret.x <- x
    ret.y <- y
    cl <- match.call()
    mf <- match.call(expand.dots = FALSE)
    m <- match(c("formula", "data", "subset", "weights", "na.action", 
        "offset"), names(mf), 0L)
    mf <- mf[c(1L, m)]
    mf$drop.unused.levels <- TRUE
    mf[[1L]] <- quote(stats::model.frame)
    mf <- eval(mf, parent.frame())
    if (method == "model.frame") 
        return(mf)
    else if (method != "qr") 
        warning(gettextf("method = '%s' is not supported. Using 'qr'", 
            method), domain = NA)
    mt <- attr(mf, "terms")
    y <- model.response(mf, "numeric")
    w <- as.vector(model.weights(mf))
    if (!is.null(w) && !is.numeric(w)) 
        stop("'weights' must be a numeric vector")
    offset <- as.vector(model.offset(mf))
    if (!is.null(offset)) {
        if (length(offset) != NROW(y)) 
            stop(gettextf("number of offsets is %d, should equal %d (number of observations)", 
                length(offset), NROW(y)), domain = NA)
    }
    if (is.empty.model(mt)) {
        x <- NULL
        z <- list(coefficients = if (is.matrix(y)) matrix(NA_real_, 
            0, ncol(y)) else numeric(), residuals = y, fitted.values = 0 * 
            y, weights = w, rank = 0L, df.residual = if (!is.null(w)) sum(w != 
            0) else if (is.matrix(y)) nrow(y) else length(y))
        if (!is.null(offset)) {
            z$fitted.values <- offset
            z$residuals <- y - offset
        }
    }
    else {
        x <- model.matrix(mt, mf, contrasts)
        z <- if (is.null(w)) 
            lm.fit(x, y, offset = offset, singular.ok = singular.ok, 
                ...)
        else lm.wfit(x, y, w, offset = offset, singular.ok = singular.ok, 
            ...)
    }
    class(z) <- c(if (is.matrix(y)) "mlm", "lm")
    z$na.action <- attr(mf, "na.action")
    z$offset <- offset
    z$contrasts <- attr(x, "contrasts")
    z$xlevels <- .getXlevels(mt, mf)
    z$call <- cl
    z$terms <- mt
    if (model) 
        z$model <- mf
    if (ret.x) 
        z$x <- x
    if (ret.y) 
        z$y <- y
    if (!qr) 
        z$qr <- NULL
    z
}
Browse[2]> 
Browse[2]> n
debug: ret.x <- x
Browse[2]> ret.x <- x
Browse[2]> n
debug: ret.y <- y
Browse[2]> ret.y <- y
Browse[2]> n
debug: cl <- match.call()
Browse[2]> cl <- match.call()
Browse[2]> n
debug: mf <- match.call(expand.dots = FALSE)
Browse[2]> mf <- match.call(expand.dots = FALSE)
Browse[2]> n
debug: m <- match(c("formula", "data", "subset", "weights", "na.action", 
    "offset"), names(mf), 0L)
# 退出调试模式
Browse[2]> Q

recover

> options(error = recover)
> read.csv("nosuchfile")
Error in file(file, "rt") : cannot open the connection
In addition: Warning message:
In file(file, "rt") :
  cannot open file 'nosuchfile': No such file or directory

Enter a frame number, or 0 to exit   

1: read.csv("nosuchfile")
2: read.table(file = file, header = header, sep = sep, quote = quote, dec = dec, f
3: file(file, "rt")

Selection: 
Selection: 0
>

小结

• R有三个主要的显示错误、提示信息的函数，分别为：message, warning, error
  • 程序员的世界只有 error ，只有 error 需要调试
• 在分析报错的函数时，首先确认问题的可重复性，明确期望结果、输出与期望结果的差异
• 交互式调试工具 traceback, debug, browser, trace 和 recover 可以帮助判断函数中错误代码
• 调试工具不能取代思考

生成随机数

R中的概率分布函数

• rnorm: 根据给定的平均值和标准偏差生成随机正态变量
• dnorm: 计算一个点或点向量的正态概率密度（根据给定的平均值/ 标准偏差）
• pnorm: 生成正态累积分布函数
• rpois: 根据给定的比率生成随机泊松分布随机数

每种概率分布函数通常有四种前缀

• d 密度
• r 生成随机数
• p 累积分布
• q 分位数函数

生成正态分布可以使用以下四个函数

dnorm(x, mean = 0, sd = 1, log = FALSE)
pnorm(q, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
qnorm(p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
rnorm(n, mean = 0, sd = 1)

如果$\Phi$是标准正态累积分布i， 那么pnorm(q) = $\Phi (q)$ ， qnorm(p) = $\phi ^{-1}(p)$。

> x <- rnorm(10) 
> x
 [1]  0.3357161 -0.4786192  0.3952794 -1.5230122 -0.6496318 -1.2714863  0.6367861
 [8] -0.8809022 -0.4377379 -0.3063769
> x <- rnorm(10, 20, 2) 
> x
 [1] 16.15086 15.89892 23.22139 20.60856 25.15596 18.85948 19.50671 20.81849
 [9] 17.82214 18.43590
> summary(x)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  15.90   17.98   19.18   19.65   20.77   25.16

使用 set.seed 函数设定随机数种子以保证结果的可重复性

> set.seed(1)
> rnorm(5)
[1] -0.6264538  0.1836433 -0.8356286  1.5952808  0.3295078
> rnorm(5)
[1] -0.8204684  0.4874291  0.7383247  0.5757814 -0.3053884
> set.seed(1)
> rnorm(5)
[1] -0.6264538  0.1836433 -0.8356286  1.5952808  0.3295078

在进行模拟时最好设置随机数种子

生成泊松分布数据

> rpois(10, 1)
 [1] 0 0 1 1 2 1 1 4 1 2
> rpois(10, 2)
 [1] 4 1 2 0 1 1 0 1 4 1
> rpois(10, 20)
 [1] 19 19 24 23 22 24 23 20 11 22
## 累积分布‘
## x <= 2的概率分布
> ppois(2, 2)  
[1] 0.6766764
## x <= 4的概率分布 
> ppois(4, 2)
[1] 0.947347
## x <= 6的概率分布
> ppois(6, 2) 
[1] 0.9954662

通过线性模型生成随机数

通过$y=\beta _{0}+\beta _{1}x+\varepsilon$模型生成随机数。 已知$\varepsilon \sim {\mathcal {N}}(0,2^{2})$​. Assume $x\sim {\mathcal {N}}(0,1^{2})$​, ​$\beta _{0}=0.5$ ， ​$\beta _{1}=2$。

线性模型生成随机数1

> set.seed(20)
> x <- rnorm(100)
> e <- rnorm(100, 0, 2)
> y <- 0.5 + 2 * x + e
> summary(y)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-6.4084 -1.5402  0.6789  0.6893  2.9303  6.5052 
> plot(x, y)

当 x 为0或1时

线性模型生成随机数2

set.seed(10)
x <- rbinom(100, 1, 0.5)
e <- rnorm(100, 0, 2)
y <- 0.5 + 2 * x + e
plot(x, y)

通过广义线性模型生成随机数

模拟泊松模型：Y ~ Poisson(μ)

log μ = ​$\beta _{0}+\beta _{1}x$

已知$\beta _{0}=0.5$和$\beta _{1}=0.3$​。

使用 rpois 函数生成随机数

set.seed(1)
x <- rnorm(100)
log.mu <- 0.5 + 0.3 * x
y <- rpois(100, exp(log.mu))
summary(y)
plot(x, y)

生成随机样本

sample 函数可以从一个指定的组合中生成随机样本。

> set.seed(1)
> sample(1:10, 4)
[1] 3 4 5 7
> sample(1:10, 4)
[1] 3 9 8 5
> sample(letters, 5)
[1] "q" "b" "e" "x" "p"
## 随机排列组合
> sample(1:10)  
 [1]  4  7 10  6  9  2  8  3  1  5
> sample(1:10)
 [1]  2  3  4  1  9  5 10  8  6  7
## 随机替换，生成可重复的样本
> sample(1:10, replace = TRUE)  ## Sample w/replacement
 [1] 2 9 7 8 2 8 5 9 7 8

小结

• r* 系列函数根据特殊的概率分布生成随机数
• 内置的正态分布：标准，泊松，二项式，指数, 伽马等
• sample 函数提取随机样本
• 通过set.seed函数生成随机数生成种子，保证结果的可重复性

程序设计

• 检查程序代码不同部分的执行时间对于代码优化很有用。
• 通常代码运行一次会表现良好，但如果你将其放如1000次循环中，运行效率是否受到影响，这时分析运行时间很有用。

优化代码

• 影响代码执行速度的最大因素是代码花费最多时间的部分
• 不通过程序执行时间分析很难完成

We should forget about small efficiencies, say about 97% of the time: premature optimization is the root of all evil（过早最优化是万恶的根源）--Donald Knuth

---

优化的一般目标

• 先设计后优化
• 过早优化是万恶的根源
• 操作（收集数据），不要猜测
• 科学家的基本素养

使用 system.time()

• 接收任意R表达是作为输入，返回运行该表达式所需要的时间
• 计算执行表达所需要的时间（秒）
  • 如果有错误，返回报错前的时间
• 返回 proc_time 类的对象
  • user time: CPU时间
  • elapsed time: 总流逝时间
• 对于直接的计算任务，user time和elapsed time很接近
• 如果CPU花费大量的等待时间，elapsed time会比user time大很多
• 如果机器拥有多线程elapsted time会比user time小很多
  • 多线程BLAS库 (vecLib/Accelerate, ATLAS, ACML, MKL)
  • 多线程包parallel

> system.time(readLines("http://www.baidu.com"))
   user  system elapsed 
  0.026   0.010   1.231 

> hilbert <- function(n) { 
        i <- 1:n
        1 / outer(i - 1, i, "+")
}
> x <- hilbert(1000)
> system.time(svd(x))
   user  system elapsed 
  2.722   0.029   2.792

计算更长表达式的运行时间

> system.time({
    n <- 1000
    r <- numeric(n)
    for (i in 1:n) {
        x <- rnorm(n)
        r[i] <- mean(x)
    }
})
   user  system elapsed 
  0.079   0.005   0.086

system.time()的局限

system.time() 允许测试特定的函数或代码块的运行时间，这是在知道问题在哪的时候，但如果不知道问题在哪呢？

---

R 分析器

Rprof() 函数可以启动分析，summaryRprof() 函数可以总结 Rprof() 的输出结果，不要同时使用 system.time() 和 Rprof() 函数

• Rprof() 有规律地追踪各个函数的调用，采样收集每个函数所花费的时间
• 默认采样间隔为0.02秒
• 如果代码已经运行很快，你就不需要使用这些分析函数

R分析的原始输出

## lm(y ~ x)

sample.interval=10000
"list" "eval" "eval" "model.frame.default" "model.frame" "eval" "eval" "lm" 
"list" "eval" "eval" "model.frame.default" "model.frame" "eval" "eval" "lm" 
"list" "eval" "eval" "model.frame.default" "model.frame" "eval" "eval" "lm" 
"list" "eval" "eval" "model.frame.default" "model.frame" "eval" "eval" "lm" 
"na.omit" "model.frame.default" "model.frame" "eval" "eval" "lm" 
"na.omit" "model.frame.default" "model.frame" "eval" "eval" "lm" 
"na.omit" "model.frame.default" "model.frame" "eval" "eval" "lm" 
"na.omit" "model.frame.default" "model.frame" "eval" "eval" "lm" 
"na.omit" "model.frame.default" "model.frame" "eval" "eval" "lm" 
"na.omit" "model.frame.default" "model.frame" "eval" "eval" "lm" 
"na.omit" "model.frame.default" "model.frame" "eval" "eval" "lm" 
"lm.fit" "lm" 
"lm.fit" "lm" 
"lm.fit" "lm"

使用 summaryRprof()

• summaryRprof() 计算每个函数的运行时间
• 两种汇总方法
  • "by.total" 使用每个函数运行时间除以总运行时间
  • "by.self" 减去调用的时间

小结

• Rprof() 评估R代码的运行
• summaryRprof() 汇总Rprof()的结果，给出每个函数运行的时间百分比（两种汇总方法）