Subscribed unsubscribe Subscribe Subscribe

R軟體再談行導向

R軟體閒話家常

進行資料分析前,資料正規化經常是一個必要的步驟。此時若能留意R 的行導向特性,當可事半功倍。例如:欲對下列三行變數做最大-最小正規化

> (x <- matrix(1:24,8,3))

       [,1]  [,2]  [,3]

[1,]    1    9    17

[2,]    2   10   18

[3,]    3   11   19

[4,]    4   12   20

[5,]    5   13   21

[6,]    6   14   22

[7,]    7   15   23

[8,]    8   16   24

首先以先前提過的apply 函數求出各行最小值,

> (y <- apply(x, 2, min))

[1]  1  9 17

再求出各行最大值,

> (z <- apply(x, 2, max))
[1] 8 16 24

y 與z 看起來像是列向量,實際上它們是行向量!因此x 須先轉置(transpose),才能將各行剪去其最小值。

>  t(x) - y

     [ ,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]

[1,]    0    1     2     3    4     5     6    7

[2,]    0    1     2     3    4     5     6    7

[3,]    0    1     2     3    4     5     6    7

讀者可試著輸入x 未轉置的指令,結果是x 減去[1 9 17 1 9 17 1 9]T,因為y 的長度不及8,所以R 會將較短向量y 中元素自動重複,直到其長度為8 為止。

> x - y
       [,1]   [,2]   [,3]

[1,]     0    -8      8

[2,]    -7     9      1

[3,]  -14     2     18

[4,]     3    -5     11

[5,]    -4    12      4

[6,]  -11     5      21

[7,]    6     -2     14

[8,]   -1    15       7

我們所要的正規化結果是將(t(x) - y)再除以各行最大值與最小值的差(z - y),

> (t(x) - y)/(z - y)

        [,1]           [,2]             [,3]             [,4]             [,5]            [,6]             [,7]  [,8] 

[1,]     0  0.1428571  0.2857143  0.4285714  0.5714286  0.7142857  0.8571429    1

[2,]     0  0.1428571  0.2857143  0.4285714  0.5714286  0.7142857  0.8571429    1

[3,]     0  0.1428571  0.2857143  0.4285714  0.5714286  0.7142857  0.8571429    1

再將它轉置就回到原始資料8 列3 行的形式了!

> t*1

                [,1]            [,2]             [,3]

[1,] 0.0000000  0.0000000  0.0000000

[2,] 0.1428571  0.1428571  0.1428571

[3,] 0.2857143  0.2857143  0.2857143

[4,] 0.4285714  0.4285714  0.4285714

[5,] 0.5714286  0.5714286  0.5714286

[6,] 0.7142857  0.7142857  0.7142857

[7,] 0.8571429  0.8571429  0.8571429

[8,] 1.0000000  1.0000000  1.0000000

 

撰文者:
鄒慶士 博士
現任:
北商資訊與決策科學所教授
中華 R 軟體學會理事長
信箱:
vince.tsou@gmail.com

