Coursera R Programming Week3 心得筆記
R 語言程序開發
約翰霍普金斯大學 公共衛生學院
R Programming
Johns Hopkins Bloomberg School of Public Health
Week 3
第三週學習筆記
2015年6月1日 - 2015年6月29日
R Programming: Week 3
We have now entered the third week of R Programming which also marks the halfway point. The lectures this week cover loop functions and the debugging tools in R. These aspects of R make R useful for both interactive work and writing longer code, and so they are commonly used in practice.
The Programming Assignment is challenging and so I encourage you to start early if you have the chance. It requires you to explore some of the more interesting aspects of the R language, including taking advantage of the scoping rules to implement state preservation in R objects.
Note that the programming assignment this week is implemented as a Peer Assessment so you will not see it listed with the other Programming Assignments. Please go to the Programming Assignment 2 section of the course to find the assignment instructions. Also, for this assignment, you will need to setup your GitHub account if you have not yet done so.
Best of luck!
Roger Peng and the Data Science Team
Mon 15 Jun 2015 8:01 AM CST
Week 3: Loop Functions and Debugging
This week is what I call "loop functions" in R, which are functions that allow you to execute loop-like behavior in a compact form. These functions typically have the word "apply" in them and are particularly convenient when you need to execute a loop on the command line when using R interactively. These functions are some of the more interesting functions of the R language. This week we also cover the debugger that comes with R and how to interpret its output to help you find problems in your programs and functions.
Learning Objectives
By the end of this week you should be able to:
- Define an anonymous function and describe its use in loop functions [see lapply]
- Describe how to start the R debugger for an arbitrary R function
- Describe what the traceback() function does and what is the function call stack
Programming
There is a graded programming assignment for this week. Please note that this assignment is graded via peer assessment.
- Programming assignment 2: Lexical Scoping
Week 3
Loop Functions - lapply
當你想要對一個或一組對象執行循環時，loop可以使你在非常少的空間內執行大量的工作，這樣你就不必使用命令行做很多輸入。當然for、while都是很有用的函數，都某種意義上，它們並不夠簡潔。在R語言裡有幾個循環函數，它們的名字裡通常帶有apply一詞：
lapply(): 主要用途是有一個對象列表，而你想要列出這個列表，並對列表裡的每一個元素運用函數。只需要用很少量的輸入，就可以完成非常強大的計算。
sapply(): lapply()的一個變化，它簡化了lapply()的結果。
apply(): 一個對數組進行行或列運算的函數，如果你想對矩陣或其他高維數組求和，這函數相當好用。
tapply(): 是table apply()的縮寫，它將函數應用於向量的子集。
mapply(): 是lapply()的多變量版本。
split(): 它不對對象進行任何操作，但它和lapply()或sapply()組合使用時效果出色。因為它會把對象分成子塊。
> lapply
function (X, FUN, ...)
{
FUN <- match.fun(FUN)
if (!is.vector(X) || is.object(X))
X <- as.list(X)
.Internal(lapply(X, FUN))
}
<bytecode: 0x0000000005875000>
<environment: namespace:base>
參數...用來給函數，就是對列表裡的每個元素做運算的那個函數傳遞參數。如果x不是列表，可能的話，它會被as.list()強制轉化為列表，如果不能強制轉換，就會報錯。
剩餘的函數在內部適用C代碼來實現。
> x <- list(a = 1:5, b = rnorm(10))
> lapply(x, mean)
$a
[1] 3
$b
[1] -0.3767029
你要對列表中的每個元素運用函數，函數會給你返回值，它可能和列表中的原始數據不一樣。所以，它可能拿一個向量做輸入，但返回一個標量的結果，函數會針對列表中的每個對象返回一個值，返回值會組成一個新的列表，這就是lapply()的返回值。
任何可能是或不是一個列表的對象，都會被強制轉化成一個列表。
> x <- list(a = 1:5, b = rnorm(10), c = rnorm(20, 1), d = rnorm(100, 5))
> lapply(x, mean)
$a
[1] 3
$b
[1] 0.5222887
$c
[1] 1.053972
$d
[1] 4.934747
> x <- 1:4
> lapply(x, runif) # runif 會隨機生成符合均勻分布的隨機變量
[[1]] ## 由1個均勻隨機變量組成的向量
[1] 0.3414567
[[2]] ## 由2個均勻隨機變量組成的向量
[1] 0.3797821 0.3360901
[[3]] ## 由3個均勻隨機變量組成的向量
[1] 0.1041617 0.7410379 0.7475099
[[4]] ## 由4個均勻隨機變量組成的向量
[1] 0.8715441 0.3762871 0.7795657 0.7038220
> x <- 1:4
> lapply(x, runif, min = 0, max = 10)
[[1]]
[1] 0.4910904
[[2]]
[1] 9.531593 9.487243
[[3]]
[1] 2.250094 8.227296 4.207857
[[4]]
[1] 3.090255 3.825546 1.223638 7.697467
> 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
例如，結果是一個每個元素長度都為1的列表，那sapply()將會返回一個包含所有元素的向量，但通常你不會想要一組每個元素都是單個數字的列表，而sapply()就會把它簡化成一個向量。
如果結果是一個列表，裡面每個元素都是等長的向量，例如，返回列表裡都是長度為5的元素，那sapply()會把這些元素放進矩陣，這列表就是長度為5的矩陣，通常這才是你所期望的結果。
> x <- list(a = 1:5, b = rnorm(10), c = rnorm(20, 1), d = rnorm(100, 5))
> lapply(x, mean)
$a
[1] 3
$b
[1] -0.1844007
$c
[1] 0.9634149
$d
[1] 4.892546
如果我得到的是一個向量，包含了所有這些數字，那就好多了，這就是sapply()所做的事情。
> sapply(x, mean)
a b c d
3.0000000 -0.1844007 0.9634149 4.8925462
> mean(x) # mean() 不適用於列表
[1] NA
Warning message:
In mean.default(x) : argument is not numeric or logical: returning NA
Loop Functions - apply
apply()可以把一個函數應用在一個數組的各個維度上，通常這個函數會是匿名函數。
它通常的應用對象是矩陣的行或列，因為矩陣作為一個二維數組，是R中最為常見的數組模型，但有時也會碰到三維數組等等。你可以對一般數組使用apply()，比如對一個矩陣數組取平均值。
也許你有時聽說在某種程度上，使用apply()要比for循環更好，或apply()的運行速度更快，在老舊的S語言中，這樣說是對的，但就R語言來說，這樣毫無根據。
使用apply()的主要原因是，它涉及的輸入相對更少。所以apply函數很有用，尤其是在使用命令行時。因為在命令行時，當我們與數據交互，做探索性分析時，我們會希望盡可能減少輸入。
> str(apply)
function (X, MARGIN, FUN, ...)
> x <- matrix(rnorm(200), 20, 10)
> apply(x, 2, mean) # MARGIN = 2 保留所有的行(第二個維度)，消除所有的列
[1] -0.13537403 -0.35390706 0.46310221 -0.47075599 0.16097141 0.01542872 0.14100634
[8] 0.14707695 0.50926972 0.14917494
> apply(x, 1, sum) # MARGIN = 1 保留所有的列(第一個維度)，消除所有的行
[1] -7.19576005 -1.39252343 2.16720554 -2.40863768 4.03809617 4.09477830 -0.22786207
[8] 6.28879396 -4.21924554 7.80916926 1.87933843 7.67595256 -2.18177460 -2.09966415
[15] 5.05770642 -2.17141137 1.21512590 0.69800676 -0.01445567 -6.49297474
rowMeans = apply(x, 1, mean)
colSums = apply(x, 2, sum)
colMeans = apply(x, 2, mean)
> x <- matrix(rnorm(200), 20, 10)
> apply(x, 1, quantile, probs = c(0.25, 0.75))
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
25% -0.4149475 -0.3712762 -0.3015529 -1.48281189 -0.4811816 -0.8465293 -0.1397921 -0.483606
75% 0.8552090 0.9579350 0.4593349 0.08718044 0.5852659 0.3419130 1.1020690 0.732006
[,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
25% -0.9324121 -0.6140571 -0.2374867 -0.9377337 -0.6591365 -0.4673306 -0.5004981 -0.3082152
75% -0.3874239 0.8506858 0.9538714 1.1309752 0.2804113 1.1842481 0.6574628 0.7964215
[,17] [,18] [,19] [,20]
25% -0.4594000 -0.7089974 -0.7692496 -0.8420129
75% 0.9451076 0.6821850 1.4894648 0.1462780
> a <- array(rnorm(2 * 2 * 10), c(2, 2, 10))
> apply(a, c(1, 2), mean)
[,1] [,2]
[1,] -0.15282526 -0.3374177
[2,] 0.05267209 -0.3290474
> rowMeans(a, dims = 2)
[,1] [,2]
[1,] -0.15282526 -0.3374177
[2,] 0.05267209 -0.3290474
Loop Functions - mapply
mapply()是lapply()和sapply()的多變量版本，它的想法是把一個函數並行應用到一組不同的參數上。
> str(mapply)
function (FUN, ..., MoreArgs = NULL, SIMPLIFY = TRUE, USE.NAMES = TRUE)
MoreArgs是當你需要給函數傳遞更多參數時才會用到。
SIMPLIFY參數和你們在sapply()和tapply()中看到的SIMPLIFY參數差不多。
list(rep(1, 4), rep(2, 3), rep(3, 2), rep(4, 1))
> 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 <- function(n, mean, sd) {
+ rnorm(n, mean, sd)
+ }
> noise(5, 1, 2)
[1] 0.08948614 -0.22371693 3.70067239 -2.29552476 2.45988019
> noise(1:5, 1:5, 2)
[1] 1.572638 1.262969 1.501788 1.334549 3.063200
> mapply(noise, 1:5, 1:5, 2)
[[1]]
[1] 2.601908
[[2]]
[1] -0.6114717 2.8404665
[[3]]
[1] -0.6755479 4.0655113 4.6672972
[[4]]
[1] 4.0951107 3.7826142 0.3919656 4.8230890
[[5]]
[1] 10.1129603 5.3705774 4.2623584 0.6501055 4.4268859
> mapply(noise, 1:5, 1:5, 2)
[[1]]
[1] 2.601908
[[2]]
[1] -0.6114717 2.8404665
[[3]]
[1] -0.6755479 4.0655113 4.6672972
[[4]]
[1] 4.0951107 3.7826142 0.3919656 4.8230890
[[5]]
[1] 10.1129603 5.3705774 4.2623584 0.6501055 4.4268859
Loop Functions - tapply
tapply()可以把一個函數應用在向量的子集上。
> str(tapply)
function (X, INDEX, FUN = NULL, ..., simplify = TRUE)
> 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.5222887 0.5097278 1.1478271
> tapply(x, f, mean, simplify = FALSE)
$`1`
[1] 0.5222887
$`2`
[1] 0.5097278
$`3`
[1] 1.147827
> tapply(x, f, range) # range 計算觀測值的範圍
$`1`
[1] -0.5408496 1.7110455
$`2`
[1] 0.1099391 0.9004820
$`3`
[1] -0.5621404 2.9441282
Loop Functions - split
> str(split)
function (x, f, drop = FALSE, ...)
> x <- c(rnorm(10), runif(10), rnorm(10, 1))
> f <- gl(3, 10)
> split(x, f)
$`1`
[1] 1.10133042 -1.36659085 -0.06014147 -0.53467743 -0.23713734 -0.02966124 -3.31999342
[8] -0.76798265 -0.86888217 -1.22670988
$`2`
[1] 0.92967802 0.69088547 0.04862117 0.77389730 0.20059063 0.01264115 0.98192982
[8] 0.43884912 0.58794070 0.28711452
$`3`
[1] 1.40815232 0.26045983 -0.04327073 -1.28753891 0.80223386 1.34721034 1.31807625
[8] 0.67190445 2.28873492 1.72441656
> lapply(split(x, f), mean)
$`1`
[1] -0.7310446
$`2`
[1] 0.4952148
$`3`
[1] 0.8490379
> 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
> str(split(x, list(f1, f2)))
List of 10
$ 1.1: num [1:2] 1.565 -0.627
$ 2.1: num(0)
$ 1.2: num [1:2] 0.207 0.51
$ 2.2: num(0)
$ 1.3: num -1.1
$ 2.3: num 0.0792
$ 1.4: num(0)
$ 2.4: num [1:2] 0.0788 1.2059
$ 1.5: num(0)
$ 2.5: num [1:2] 0.835 0.186
> str(split(x, list(f1, f2), drop = TRUE))
List of 6
$ 1.1: num [1:2] 1.565 -0.627
$ 1.2: num [1:2] 0.207 0.51
$ 1.3: num -1.1
$ 2.3: num 0.0792
$ 2.4: num [1:2] 0.0788 1.2059
$ 2.5: num [1:2] 0.835 0.186
Debugging Tools - Diagnosing the Problem
- message: 作為一個診斷訊息，它告知有事發生，但其實可能什麼都沒發生。因為訊息不會阻止函數的執行，只是會有一條訊息顯示到螢幕上。
- warning: 警告是一種提示，意味著發生了某些出乎意料的事情，它不一定是個問題，很多時候你會想忽略它。函數期望的是一個東西，但實際得到的確不太一樣。
> log(-1)
[1] NaN
Warning message:
In log(-1) : NaNs produced- error: 錯誤會終止函數的執行，錯誤訊息是通過stop函數產生的。
- condition: 上述三種提示其實都是條件。你可能會想觸發另一種條件，它既不是錯誤、警告，也不是訊息，你可以自己創造這種條件。
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)
}
> printmessage2(NA)
[1] "x is a missing value"
> x <- log(-1)
Warning message:
In log(-1) : NaNs produced
> printmessage2(x)
[1] "x is a missing value"
你要怎麼知道你的函數哪裡有錯誤？通常你要先回答下列問題：
- 你的輸入是什麼？
- 你把什麼放進這個函數裡？
- 你怎麼調用這個函數的？
- 你給出的參數是什麼？
- 你期望得到什麼？
- 你實際得到什麼？
- 實際得到的和你期望得到的，兩者之間有什麼區別？
- 你最初的期望是不是正確的？
- 你可以重塑問題嗎？（seed ?）
Debugging Tools - Basic Tools
- traceback: 它會印出function call stack，也就是說traceback()會告訴你一共調用了幾個函數，以及錯誤發生在哪。
- debug: 你要給它傳遞一個函數作為參考，它標記這函數，進入debug模式。debug就是每當執行到這函數時，都會暫停執行，停在這函數的第一行。
- browser: 你能在代碼的任何地方調用browser()。debug()總是從函數一開始就一行行地運行，但有時你會想讓函數先執行一段，然後在某處停下，browser()讓你能在代碼的任何地方開始，函數會一直運行到哪個節點才會停止。
- trace: 它允許你在函數中插入調試代碼，這樣做避免編輯函數本身，你在調試別人代碼的時候使用它會很方便。
- recover: 通常報錯時，你會看到一條訊息告訴你錯誤是什麼，然後回到控制台，你樣你就能鍵入命令，做任何你想做的事情，但那個函數被中斷執行並回到了控制台。你能通過error handler來改變這種默認的行為，recover()就是一個錯誤處理函數，無論函數在什麼地方發生錯誤，R會停止執行但不會直接回到控制台，它會停在函數出錯的地方，並輸出function call stack讓你瀏覽。
Debugging Tools - Using the Tools
> mean(x)
Error in mean(x) : object 'x' not found
> traceback()
1: mean(x)
> lm(y ~ x)
Error in eval(expr, envir, enclos) : object 'y' not found
> traceback()
7: eval(expr, envir, enclos)
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(expr, envir, enclos)
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(, 0,
3) 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]>
實際上，還有其他參數存在lm()，只是使用了預設值沒有被列出來，所以我不能說只有一條公式，因為還有其他參數。
Browse[2]> n
debug: ret.x <- x
Browse[2]> n
debug: ret.y <- y
Browse[2]> n
debug: cl <- match.call()
Browse[2]> n
debug: mf <- match.call(expand.dots = FALSE)
Browse[2]> n
debug: m <- match(c("formula", "data", "subset", "weights", "na.action",
"offset"), names(mf), 0L)
> 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 =
3: file(file, "rt")
Selection:
在上面這個例子，你可以選擇1~3瀏覽函數環境。
- 基本上有三種關於問題或條件的提示：message、warning和error。
- 三種之中只有error會使函數停止運行。
- 當你分析一個你覺得有問題的函數時，確保你可以重現那個錯誤，並且能清楚表述你的期望是什麼、結果是什麼以及結果與你的期望有什麼不同。
- 另外有很多交互工具可以使用，例如traceback()、debug()、browser()、trace()以及recover()，重點在於＂交互＂，這些工具有用就是因為它們是交互的，你可以在控制台上進行一些操作。
- 當然，調試工具不能取代你的思考！你應該多思考代碼寫得怎樣，而不是隨便丟給調適工具，等它幫你找出問題。