標籤

Coursera R Programming Week4 心得筆記

R 語言程序開發
約翰霍普金斯大學 公共衛生學院
R Programming
Johns Hopkins Bloomberg School of Public Health

Week 4
第四週學習筆記
2015年6月1日 - 2015年6月29日

Week 4: Simulation and Profiling

This week covers how to simulate data in R, which serves as the basis for doing simulation studies. We also cover the profiler in R which lets you collect detailed information on how your R functions are running and to identify bottlenecks that can be addressed. The profiler is a key tool in helping you optimize your programs. Finally, we cover the str function, which I personally believe is the most useful function in R.

Learning Objectives

By the end of this week you should be able to:
  • Call the str function on an arbitrary R object
  • Describe the difference between the "by.self" and "by.total" output produced by the R profiler
  • Simulate a random normal variable with an arbitrary mean and standard deviation
  • Simulate data from a normal linear model

Programming

There is a graded programming assignment for this week.
  • Programming assignment 3: Hospital Quality

Week 4
The str Function

srt()可以緊湊地顯示R對象的內部結構，srt代表的就是structure。
  • 它是一個非常通用的診斷函數，你可以把它當作是summary()的替代品。
  • 它特別適用於緊湊地顯示大型列表，例如那些嵌套有其他列表的列表。
  • 它的目標是對每個基本對象，輸出大約一行的內容。如果你傳遞給它一個簡單的對象，例如向量，它將返回給你一行的輸出備份，並將結果印在控制台上。
srt()的基本目的是回答這樣一個問題：對象裡有什麼？
> str(str)
function (object, ...)

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

> str(ls)
function (name, pos = -1L, envir = as.environment(pos), all.names = FALSE, 
    pattern, sorted = TRUE) 

> x <- rnorm(100, 2, 4)
> summary(x)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-11.2800  -0.4629   2.4680   2.0820   4.9030  12.2300 
> str(x)
 num [1:100] 0.179 -0.447 7.401 -4.591 4.92 ...

> f <- gl(40 ,10)
> str(f)
 Factor w/ 40 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
> summary(f)
 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 
10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 
26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 
10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 

> 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
> str(airquality)
'data.frame': 153 obs. of  6 variables:
 $ Ozone  : int  41 36 12 18 NA 28 23 19 8 NA ...
 $ Solar.R: int  190 118 149 313 NA NA 299 99 19 194 ...
 $ Wind   : num  7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
 $ Temp   : int  67 72 74 62 56 66 65 59 61 69 ...
 $ Month  : int  5 5 5 5 5 5 5 5 5 5 ...
 $ Day    : int  1 2 3 4 5 6 7 8 9 10 ...

> m <- matrix(rnorm(100), 10, 10)
> str(m)
 num [1:10, 1:10] -0.833 -0.555 -0.886 0.215 0.858 ...
> m[, 1]
 [1] -0.83330357 -0.55521696 -0.88617367  0.21543180  0.85833214 -0.03012694
 [7] -0.27404048  0.20569141 -0.45482943  0.55274505

