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()的結果，給出每個函數運行的時間百分比（兩種匯總方法）