關卡 1

這門課程是要介紹R的線性代數運算系統。線性代數在許多現代的統計方法、DataMining方法上都很重要。所以理解R的線性代數系統,對於撰寫自己的演算法,以及了解OpenSource的演算法,是非常重要的。

關卡 2

我們先來複習R中的向量。請同學建立一個型態為numeric,由1到18的向量。

1:18

關卡 3

我們可以利用matrix這個函數來建立一個矩陣。舉例來說,matrix(1:18,6,3)就可以建立一個6乘3的矩陣。請同學試試看,並把這個矩陣寫入變數x。

x <- matrix(1:18, 6, 3)

關卡 4

要了解matrix,我們可以先打開matrix的說明頁面。

help(matrix)

關卡 5

根據matrix的說明頁面,請問以下哪些不是matrix的參數?

dim

關卡 6

數學上來說,一個矩陣除了值之外,需要的就是維度的定義。而matrix這個函數,data的參數代表矩陣的值,nrow的部份代表這個矩陣有多少列,而ncol的部份代表這個矩陣有多少行。所以matrix(1:18,6,3)會產生一個6列3行(簡稱6乘3)的矩陣。請同學輸入x看一下R是如何印出x這個6乘3的矩陣。

x

關卡 7

我們可以看看x的屬性。

attributes(x)

關卡 8

我們看到有一個名字叫做“dim”的屬性,值是c(6,3),代表這是一個6乘3的矩陣。

關卡 9

重要的屬性,R都會提供內建函數來方便存取。dim就是可以存取矩陣維度的函數。請同學試試看:dim(x)

dim(x)

關卡 10

也可以利用dim(x)<-c(3,6)來更改x的維度。請同學試試看。

dim(x) <- c(3, 6)

關卡 11

我們可以把x印出來看看。

x

關卡 12

R的除了矩陣(matrix)之外,還有高維度的陣列。舉例來說,我們可以直接更改x的維度到更高維:dim(x)<-c(3,3,2)。請同學試試看。

dim(x) <- c(3, 3, 2)

關卡 13

請同學把x印出來看看。

x

關卡 14

我們也可以用array(1:18,c(3,3,2))來建立一個高維矩陣。這裡很清楚可以看到:array的第二個參數就是維度。我們請同學試試看輸入array(1:18,c(3,3,2))

array(1:18, c(3, 3, 2))

關卡 15

在RBasic-02-Data-Structure-Vectors的課程中,我們有提過兩種從向量中挑選一部分的值得方式。一種是用座標,另外一種則是使用邏輯向量。在矩陣和陣列上,還是可以用這兩種方式從矩陣或陣列中挑選一部分的值。

關卡 16