> s <- split(airquality, airquality$Month)
> str(s)
List of 5
 $ 5:'data.frame': 31 obs. of  6 variables:
  ..$ Ozone  : int [1:31] 41 36 12 18 NA 28 23 19 8 NA ...
  ..$ Solar.R: int [1:31] 190 118 149 313 NA NA 299 99 19 194 ...
  ..$ Wind   : num [1:31] 7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
  ..$ Temp   : int [1:31] 67 72 74 62 56 66 65 59 61 69 ...
  ..$ Month  : int [1:31] 5 5 5 5 5 5 5 5 5 5 ...
  ..$ Day    : int [1:31] 1 2 3 4 5 6 7 8 9 10 ...
 $ 6:'data.frame': 30 obs. of  6 variables:
  ..$ Ozone  : int [1:30] NA NA NA NA NA NA 29 NA 71 39 ...
  ..$ Solar.R: int [1:30] 286 287 242 186 220 264 127 273 291 323 ...
  ..$ Wind   : num [1:30] 8.6 9.7 16.1 9.2 8.6 14.3 9.7 6.9 13.8 11.5 ...
  ..$ Temp   : int [1:30] 78 74 67 84 85 79 82 87 90 87 ...
  ..$ Month  : int [1:30] 6 6 6 6 6 6 6 6 6 6 ...
  ..$ Day    : int [1:30] 1 2 3 4 5 6 7 8 9 10 ...
 $ 7:'data.frame': 31 obs. of  6 variables:
  ..$ Ozone  : int [1:31] 135 49 32 NA 64 40 77 97 97 85 ...
  ..$ Solar.R: int [1:31] 269 248 236 101 175 314 276 267 272 175 ...
  ..$ Wind   : num [1:31] 4.1 9.2 9.2 10.9 4.6 10.9 5.1 6.3 5.7 7.4 ...
  ..$ Temp   : int [1:31] 84 85 81 84 83 83 88 92 92 89 ...
  ..$ Month  : int [1:31] 7 7 7 7 7 7 7 7 7 7 ...
  ..$ Day    : int [1:31] 1 2 3 4 5 6 7 8 9 10 ...
 $ 8:'data.frame': 31 obs. of  6 variables:
  ..$ Ozone  : int [1:31] 39 9 16 78 35 66 122 89 110 NA ...
  ..$ Solar.R: int [1:31] 83 24 77 NA NA NA 255 229 207 222 ...
  ..$ Wind   : num [1:31] 6.9 13.8 7.4 6.9 7.4 4.6 4 10.3 8 8.6 ...
  ..$ Temp   : int [1:31] 81 81 82 86 85 87 89 90 90 92 ...
  ..$ Month  : int [1:31] 8 8 8 8 8 8 8 8 8 8 ...
  ..$ Day    : int [1:31] 1 2 3 4 5 6 7 8 9 10 ...
 $ 9:'data.frame': 30 obs. of  6 variables:
  ..$ Ozone  : int [1:30] 96 78 73 91 47 32 20 23 21 24 ...
  ..$ Solar.R: int [1:30] 167 197 183 189 95 92 252 220 230 259 ...
  ..$ Wind   : num [1:30] 6.9 5.1 2.8 4.6 7.4 15.5 10.9 10.3 10.9 9.7 ...
  ..$ Temp   : int [1:30] 91 92 93 93 87 84 80 78 75 73 ...
  ..$ Month  : int [1:30] 9 9 9 9 9 9 9 9 9 9 ...
  ..$ Day    : int [1:30] 1 2 3 4 5 6 7 8 9 10 ...

Simulation - Generating Random Numbers

有一些函數，可以用於模擬已知機率分布中的數字或變量：
  • rnorm: 根據給定的均值和標準差生成隨機常態分布。
  • dnorm: 計算常態機率密度（Normal probability density）。
  • pnorm: 計算常態分布的累積分布函數（cumulative distribution function）。
  • rpois: 根據已知平均發生率（rate）生成隨機泊松分佈（Poisson distribution）。
有四種基本的函數與此相關：
  • d: 密度。
  • r: 隨機生成。
  • p: 累積分布。
  • q: 分位數函數。 
