Coursera R Programming Week2 心得筆記
R 語言程序開發
約翰霍普金斯大學 公共衛生學院
R Programming
Johns Hopkins Bloomberg School of Public Health
Week 2
第二週學習筆記
2015年6月1日 - 2015年6月29日
R Programming: Week 2
Today marks the beginning of Week 2 of R Programming. This week we take the gloves off and the lectures cover key topics like control structures and functions. We also introduce the first programming assignment for the course, which is due at the end of the week.
A few notes about the Programming Assignment:
Each part of the Assignment can be submitted an infinite number of times---there is no limit on the number of submissions.
For each part, we take your maximum score over all of your submissions.
There is a submission script that you will have to download to submit your assignment.
Roger Peng and the Data Science Team
Mon 8 Jun 2015 8:01 AM CST
Week 2: Programming with R
This week is all about functions and about controlling the flow of an R program. We start with control structures (like if-else, and for loops) and then move on to writing functions. Next, we discuss the lexical scoping features of the language and how they can be used in interesting ways, particularly for statistical applications.
Learning Objectives
By the end of this week you should be able to:
- Write an if-else expression
- Write a for loop, a while loop, and a repeat loop
- Define a function in R and specify its return value [see Functions Part 1 and Part 2]
- Describe how R binds a value to a symbol via the search list
- Define what lexical scoping is with respect to how the value of free variables are resolved in R
- Describe the difference between lexical scoping and dynamic scoping rules
- Convert a character string representing a date/time into an R datetime object. [see Dates and Times]
Programming
There is a graded programming assignment for this week.
- Programming assignment 1: Air Pollution
This programming assignment has multiple parts and you will submit your answers using the submit script described in the instructions.
Week 2
Control Structures - Introduction
控制結構
if/else: 測試邏輯條件
for: 執行固定次數循環（loop）
while: 當某個條件成立時，執行循環（loop）
repeat: 直接地執行無限循環（loop）
break: 從任意循環跳出
next: 忽略循環中的某次迭代
return: 退出函數
Control Structures - If-else
if(<condition>) {
## do something
} else {
## do something else
}
if(<condition>) {
## do something
} else if(<condition>) {
## do something different
} else {
## do something different
}
if(x > 3) {
y <- 10
} else {
y <- 0
}
y <- if(x > 3) {
10
} else {
0
}
if(<condition1>) {
## do something
}
if(<condition2>) {
## do something
}
Control Structures - For loops
for(i in 1:10) {
print(i)
}
x <- c("a", "b", "c", "d")
for(i in 1:10) {
print(i)
}
for(i in seq_aling(x)) {
print(x(i))
}
for(letter in x) {
print(letter)
}
for(i in 1:4) print(x(i))> x <- matrix(1:6, 2, 3)
> x
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
> seq_len(ncol(x))
[1] 1 2 3
> seq_len(nrow(x))
[1] 1 2
> for(i in seq_len(nrow(x))) {
for(j in seq_len(ncol(x))) {
print(x[i, j])
}
}
for(i in 1:4) print(x(i))
[1] 1
[1] 3
[1] 5
[1] 2
[1] 4
[1] 6Control Structures - While loops
count <- 0
while(count < 10) {
print(count)
count <- count + 1
}
z <- 5
while(z >= 3 && z <= 10) {
print(z)
coin <- rbinom(1, 1, 0.5) ## Binomial distribution, rbinom(n, size, prob)
if(coin == 1) { ## random walk
z <- z + 1
} else {
z <- z - 1
}
}
接著進入while，程式碼第一行會先印出出5分，然後丟硬幣，丟硬幣這過程符合二項分配，正反機率各為0.5。
進入if，當硬幣為正面（1），給z加1分，然後回到while的第一行印出6分，再次進入if迴圈。
若硬幣為反面（0），給z減1分，然後回到while的第一行印出4分，再次進入if迴圈。
那什麼時候會退出迴圈呢？只要傳回while第一行的值在 z < 3 或 z > 10 的時候，while迴圈就會結束！
我試了3次，結果如下：
[1] 5
[1] 4
[1] 3
[1] 5
[1] 4
[1] 3
[1] 4
[1] 3
[1] 4
[1] 3
[1] 5
[1] 4
[1] 5
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 8
[1] 9
[1] 10
[1] 9
[1] 8
[1] 7
[1] 8
[1] 9
[1] 10
[1] 9
[1] 8
[1] 9
[1] 8
[1] 7
[1] 6
[1] 7
[1] 6
[1] 5
[1] 6
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
Control Structures - Repeat, Next, Break
x0 <- 1
tol <- 1e-8
repeat {
x1 <- computeEstimate() ## computeEstimate() is not a real function
if(abs(x1 - x0) < tol) { ## abs(x) computes the absolute value of x
break
} else {
x0 <- x1
}
}
for(i in 1:100) {
if(i <= 20) {
## Skip the first 20 iterations
next
}
## Do something here
}
Your First R Function
add2 <- function(x, y) {
x + y
}
above10 <- function(x) {
use <- x > 10
x[use]
}
above <- function(x, n = 10) {
use <- x > n
x[use]
}
columnmean <- function(y, removeNA = TRUE) {
nc <- ncol(y)
means <- numeric(nc)
for(i in 1:nc) {
means[i] <- mean(y[, i], na.rm = removeNA)
}
means
}
Functions (part 1)
f <- function(<arguments>) {
## Do Something interesting
}
Function Arguments
- formal arguments：包含在函數定義裡的參數。
- formals() 會將一個函數作為輸入，並返回函數內的一組包含所有formal arguments的列表。
- 在R，不是每次調用函數都會用到所有的formal arguments。也就是說，如果一個函數裡有10個不同的參數，你不需要指定每個參數的值，當使用者沒有明確給參數賦值時，它們的值為缺省值（defaults）。
Arguments Matching
> mydata <- rnorm(100)
> sd(mydata)
> sd(x = mydata)
> sd(x = mydata, na.rm = FALSE)
> sd(na.rm = FALSE, x = mydata)
> sd(na.rm = FALSE, mydata) ## 儘管結果不變，但調換順序容易引起誤解> atgs(lm)
function (formula, data, subset, weights, na.action, method = "qr", model = TRUE,
x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE, contrasts = NULL,
contrasts = NULL, offset, ...)
lm(data = mydata, y ~ x, model = FALSE, 1:100) ## y ~ x <- formula, 1:100 <- subset
lm(y ~ x, mydata, 1:100, model = FALSE)如果你只想指名剩餘參數的其中一個，直接通過命名的方式是最快的。
所有函數參數可以進行部分匹配，這功能的好處會體現在交互式工作上。R的操作順序如下：
- 檢查完全分配。如果你指定一個參數名，R會檢查是否有參數跟該命名完全匹配。
- 如果找不到，它會檢查是否有部分匹配（partial match）。
- 如果都沒有，則進行位置匹配（positional match）。
Functions (part 2)
f <- function(a, b = 1, c = 2, d = NULL) {
}
Lazy Evaluation
f <- function(a, b) {
a ^ 2
}
f(2)
## [1] 4
f <- function(a, b) {
print(a)
print(b)
}
f(45)
## [1] 45
## Error: argument "b" is missing, with no default
The “...” Argument
當你在擴展其他函數，且你希望複製整個參數時，你可以使用它。
myplot <- function(x, y, type = "1", ...) {
plot(x, y, type = type, ...)
}
> mean
function (x, ...)
UseMethod("mean")paste()的作用是將一組字符串連起來，來新建一個字符串或是字符向量，這函數會用到許多參數且無法預先注明有多少參數需要連接。
cat()的功能與paste()相似，它將一組字符串連起來，然後將連接後的字符串輸出到控制台或是文件。
> 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 = "i")
[1] "a:b"
> paste("a", "b", se = "i")
[1] "a b :"Scoping Rules - Symbol Binding
> search()
[1] ".GlobalEnv" "tools:rstudio" "package:stats" "package:graphics"
[5] "package:grDevices" "package:utils" "package:datasets" "package:methods"
[9] "Autoloads" "package:base"
當你用library()時，你載入的包會插入在".GlobalEnv"後面，也就是第二位置。
需注意的是，函數和非函數在R中有不同的命名空間。
Scoping Rules
這是R語言與原始的S語言最主要區別，作用域規則決定了一個值如何與函數中的自由變量綁定。在一個函數中，有兩種類型的變量：其中一種是函數的參數，它通過函數定義傳入函數；另一種是存在於函數中，但非參數的變量或符號。
問題是，你如何為這些符號賦值？
R用的是詞法作用域（Lexical Scoping），或稱靜態作用域（Static Scoping），它們是除動態作用域（Dynamic Scoping）之外最常見的選擇，與作用域規則相關的是：R用搜索列表給符號賦值的方法。
詞法作用域的優點之一是：它能簡化運算，尤其是在統計方面上十分有效。
f <- function(x, y) {
x ^ 2 + y / z
}
在定義函數的環境中搜索自由變量的值，什麼是環境呢？
環境是符號值對（pairs）的集合，你可以把R中所有的東西都想像成成對的符號和值。
每個環境，也就是成對的符號與值的組合，都有一個父環境，有點像一個環境的上層環境，這個環境繼承自它的父環境，面對一個環境來說，它可能有多個子環境。
只有一個環境沒有父環境，這環境是空環境。
R使用了很多這類環境，你可以把全局環境，也就是你的工作空間，看成一系列符號值對。
A function + an environment = a closure of function closure在R中，如果你將一個函數與環境聯繫起來，就創建了一個閉包（closure），或稱函數閉包。這些閉包是R中各種有趣操作的關鍵。
如果你在函數裡遇到一個自由變量，會怎樣呢？
你首先要找的是，這個函數是在哪個環境中被定義的。
舉例，如果我是在全局環境中定義了這個函數，那它自然就是在全局環境中被定義。一旦我在這個函數內無法得出它的值，我就會在全局環境中搜索，如果什麼都找不到，那會繼續在全局環境的父環境中搜索。
在全局環境外定義函數也是有可能的，一般而言，一個函數會在定義它的環境中搜索一個值，如果在那沒有找到，就會去它的父環境搜索，以此類推，直到所謂的頂層環境。
通常，頂層環境就是全局環境，但如果函數是在包中定義的，那頂層環境就是那個包的命名空間。如果在那仍沒有找到值，搜索會沿著搜索列表直到空環境。如果在這些環境中都找不到符號，就會報錯。
Scoping Rules - R Scoping Rules
為什麼作用域規則如此重要呢？
通常，在全局環境中定義了一個函數，這樣在用戶的工作區就能找到自由變量的值。
但在R中，你能在函數裡再定義其他函數，一個函數的返回值可以是一個函數。而作為返回值的這個函數，定義它的環境就不是全局環境，而是在另一個函數內部。
make.power <- function(n) {
pow <- function(x) {
x ^ n ## 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] 2y <- 10
f <- function(x) {
y <-2
y ^ 2 + g(x)
}
g <- function(x) {
x * y
}有一種情況下，詞法作用域和動態作用域的結果看起來是一樣的，就是在全局環境中定義函數，並也是在全局環境中調用該函數。
> g <- function(x) {
+ a <-3
+ x + a + y
+ }
> g(2)
Error in g(2) : object "y" not found
> y <- 3
> g(2)
[1] 8還有很多語言支持詞法作用域：
- Python
- Perl
- Scheme
- Common Lisp (all languages converge to Lisp)
因為詞法作用愈的本質和環境的複雜性，以及它們的關聯方式。很難在物理內存以外的地方實現這樣的模型。
每個函數都有一個指針指向它的定義環境，而這個定義環境可以是任何地方，因為函數裡可再定義函數，如果一個函數的返回值是另一個函數，就要有個指針指向儲存這個環境的那部份內存，這讓模型更複雜，但也更有用。
Scoping Rules - Optimization Example (OPTIONAL)
R中有一些用於優化的函數，optim()、nim()和optimize()，向這些函數傳入一個目標函數，這目標函數的參數是向量。
但在統計學中，我們想要最小化或最大化的這目標函數，例如log-likelihood，往往和其他東西有關，而不僅僅依賴於那些你在最大化時能調整的參數值。它尤其依賴數據之類的東西。
另外，當你做這些優化時，大多數的情況，固定某些參數值是非常有用的。
你構造一個目標函數，把目標函數所依賴的所有數據及其他東西，都包含在目標函數的定義環境中，這樣你就不必每次調用這函數的時候都要對這些東西命名，你唯一要做的就是命名參數的值。
下面是一個構造函數，構造一個負的log-likelihood。
※R中大多數的優化函數，像optim()、nim()和optimize()，都默認求函數的最小值，因此如果你想取最大值，就應該取它們的負值。
make.NegLogLik <- function(data, fixed = c(FALSE, FALSE)) { ## fixed 決定是否要固定某些參數
params <- fixed
function(p) { ## 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)
}
}> set.seed(1); normals <- rnorm(100, 1, 2) ## set.seed() 用來設定同一組亂數
> 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)
}
<environment: 0x165bla4> # pointer
> ls(environment(nLL))
[1] "data" "fixed" "params"> optim(c(mu = 0, sigma = 1), nLL)$par
mu sigma
1.218239 1.787343
> nLL <- make.NegLogLik(normals, c(FALSE, 2)) ## Fixed σ = 2
> optimize(nLL, c(-1, 3))$minimum ## optimize() 只能最小化單變量函數
[1] 1.217775
> nLL <- make.NegLogLik(normals, c(1, FALSE)) ## Fixed μ = 1
> optimize(nLL, c(1e-6, 10))$minimum
[1] 1.800596nLL <- 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
Coding Standards
如果你是第一次編寫程式，有幾個一貫的標準是我們會去注意的：
- 永遠使用文本（text files）／文字編輯器（text editor）。
- 縮排你的編碼，重要的是使你的編碼更容易讓別人理解。
- 限制你編碼的寬度，同第二點，這麼做是為了使別人更容易理解。
- 限制函數的長度，盡量將一個函數的功能限制在一種基本操作。
在Rstudio中，你可以設定你縮排（Tab）的長度。點選功能表列的Tools，點擊Global Options，在左欄點擊Code。
第一個選項你可以改變你的縮排距離（Tab width）。
在第四個選項你可以限制你編碼的寬度（Margin column）為幾行。
Dates and Times
Date: 表示日期。內存是以1970-01-01至今的天數來儲存的。
POSIXct / POSIXit: 表示時間。內存是以1970-01-01至今的秒數來儲存的。
POSIXct: 時間用非常大的整數來表示的，如果你想把時間存在數據框之類的東西，這種類是很有用的。
POSIXit: 把時間當作列表來儲存，它還存有一些跟給定時間相關的訊息，例如：這個時間是星期幾？是一年中的第幾天？是幾月？是幾號？
另外還有一些泛用函數有關日期和時間的：
weekdays(): 告訴你給定的時間或日期是星期幾。
months(): 告訴你給定的時間或日期是在幾月。
quarters(): 告訴你給定的時間或日期楚處於第幾季度（Q1, Q2, Q3, Q4）。
x <- as.Date("1970-01-01")
x
## [1] "1970-01-01"
unclass(x)
## [1] 0
unclass(as.Data("1970-01-01"))
## [1] 1
x <- Sys.time()
x
## [1] "2015-06-13 19:30:47 CST"
p <- as.POSIXlt(x)
names(unclass(p))
## [1] "sec" "min" "hour" "mday" "mon" "year" "wday" "yday" "isdst"
## [10] "zone" "gmtoff"
unclass(p)
## $sec
## [1] 47.73445
## $min
## [1] 30
## $hour
## [1] 19
## $mday
## [1] 13
## $mon
## [1] 5
## $year
## [1] 115
## $wday // 禮拜幾
## [1] 6
## $yday // 今年的第幾天
## [1] 163
## $isdst // 夏令時間
## [1] 0
## $zone
## [1] "CST"
## $gmtoff
## [1] 28800
## attr(,"tzone")
## [1] "" "CST" "CDT"
x <- Sys.time()
x ## Already in 'POSIXct' format
## [1] "2015-06-13 19:30:47 CST"
unclass(x)
## [1] 1434195048 // 它是自1970年1月1日至今的秒數
> x <- as.Date("2012-03-01")
> y <- as.Date("2012-02-28")
x - y
## Time difference of 2 days
> x <- as.POSIXct("2012-10-25 01:00:00")
> y <- as.POSIXct("2012-10-25 01:00:00", tz = "GMT")
y - x
## Time difference of 8 hours