我們可以利用中括號[來取出矩陣或陣列中的值。舉例來說,如果x是矩陣,x[2,1]會取出x第二列第一行的值。陣列的狀況,也是類似的。x[1,2,3]會拿出第一個維度的第1個方向,第二個維度的第2個方向,第三個維度的第3個方向。這就是以座標的方式做挑選。

關卡 17

我們來測測看同學有沒有理解剛剛的說明。請問目前這個相當於array(1:18,c(3,3,2))的x,x[1,2,2]的值會是多少呢?

13

關卡 18

在R中,矩陣的資料順序是:c(x[1,1],x[2,1],x[3,1],...)。也就是說,如果我們下指令:matrix(1:18,3,6),則x[1,1]會是1,x[2,1]會是2,x[3,1]是3,x[1,2]是4,以此類推。高維度的矩陣也是類似。我們來試試看,目前的x應該是一個高維陣列,而請問大家:x[3,1,1]的值會是多少呢?

3

關卡 19

如果給一個陣列,我們只是要取一個方向的向量,也可以適當的使用[。先讓我們再印出x一次。

x

關卡 20

我們可以看到,螢幕上的,,2之後,顯示出一個3乘3的矩陣。這個矩陣,就是我們沿著第三個維度的第2個方向,所切面出來的矩陣。這也就是螢幕上的,,2的意思。而在R中,我們可以利用x[,,2]取得這個矩陣。請同學取得這個矩陣,並把它存到變數x2中。

x2 <- x[,,2]

關卡 21

我們先印出x2來看看。

x2

關卡 22

如果我們要選出第一列,只要輸入x2[1,]就可以了。同學可以看一下x2的輸出,在第一列的左方有[1,]的圖示,就是提示使可以用x2[1,]來選第一列。我們請同學試試看選出x2[1,]

x2[1,]

關卡 23

如果我們要直接從x中呼叫,,2之下的矩陣裡面的第一列,我們可以使用:x[1,,2]。這裡的矩陣,在左方有[1,][2,][3,],就代表要取出對應的列,我們要分別使用:x[1,]x[2,]x[3,]這樣的符號。由於目前x是一個三個維度的陣列,所以最後還要加上上面的:,,2,所以就變成:x[1,,2]。請同學試試看。

x[1,,2]

關卡 24

回答下一個問題之前,我們先印出x

x

關卡 25

請問x[,3,1]的值會是什麼呢?

c(7,8,9)

關卡 26

我們也可以利用中括號[]搭配邏輯向量取出矩陣或陣列中的值。R會根據邏輯向量在[]中的位置,選擇該維度,只挑出該邏輯向量為TRUE的座標。

關卡 27

請問同學,依照上述的說明,x[c(TRUE,TRUE,FALSE),,]的效果就等同於以下哪一個程式碼?

x[1:2,,]

關卡 28

接著我們介紹如何修改矩陣和陣列的元素。

關卡 29

如果我們要修改矩陣的值,就是搭配[]<-就可以了。首先,如果我們要更改x[1,1,1]的值為2,我們可以使用x[1,1,1]<-2。請同學試試看。

x[1,1,1] <- 2

關卡 30

我們可以印出x[1,1,1](不是印出全部的x)看看值有沒有被成功地更動。

x[1,1,1]

關卡 31

如果我們想要更改一排的資料,要領也是類似。舉例來說,x[,1,1]<-3可以把,,1底下矩陣的[,1]那行改成c(3,3,3)向量。這樣的運算也是向量式,R會自動重複3直到把x[,1,1]填滿。

x[,1,1] <- 3

關卡 32

我們可以把x[,,1]印出來看看更改後的效果。

x[,,1]

關卡 33

希望這樣的說明可以讓同學了解R如何存取矩陣和陣列的資料。

關卡 34

接下來我們來說明陣列和矩陣的向量式運算。我們先複習一下在向量的狀況下,R的向量式運算。

關卡 35

請問c(1,2,3)+c(2,4,6)的結果會是?

c(3,6,9)

關卡 36

請問c(1,2,3)+2的結果會是?

c(3,4,5)

關卡 37

我們先說明矩陣的向量式運算。它和向量的運算也非常類似,也是會自動比對位置,並且在相同運算的位置上做運算。請同學先印出x[,,1]

x[,,1]

關卡 38

再請同學印出x[,,2]

x[,,2]

關卡 39

再請同學試試看:x[,,1]+x[,,2]

x[,,1] + x[,,2]

關卡 40

這裡相加的結果,就是把x[,,1]x[,,2]對應的位置做相加。同學可以比對前面的輸出結果。

關卡 41

而當維度不相同的時候,R會自動重複比較短的那邊。當一邊是矩陣,一邊是單值的時候,很簡單。請同學試試看:x[,,1]+1

x[,,1] + 1

關卡 42

但是當兩邊都是矩陣或陣列的時候,長度不一樣時,就不容易了。這裡我們先複習一下,陣列的值得排列方式。如果y<-array(1:8,c(2,2,2)),那y[1,1,1]是1、y[2,1,1]是2。R會把1:8的值依照第一個維度(列)、第二個維度(行)到第三個維度,一個一個填入。以此類推,所以y[1,2,1]的值會是?

3

關卡 43

當陣列和向量相加的時候,R會一直重複向量,到它的長度和陣列需要的長度相符合。舉例來說,一個3乘2乘2的陣列需要12個值。接著,R會把這個向量轉換成和陣列一樣的維度,然後在對應的位置做相加。請同學試試看:matrix(1:4,2,2)+1:2

matrix(1:4, 2, 2) + 1:2

關卡 44

另外我們提一個重要的事情,是關於上一個課程中有提到屬性的概念。現在的變數x也是個R物件,也有屬性。請同學查一下它的屬性。

attributes(x)

關卡 45

同學有沒有看到名稱叫做dim的一個整數向量呢?這個dim屬性非常重要,所有的矩陣和陣列,就是一般的向量加上dim這個屬性而已。R也提供了dim這個函數讓使用者可以存取dim屬性。舉例來說,dim(x)會印出x的維度。而我們可以透過dim(x)<-NULL來移除x的dim屬性。請同學試試看dim(x)<-NULL

dim(x) <- NULL

關卡 46

請同學印出x來看看。

x

關卡 47

是不是就變成向量了呢?

關卡 48

R在陣列和向量的向量式運算,也可以回到兩個向量的運算。差別只是在它們多了維度的屬性,所以當維度差異太大的時候,R會認為向量式運算無效。所以一般來說,我們只會拿單值或向量去和陣列做運算,或是維度相同的矩陣或陣列做運算。

關卡 49

基本上,+-*/都會使用向量式運算。

關卡 50

在R中,使用cbindrbind則可以合併兩個矩陣。舉例來說,cbind(matrix(1:4,2,2),matrix(1:4,2,2))會把兩個矩陣的行合併,所以會變成2乘4的矩陣。請同學試試看。

cbind(matrix(1:4, 2, 2), matrix(1:4, 2, 2))

關卡 51

使用rbind(matrix(1:4,2,2),matrix(1:4,2,2))則會把兩個矩陣的列合併,所以會變成4乘2的矩陣。請同學試試看。

rbind(matrix(1:4, 2, 2), matrix(1:4, 2, 2))

關卡 52

最後,我們介紹一些和線性代數相關的運算子和函數。如果對線性代數不熟悉的同學,可以在看過反矩陣之後跳過(使用skip())後面的課程。

關卡 53

在R中兩個矩陣要作矩陣乘法,就是使用%*%這個運算符號。請同學試試看:matrix(1:6,2,3)%*%matrix(3:8,3,2)

matrix(1:6, 2, 3) %*% matrix(3:8, 3, 2)

關卡 54

當兩個向量使用%*%做運算,會得到他們的內積。舉例來說:1:3%*%1:3會得到什麼呢?請同學試試看。

1:3 %*% 1:3

關卡 55

當矩陣乘向量的時候,單一向量會被視為行向量,最後會得到矩陣每個列和單一向量的內積,形成一個新的行向量。舉例來說:matrix(1:9,3,3)%*%1:3會得到什麼呢?請同學試試看。

matrix(1:9,3,3) %*% 1:3

關卡 56

R上的線性代數運算的底層是透過BLAS等函式庫做運算的,所以效能遠勝過我們自己用C寫的線性代數運算。另外R預設是使用比較被廣泛驗證過正確性的BLAS庫,而不是效能比較快,但是還比較年輕的BLAS函式庫,例如OpenBLAS。這是因為,RcoreTeam認為正確性比較重要,所以目前是採用比較舊,但是也比較可靠的BLAS庫。

關卡 57

在R中,要對一個矩陣做轉置,可以用t這個函數。請同學先印出matrix(1:4,2,2)

matrix(1:4, 2, 2)

關卡 58

再印出:t(matrix(1:4,2,2))

t(matrix(1:4, 2, 2))

關卡 59

我們可以用diag快速建構對角化的矩陣。舉例來說:diag(1,3)會建立出3乘3的單位矩陣。請試試看。

diag(1, 3)

關卡 60

如果已知A%*%x=b,給定Ab我們可以用solve解出x。舉例來說,若Amatrix(1:4,2,2)、b是c(3,8)solve(A,b)就會給出x。變數Ab已經存在了,所以請同學直接解出x並且把它存到x之中。

x <- solve(A, b)

關卡 61

我們來驗算一下,A%*%x看看是不是真的是b

A %*% x

關卡 62

如果直接使用solve(A),則R會算出A的反矩陣。請同學試試看。

solve(A)

關卡 63

使用eigen(A)則可以解出A的特徵值(eigenvalues)和特徵向量(eigenvectors)。我們可以試試看:A%*%eigen(A)$vectors[,1]/eigen(A)$vectors[,1]是不是相除後,會拿到一模一樣的值呢?這裡的$用法,等到下一個課程,我們會有詳細的說明。

A %*% eigen(A)$vectors[,1] / eigen(A)$vectors[,1]

關卡 64

最後,我們請同學用所學到的技巧,做一個簡短的練習。請同學在完成之後存檔,並輸入submit()來檢查結果是否符合預期。如果同學在檔案中看到亂碼,請使用Rstudio左上角的File->ReopenWithEncoding…->選取:UTF-8

X <- cbind(x1 = 1, x2 = 1:10, x3 = sin(1:10))
beta <- c(0.5, -1, 4.3)

y <- X %*% beta

#' epsilon 是一個隨機產生的雜訊向量
epsilon <- c(-1.24462014500259, 0.146172987456978, 1.56426869006839, -0.856920339050681, 
    -1.15277300953772, 0.717919832604741, -0.270623615316431, -1.66281578024014, 
    -1.15557078461633, -0.730253254897595)
#' 我們讓y 參雜了雜訊
y <- y + epsilon

beta.hat <- solve(t(X)%*%X, t(X)%*%y)

#' 同學可以比較一下beta.hat和beta,體驗看看迴歸分析的方法,是不是真的有道理。