對於常態分佈來說，如果你不指定均值和標準差，它的預設會是標準常態分布，也就是mean = 0, sd = 1。
dnorm(x, mean = 0, sd = 1, log = FALSE)
pnorm(q, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)  ## lower.tail: 左尾
qnorm(p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
rnorm(n, mean = 0, sd = 1)

如果Φ是標準常態分布的累積分布函數，那pnorm()就等於Φ，qnorm()就等於Φ的反函數（inverse）。
> x <- rnorm(10)
> x
 [1] -0.5265158  1.7975343 -0.3502657  0.1656121  0.1395387 -1.3202982
 [7] -0.4113113  0.9013153  1.5222450 -1.3703012
> x <- rnorm(10, 20, 2)
> x
 [1] 20.50187 18.50558 23.53972 18.33857 19.88685 22.79343 21.88073 17.99585
 [9] 20.04389 21.91555
> summary(x)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  18.00   18.85   20.27   20.54   21.91   23.54 


> 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
> ppois(2, 2)  ## Cumulative distribution
[1] 0.6766764  ## Pr(x <= 2)
> ppois(4, 2)
[1] 0.947347   ## Pr(x <= 4)
> ppois(6, 2)
[1] 0.9954662  ## Pr(x <= 6)

Simulation - Simulating a Linear Model
> 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.4080 -1.5400  0.6789  0.6893  2.9300  6.5050 
> plot(x, y)

> set.seed(10)
> x <- rbinom(100, 1, 0.5)  ## binary random variable (n = 1, p = 0.5)
> e <- rnorm(100, 0, 2)
> y <- 0.5 + 2 * x + e
> summary(y)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-3.4940 -0.1409  1.5770  1.4320  2.8400  6.9410 
> plot(x, y)

> set.seed(1)
> x <- rnorm(100)
> log.mu <- 0.5 + 0.3 * x
> y <- rpois(100, exp(log.mu))
> summary(y)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00    1.00    1.00    1.55    2.00    6.00 
> plot(x, y)
Simulation - Random Sampling

 smaple()允許你從指定的一組對象集合中隨機抽樣，如果你給它一個數字向量，它允許你從這一組向量隨機抽取一個樣本。通過指定某類對象的向量，你幾乎可以創建任意你想要的分布，然後從中抽樣。
> set.seed(1)
> sample(1:10, 4)  ## without replacement
[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)     ## permutation
 [1]  4  7 10  6  9  2  8  3  1  5
> sample(1:10)
 [1]  2  8  6  9  1  4 10  3  5  7
> sample(1:10, replace = TRUE)  ## Sample w/replacement
 [1] 2 3 4 1 4 9 4 5 6 5

總結
  • 你可以使用r*函數從指定的機率分布中隨機取樣。像是rnorm()、rpois()、rbinom()。
  • 所以可能用到的標準分布都已經內嵌在R裡面。Normal, Poisson, Binomial, Exponential, Gamma, etc.
  • sample()用來從任意向量中隨機取樣。
  • 只要你想模擬，都一定要記得設置種子，這樣你才可以重塑你之前得到的結果。

R Profiler

當你在開發一個規模較大的程序，或者做很大型的資料分析時，R裡面的分析器是一個非常方便的工具。基本上，運行R的代碼通常需要很時間，分析器（profiler）就是一個可以確切分析為何耗時如此之長的方便工具，並提供一些建議來解決你的問題。
“We should forget about small efficiencies, say about 97% of the time: premature optimization is the root of all evil.”
──Donald Knuth
基本的是，你應該先設計好代碼，讓代碼易懂，在你寫出可以運行的代碼之後，再嘗試去優化它。

 Using system.time()

 system.time()的作用是，取出並分析任意一條R語句，然後函數會告訴你運行那個語句所需的時間。當你在計算機上執行語句時，有兩個關於時間的重要概念：
  • user time: 運行這行語句需要占用中央處理處（CPU(s)）的時間。
  • elapsed time: 有時也叫做wall clock time，是你所經歷的時間，所以即使你就是用戶，你經歷的也不是user time，你經歷的是elapsed time。
> system.time(readLines("http://www.jhsph.edu"))
   user  system elapsed 
   0.004  0.002   0.431    ## Elapsed time > user time
> hilbert <- function(n) {
          i <- 1:n
          1 / outer(i - 1, i, "+")
  }
> x <- hilbert(1000)
> system.time(svd(x))      ## svd() singular value decomposition,
                           ## use the accelerate framework on the Mac.
   user  system elapsed 
  1.605   0.094   0.742    ## Elapsed time < user time


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

但如果你不知道問題可能在哪裡？從哪開始檢查時，怎麼辦？
  • Rprof()
    • R必須再分析器的支持下進行編輯，不過並不是所有情況下都這樣，但基本上絕大部分的情況都是如此。
  • summaryRprof()
    • 它會取出分析器的結果，並總結成可以閱讀的方式，因為分析器的原始輸出基本上是用不了的。
需要注意的是，你不能同時使用system.time()和Rprof()，它們本來就不是為了共同使用而被設計的。

Rprof()可以在規律的樣本區間內紀錄追蹤函數調用線（function call stack），也就是說當你的函數運行時，它會查詢函數調用線，看有多少函數調用了其他函數，而這些函數又調用了那些其他的函數，它會把它們都印出。基本上，它所做的就是在非常短的區間內打印出輸出函數的調用線。
每過0.02秒，它就會打印出輸出函數的調用線。
如果你的函數運行少於0.02秒，那R的分析器就沒什麼用了。 
> system.time({
     n <- 1000
     r <- numeric(n)
     for (i in 1:n) {
         x <- rnorm(n)
         r[i] <- mean(x)
     }
  })
   user  system elapsed 
  0.097   0.002   0.099
summaryRprof()
當你看到函數調用線，你就知道每一行間相隔了0.02秒，這是獲得取樣的頻率。既然你可以計算每個函數需要的時間，如果它出現在函數調用線裡，那麼你必然在其上花費了一些時間。這裡有兩種方法，可以將你從R分析器中得到的數據規範化：
  • by.total: 把每個函數使用的時間依照總運行時間排序。
  • by.self: 作用相同，但它會先減去調用線中之前的函數耗費的時間。
基本思路就是：通過規範化（normalizing）在某個函數上耗費的總時間，你就可以了解時間的耗費，以及函數被調用的次數。
假設你調用的函數是lm()，那你便把100%的時間都用在這個函數上，因為它位於頂層（top level），但事實上，頂層函數通常不會真正做任何重要的事情，它們只是會調用那些真正做事的輔助函數。
因此，知道頂層函數使用了多少時間並沒有用，因為那並非實際在工作的函數。大多數的時候，你會發現當你排除所有被調用的低級函數之後，頂層函數只會耗費非常少的時間，因為所有的工作及計算，都由低級函數完成的。
因此我認為by.self()的格式是最有用的，它可以告訴你在減去所有被調用的低級函數耗費的時間之後，一個給定的函數所耗費的時間。

標籤：

熱門文章