*1:t(x) - y)/(z - y

巨人的肩膀與 plyr R軟體套件

R軟體閒話家常

瞎忙了一陣子,驚覺許久未動筆撰寫閒話家常了,雖然每天都會操作或閱讀一些 R 相關文件,但總覺不能脫稿太久。回想起常在課堂上對學生或學員們說:「 R 海無涯,為勤是岸。」不但是鼓勵學員,也算是一種自我惕勵。但這樣說來,好像學 R 是件令人可懼的事,其實讀者們不用太擔心!因為 R 玩家中有許多強大的巨人會引領我們向前航行。比如說今天的主角<br />plyr 套 件 ,就 是由多 產的 美國 Rice 大學統 計系助理 教授,紐西蘭裔的 Hadley Wickham(http://had.co.nz/)所開發的,他是資料視覺化(data visualization)的專家,也是知名的資料整理(data manipulation)套件 reshape 與繪圖套件 ggplot2 的主要開發者。

plyr 的讀音同英文的老虎鉗(pliers)去掉字尾 s,顧名思義它是運用分割-應用-組合(split-apply-combine)的策略於資料分析的工作上,也就是說將一大的問題,切割為可管控的子集,個別進行問題解決後,再重新組合起來。舉例來說,先建一個包含三個數值變數與兩個因子的72 筆資料,

> library(plyr)

> dd <- data.frame(matrix(round(rnorm(216), 2), 72, 3), c(rep("A", 24), rep("B", 24), rep("C", 24)), c(rep("J", 36), rep("K", 36)))

> colnames(dd) <- c("v1", "v2", "v3", "dim1", "dim2")

> head(dd)

          v1        v2       v3   dim1   dim2
1     0.32    -0.55    2.52         A        J
2    -1.13     1.27   -0.19         A        J
3     0.40     1.51   -0.30         A        J
4     0.40     0.65    1.37         A        J
5    -0.93     0.67   -0.09         A       J
6    -0.05     0.17    0.71         A       J

plyr 套件中的 ddply 函數可以不用撰寫迴圈,迅速地計算出兩因子 dim1 與 dim2 各種可能組合下,數值變數 v1 的平均數。

> ddply(dd, c("dim1", "dim2"), function(df) mean(df$v1))

    dim1   dim2                 V1
1        A        J    -0.05333333
2        B        J    -0.15416667
3        B        K     0.14666667
4        C        K     0.17041667

其實它的功能如同 SQL 中的某些指令,或是可看成 EXCEL 的樞紐分析。這對彈性強大的 R來說不是一件新鮮事,R 玩家常掛在嘴邊的一句話就是:「幾乎沒有 R 完成不了的資料分析工作!」上述計算的另一種 ddply 函數寫法如下:

> ddply(dd, c("dim1", "dim2"), summarise, V1 = mean(v1))

    dim1   dim2                 V1
1        A        J    -0.05333333
2        B        J    -0.15416667
3        B        K     0.14666667
4        C        K     0.17041667

ddply 的頭兩個英文字母分別代表函數的輸入為資料框(data frame),輸出結果亦為資料框。當然,plyr 套件還有其它的函數,ldply 表輸入為串列(list),輸出為資料框;alply 表輸入為陣列(array),輸出為串列;d_ply 則不傳回任何物件,通常此函數輸出到圖形裝置或檔案,更多的資訊請參考 plyr 手冊或下面的參考資料。

 

參考資料:Wickham, H. (2011), &ldquo;The split-apply-combine strategy for data analysis&rdquo;, Journal of Statistical Software, Vol. 40, Issue 1, pp. 1-29.

  

撰文者:
鄒慶士 博士
現任:
北商資訊與決策科學所教授
中華 R 軟體學會理事長
信箱:
vince.tsou@gmail.com

 

 

R軟體應用於機率密度曲線與累積分配曲線的繪製

R軟體閒話家常

機率與統計的教科書都會有機率密度函數(probability density function, pdf) 與累積分配函數(cumulative distribution function, cdf)的圖形,以說明累積機率值與分位點(quantile)的對應關係。在R 中,讀者可以pnorm 函數求出某一分位點的累積幾率值。例如:標準常態分位點z = 1.5 的累積機率值如下。

> pnorm(1.5, 0, 1)

[1] 0.9331928

 

反過來說,qnorm 可以求出累積機率為0.9331928 所對應的標準常態分位點。

> qnorm(0.9331928, 0, 1)

[1] 1.5

 

為了讓讀者深刻瞭解前述的對應關係,先以R 中的curve 函數繪製標準常態的密度曲線:

> curve(dnorm(x, 0, 1), xlim = c(-3, 3), main = "Area to left of 1.5 = pnorm(1.5, 0, 1)")

再利用多邊形繪製函數polygon,連接密度相當高的點(橫座標從-3 到1.5)形成灰色陰影區域,並在適當的位置加上說明文字。

> cord.x <- c(-3, seq(-3, 1.5, 0.01), 1.5)

> cord.y <- c(0, dnorm(seq(-3, 1.5, 0.01)), 0)

> polygon(cord.x, cord.y, col = "grey")

> text(-0.5, 0.1, "area=pnorm(1.5,0,1)")

>text(1.5, 0.02, "1.5=qnorm(0.9331928, 0 ,1)")

f:id:DataScience:20140817130331p:plain

 

同樣地,R 可繪製累積分配函數的圖形,再利用abline 函數加上適當的水平與垂直線、以及說明文字。走筆至此,冀望讀者能感受到R 強大且高品質的繪圖輸出,更多的繪圖功能容後再敘。

> curve(pnorm(x, 0, 1), xlim = c(-3,3), main = "F(z)=P(Z<=z)")

> abline(h = 0.9331928, lty = 2)

> abline(v = 1.5, lty = 3)

> text(-2, 0.85, "pnorm(1.5) is 0.9331928")

> text(1.5, 0.02, "qnorm( 0.9331928) is 1.5")

f:id:DataScience:20140817131027p:plain

 參考資料:

Verzani, J. (2001), simpleR &ndash; Using R for Introductory Statistics, The CSI MathDepartment, City University of New York.

 

撰文者:
鄒慶士 博士
現任:
北商資訊與決策科學所教授
中華 R 軟體學會理事長
信箱:
vince.tsou@gmail.com

R軟體的寬資料與長資料

R軟體閒話家常

當我們對單一資料蒐集對象測量多個變數值時,寬資料是指同一對象的所有測量值都排在同一列;長資料則是各個測量值單獨成一列,並標明其是哪個變數的測量值。某些統計分析須使用寬資料,例如相關分析;也有些分析偏好長資料,例如變異數分析。

 

首先,我們以data.frame 函數建一個五個蒐集對象、三個變數的寬資料:

> mydata <- data.frame(var1 = c(12, 15, 19, 22, 15), var2 = c(18, 12, 42, 29, 44), var3 = c(8, 17, 22, 19, 31))

>mydata

f:id:DataScience:20140817123740p:plain

運用stack 函數可將資料轉成長資料,

> sdata <- stack(mydata)

> sdata

f:id:DataScience:20140817123957p:plain

stack 函數還可使用select 參數選擇要轉換成長資料的變數,此處我們選var1 和var2。

> sdata <- stack(mydata, select = c(var1, var2))

讀者可利用下列指令檢查長資料十筆數據中均無var3!

>sdata$in == "var3"

[1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 

 

反過來說,unstack 函數可將長資料轉換為寬資料,要注意的是必須使用R 中的公式符號(formula notation) 'values ~ ind',此處解讀為依不同的ind 值將values 歸類。最後,對於更複雜的長寬資料格式轉換,讀者可運用R 中的reshape 套件,這個部分容後再敘。

 

> mydata <- unstack(sdata, values ~ ind)

> head(mydata)

f:id:DataScience:20140817124601p:plain

 

 參考資料:Spector, P. (2008), Data Manipulation with R, Springer.

撰文者:
鄒慶士 博士
現任:
北商資訊與決策科學所教授
中華 R 軟體學會理事長
信箱:
vince.tsou@gmail.com

R軟體for迴圈的迷思

R軟體閒話家常

許多人認為學會過了程式設計,就能夠駕馭R 語言!殊不知R 語法有其特殊性,更遑論背後的統計模型與繪圖專業了。舉例來說,下面就是標準的不良R 程式(canonical bad Rprogram),用來計算a 與b 兩向量的內積:

> (a <- 1:5) 

[1] 1 2 3 4 5

> (b <- 5:1)

[1] 5 4 3 2 1

 

> d <- NULL

> for ( i in 1:length(a)) {

+ d[i] <- a[i] * b[i]

+ }

>d

[1] 5 8 9 8 5

 

> s <- 0

> for (i in 1:length(d)) {

+ s <- s+d[i]

+ }

> s

[1] 35

傳統的兩次迴圈設計,當資料量大時,會大幅影響執行時間。其實善用R 的向量化(vectorization)特性,一到兩行的指令即可快速輕易地算出向量內積了。(註:因資料量很小,此例感受不到速度的差異!)

> (d <- a * b)

[1] 5 8 9 8 5

> (s <- sum(a * b)

[1] 35

 

以程式執行效率的觀點來看,R的世界非不得已 (指當無法運用apply系列函數與vectorization特性時) 才會使用for迴圈。

 

撰文者:
鄒慶士 博士
現任:
北商資訊與決策科學所教授
中華 R 軟體學會理事長
信箱:
vince.tsou@gmail.com

R軟體apply系列基本函數的運用

R軟體閒話家常

資料分析經常需要運用for 迴圈反覆執行某項工作,然而R 程式中卻不可用太多的迴圈,否則會大大降低程式執行的效率!以紐約市空氣品質資料為例,其為包含六個變數的154 筆資料,以apply 函數可以輕易地計算出各行變數的平均值。

> head(airquality)

f:id:DataScience:20140810131628p:plain

 

> apply(airquality, 2, mean) # apply(資料及名稱, 2表逐行,套用函數名稱

f:id:DataScience:20140810131917p:plain

 

結果中出現NA 的是因為該變數原始資料有遺缺值,此時可將遺缺值NA 排除計算(i.e. na.rm = T)即可。

> apply(airquality, 2, mean, na.rm=T)

f:id:DataScience:20140810132132p:plain

 

同理,讀者可以 apply(airquality, 1, mean)逐列求出154 筆觀測值的平均值(此統計值有意義嗎?)

apply 系列的基本函數還有lapply 與sapply,lapply 將airquality 資料視為六個元素的串列(list),逐一元素加總遺缺值個數,故無須加入逐列或逐行的參數值。

 

> lapply(airquality, function(x) sum(is.na(x)))

   f:id:DataScience:20140810132626p:plain

sapply 可將上面結果以較簡單的矩陣形式展現,此即為函數名字首的含義。

> sapply(airquality, function(x) sum(is.na(x)))

f:id:DataScience:20140810132851p:plain

 

sapply 函數也可以判斷哪些變數是因子,或加入自定函數計算因子變數的水準數,讀者當可想出apply系列基本函數的其它應用。最後,善用apply 方可將R 程式技巧提升到更高的層次。

>sapply(InsectSprays, is.factor)

f:id:DataScience:20140810133424p:plain

> sapply(insectSprays, sunction(x) if (!is.factor(x)) return(0) else length(levels(x))) 

f:id:DataScience:20140810133628p:plain

 

 參考資料:

 Maindonald, J.H. (2004), Using R for Data Analysis and Graphics: Introduction, Code,and Commentary, http://www.stats.uwo.ca/DAAG/.

 

撰文者:
鄒慶士 博士
現任:
北商資訊與決策科學所教授
中華 R 軟體學會理事長
信箱:
vince.tsou@gmail.com

R軟體如何儲存隨機種子?

R軟體閒話家常

R是一個很好的模擬工具,有許多內建函數可資運用。使用者有時希望獲得相同的模擬結果,則可將隨機種子(random seed)儲存起來重複利用。例如:

> set.seed(12)
> seed <- .Random.seed
> x <- runif(10)
> .Random.seed <- seed
> y <- runif(10) x-y

[1] 0 0 0 0 0 0 0 0 0 0

 

讀者會發現兩次模擬的十筆均勻分配資料是一樣的。

如果要在兩次 R sessions中獲得相同的模擬結果,首先請輸入:

> # seesion 1

> set.seed(10)

> seed <- .Random.seed

> x<-runif(10)

 

將第一個 session 的結果儲存起來後離開 R:

> write.table(x, "x")

>write.table(seed, "seed")

> q( )

 

重新開啓 R,在第二個session 輸入:

> # session 2

> x <- read.table("x")

> seed <- read.table("seed")

> .Random.seed <- t(seed)

> y <- runif(10)

> x - y

                    x

1       3.330669e-16

2       3.330669e-16

3       -2.220446e-16

4       -4.440892e-16

5       2.775558e-17

6       -1.110223e-16

7       2.220446e-16

8       -1.665335e-16

9       -3.330669e-16

10     1.110223e-16

讀者會發現兩次 sessions 模擬的結果幾乎沒有差異!

 

參考資料:http://lists.cs.wisc.edu/archive/stat-forum/2007-March/msg00004.shtm/

撰文者:
鄒慶士 博士
現任:
北商資訊與決策科學所教授
中華 R 軟體學會理事長
信箱:
vince.tsou@gmail.com