多維數(shù)組和矩陣

數(shù)組（array）和矩陣(matrix)

數(shù)組（array）可以看成是帶多個(gè)下標(biāo)的類型相同的元素的集合，常用的是數(shù)值型的數(shù)組如矩陣，也可以有其它類型（如字符型、邏輯型、復(fù)型數(shù)組）。S可以很容易地生成和處理數(shù)組，特別是矩陣（二維數(shù)組）。

數(shù)組有一個(gè)特征屬性叫做維數(shù)向量（dim屬性），維數(shù)向量是一個(gè)元素取正整數(shù)值的向量，其長(zhǎng)度是數(shù)組的維數(shù)，比如維數(shù)向量有兩個(gè)元素時(shí)數(shù)組為二維數(shù)組（矩陣）。維數(shù)向量的每一個(gè)元素指定了該下標(biāo)的上界，下標(biāo)的下界總為1。

一組值只有定義了維數(shù)向量（dim屬性）后才能被看作是數(shù)組。比如：
> z <- 1:1500
> dim(z) <- c(3, 5, 100)

這時(shí)z已經(jīng)成為了一個(gè)維數(shù)向量為c(3,5,100)的三維數(shù)組。也可以把向量定義為一維數(shù)組，例如：

> dim(z) <- 1500

數(shù)組元素的排列次序缺省情況下是采用FORTRAN的數(shù)組元素次序（按列次序），即第一下標(biāo)變化最快，最后下標(biāo)變化最慢，對(duì)于矩陣（二維數(shù)組）則是按列存放。例如，假設(shè)數(shù)組a的元素為1:24，維數(shù)向量為c(2,3,4)，則各元素次序?yàn)閍[1,1,1], a[2,1,1], a[1,2,1], a[2,2,1], a[1,3,1], ..., a[2,3,4]。

用函數(shù)array()或matrix()可以更直觀地定義數(shù)組。array()函數(shù)的完全使用為array(x, dim=length(x), dimnames=NULL)，其中x是第一自變量，應(yīng)該是一個(gè)向量，表示數(shù)組的元素值組成的向量。dim參數(shù)可省，省略時(shí)作為一維數(shù)組（但不同于向量）。dimnames屬性可以省略，不省略時(shí)是一個(gè)長(zhǎng)度與維數(shù)相同的列表（list，見后面），列表的每個(gè)成員為一維的名字。例如上面的z可以這樣定義：

> z <- array(1:1500, dim=c(3,5,100))
	函數(shù)matrix()用來定義最常用的一種數(shù)組：二維數(shù)組，即矩陣。其完全格式為
matrix(data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL)

其中第一自變量data為數(shù)組的數(shù)據(jù)向量（缺省值為缺失值NA），nrow為行數(shù)，ncol為列數(shù)，byrow表示數(shù)據(jù)填入矩陣時(shí)按行次序還是列次序，一定注意缺省情況下按列次序，這與我們寫矩陣的習(xí)慣是不同的。dimnames缺省是空值，否則是一個(gè)長(zhǎng)度為2的列表，列表第一個(gè)成員是長(zhǎng)度與行數(shù)相等的字符型向量，表示每行的標(biāo)簽，列表第二個(gè)成員是長(zhǎng)度與列數(shù)相同的字符型向量，表示每列的標(biāo)簽。例如，定義一個(gè)3行4列，由1:12按行次序排列的矩陣，可以用：

> b <- matrix(1:12, ncol=4, byrow=T)
> b
     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    5    6    7    8
[3,]    9   10   11   12

注意在有數(shù)據(jù)的情況下只需指定行數(shù)或列數(shù)之一。指定的數(shù)據(jù)個(gè)數(shù)允許少于所需的數(shù)據(jù)個(gè)數(shù)，這時(shí)循環(huán)使用提供的數(shù)據(jù)。例如：

> b <- matrix(0, nrow=3, ncol=4)

數(shù)組下標(biāo)

要訪問數(shù)組的某個(gè)元素，只要象上面那樣寫出數(shù)組名和方括號(hào)內(nèi)用逗號(hào)分開的下標(biāo)即可，如a[2,1,2]。

更進(jìn)一步我們還可以在每一個(gè)下標(biāo)位置寫一個(gè)下標(biāo)向量，表示對(duì)這一維取出所有指定下標(biāo)的元素，如a[1, 2:3, 2:3]取出所有第一下標(biāo)為1，第二下標(biāo)為2或3，第三下標(biāo)為2或3的元素。注意因?yàn)榈谝痪S只有一個(gè)下標(biāo)所以退化了，得到的是一個(gè)維數(shù)向量為c(2,2)的數(shù)組。

另外，如果略寫某一維的下標(biāo)，則表示該維全選。例如，a[1, , ]取出所有第一下標(biāo)為1 的元素，得到一個(gè)形狀為c(3,4)的數(shù)組。a[ , 2, ]取出所有第二下標(biāo)為2的元素得到一個(gè)形狀為c(2,4)的數(shù)組。a[1,1, ]則只能得到一個(gè)長(zhǎng)度為4的向量，不再是數(shù)組（dim(a[1,1, ]) 值為NULL）。a[ , , ]或a[]都表示整個(gè)數(shù)組。比如

a[] <- 0

還有一種特殊下標(biāo)是對(duì)于數(shù)組只用一個(gè)下標(biāo)向量（是向量，不是數(shù)組），比如a[3:4] ，這時(shí)忽略數(shù)組的維數(shù)信息，把下標(biāo)表達(dá)式看作是對(duì)數(shù)組的數(shù)據(jù)向量取子集。

不規(guī)則數(shù)組下標(biāo)

在S中甚至可以把數(shù)組中的任意位置的元素作為一組訪問，其方法是用一個(gè)矩陣作為數(shù)組的下標(biāo)，矩陣的每一行是一個(gè)元素的下標(biāo)，數(shù)組有幾維下標(biāo)矩陣的每一行就有幾列。例如，我們要把上面的形狀為c(2,3,4)的數(shù)組a的第[1,1,1], [2,2,3], [1,3,4], [2,1,4]號(hào)共四個(gè)元素作為一個(gè)整體訪問，先定義一個(gè)包含這些下標(biāo)作為行的二維數(shù)組:

> b <- matrix(c(1,1,1, 2,2,3, 1,3,4, 2,1,4), ncol=3, byrow=T)
> b
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    2    2    3
[3,]    1    3    4
[4,]    2    1    4
> a[b]
[1]  1 16 23 20

> a[b] <- c(101, 102, 103, 104)
> a
或
> a[b] <- 0
> a

數(shù)組四則運(yùn)算

數(shù)組可以進(jìn)行四則運(yùn)算（+，－, *, /，^），解釋為數(shù)組對(duì)應(yīng)元素的四則運(yùn)算，參加運(yùn)算的數(shù)組一般應(yīng)該是相同形狀的（dim屬性完全相同）。例如，假設(shè)A, B, C是三個(gè)形狀相同的數(shù)組，則

> D <- C + 2*A/B

計(jì)算得到的結(jié)果是A的每一個(gè)元素除以B的對(duì)應(yīng)元素加上C的對(duì)應(yīng)元素乘以2得到相同形狀的數(shù)組。四則運(yùn)算遵循通常的優(yōu)先級(jí)規(guī)則。

形狀不一致的向量和數(shù)組也可以進(jìn)行四則運(yùn)算，一般的規(guī)則是數(shù)組的數(shù)據(jù)向量對(duì)應(yīng)元素進(jìn)行運(yùn)算，把短的循環(huán)使用來與長(zhǎng)的匹配，并盡可能保留共同的數(shù)組屬性。例如：

> x1 <- c(100, 200)
> x2 <- 1:6
> x1+x2
[1] 101 202 103 204 105 206
> x3 <- matrix(1:6, nrow=3)
> x3
     [,1] [,2]
[1,]    1    4
[2,]    2    5
[3,]    3    6
> x1+x3
     [,1] [,2]
[1,]  101  204
[2,]  202  105
[3,]  103  206

矩陣運(yùn)算

矩陣是二維數(shù)組，但因?yàn)槠鋺?yīng)用廣泛所以對(duì)它定義了一些特殊的運(yùn)算和操作。

函數(shù)t(A)返回矩陣A的轉(zhuǎn)置。nrow(A)為矩陣A的行數(shù)，ncol(A)為矩陣A的列數(shù)。

矩陣之間進(jìn)行普通的加減乘除四則運(yùn)算仍遵從一般的數(shù)組四則運(yùn)算規(guī)則，即數(shù)組的對(duì)應(yīng)元素之間進(jìn)行運(yùn)算，所以注意A*B不是矩陣乘法而是矩陣對(duì)應(yīng)元素相乘。

要進(jìn)行矩陣乘法，使用運(yùn)算符%*%，A%*%B表示矩陣A乘以矩陣B（當(dāng)然要求A的列數(shù)等于B的行數(shù)）。例如：

> A <- matrix(1:12, nrow=4, ncol=3, byrow=T)
> B <- matrix(c(1,0), nrow=3, ncol=2, byrow=T)
> A
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[3,]    7    8    9
[4,]   10   11   12
> B
     [,1] [,2]
[1,]    1    0
[2,]    1    0
[3,]    1    0
> A %*% B
     [,1] [,2]
[1,]    6    0
[2,]   15    0
[3,]   24    0
[4,]   33    0
>

另外，向量用在矩陣乘法中可以作為行向量看待也可以作為列向量看待，這要看哪一種觀點(diǎn)能夠進(jìn)行矩陣乘法運(yùn)算。例如，設(shè)x是一個(gè)長(zhǎng)度為n的向量，A是一個(gè) 矩陣，則“x %*% A %*% x”表示二次型 。但是，有時(shí)向量在矩陣乘法中的地位并不清楚，比如“x %*% x”就既可能表示內(nèi)積 也可能表示 陣 。因?yàn)榍罢咻^常用，所以S選擇表示前者，但內(nèi)積最好還是用crossprod(x)來計(jì)算。要表示 ，可以用“cbind(x) %*% x”或“x %*% rbind(x) ”。

函數(shù)crossprod(X, Y)表示一般的交叉乘積（內(nèi)積） ，即X的每一列與Y的每一列的內(nèi)積組成的矩陣。如果X和Y都是向量則是一般的內(nèi)積。只寫一個(gè)參數(shù)X的crossprod(X)計(jì)算X自身的內(nèi)積 。

其它矩陣運(yùn)算還有solve(A,b)解線性方程組 ，solve(A)求方陣A的逆矩陣，svd()計(jì)算奇異值分解，qr()計(jì)算QR分解，eigen()計(jì)算特征向量和特征值。詳見隨機(jī)幫助，例如：

> ?qr

函數(shù)diag()的作用依賴于其自變量。diag(vector)返回以自變量（向量）為主對(duì)角元素的對(duì)角矩陣。diag(matrix)返回由矩陣的主對(duì)角元素組成的向量。diag(k)（k為標(biāo)量）返回k階單位陣。

矩陣合并與拉直

函數(shù)cbind()把其自變量橫向拼成一個(gè)大矩陣，rbind()把其自變量縱向拼成一個(gè)大矩陣。cbind()的自變量是矩陣或者看作列向量的向量，自變量的高度應(yīng)該相等（對(duì)于向量，高度即長(zhǎng)度，對(duì)于矩陣，高度即行數(shù)）。rbind的自變量是矩陣或看作行向量的向量，自變量的寬度應(yīng)該相等（對(duì)于向量，寬度即長(zhǎng)度，對(duì)于矩陣，寬度即列數(shù)）。如果參與合并的自變量比其它自變量短則循環(huán)補(bǔ)足后合并。例如：

> x1 <- rbind(c(1,2), c(3,4))
> x1
     [,1] [,2]
[1,]    1    2
[2,]    3    4
> x2 <- 10+x1
> x3 <- cbind(x1, x2)
> x3
     [,1] [,2] [,3] [,4]
[1,]    1    2   11   12
[2,]    3    4   13   14
> x4 <- rbind(x1, x2)
> x4
     [,1] [,2]
[1,]    1    2
[2,]    3    4
[3,]   11   12
[4,]   13   14
> cbind(1, x1)
     [,1] [,2] [,3]
[1,]    1    1    2
[2,]    1    3    4

因?yàn)閏bind()和rbind()的結(jié)果總是矩陣類型（有dim屬性且為二維），所以可以用它們把向量表示為 矩陣（用cbind(x)）或 矩陣（用rbind(x)）。

設(shè)a是一個(gè)數(shù)組，要把它轉(zhuǎn)化為向量（去掉dim和dimnames屬性），只要用函數(shù)as.vector(a) 返回值就可以了（注意函數(shù)只能通過函數(shù)值返回結(jié)果而不允許修改它的自變量，比如t(X)返回X的轉(zhuǎn)置矩陣而X本身并未改變）。另一種由數(shù)組得到其數(shù)據(jù)向量的簡(jiǎn)單辦法是使用函數(shù)c() ，例如c(a)的結(jié)果是a的數(shù)據(jù)向量。這樣的另一個(gè)好處是c()允許多個(gè)自變量，它可以把多個(gè)自變量都看成數(shù)據(jù)向量連接起來。例如，設(shè)A和B是兩個(gè)矩陣，則c(A,B)表示把A按列次序拉直為向量并與把B按列次序拉直為向量的結(jié)果連接起來。一定注意拉直時(shí)是按列次序拉直的。

數(shù)組的維名字

數(shù)組可以有一個(gè)屬性dimnames保存各維的各個(gè)下標(biāo)的名字，缺省時(shí)為NULL（即無此屬性）。我們以矩陣為例，它有兩個(gè)維：行維和列維。比如，設(shè)x為2行3列矩陣，它的行維可以定義一個(gè)長(zhǎng)度為2的字符向量作為每行的名字，它的列維可以定義一個(gè)長(zhǎng)度為3的向量作為每列的名字，屬性dimnames是一個(gè)列表，列表的每個(gè)成員是一個(gè)維名字的字符向量或NULL。例如：

> x <- matrix(1:6, ncol=2,
+ dimnames=list(c("one", "two", "three"), c("First", "Second")),
+ byrow=T)
> x
      First Second
one       1      2
two       3      4
three     5      6
我們也可以先定義矩陣x然后再為dimnames(x)賦值。
對(duì)于矩陣，我們還可以使用屬性rownames和colnames來訪問行名和列名。如：
> x <- matrix(1:6, ncol=2,byrow=T)
> colnames(x) <- c("First", "Second")
> rownames(x) <- c("one", "two", "three")
在定義了數(shù)組的維名后我們對(duì)這一維的下標(biāo)就可以用它的名字來訪問，例如：
> x[c("one","three"),]
      First Second
one       1      2
three     5      6

數(shù)組的外積

兩個(gè)數(shù)組a和b的外積是由a的每一個(gè)元素與b的每一個(gè)元素搭配在一起相乘得到一個(gè)新元素，這樣得到一個(gè)維數(shù)向量等于a的維數(shù)向量與b的維數(shù)向量連起來的數(shù)組，即若d為a和b的外積，則dim(d)=c(dim(a), dim(b))。

> d <- a %o% b

> d <- outer(a, b, '*')

注意outer(a, b, f)是一個(gè)一般性的外積函數(shù)，它可以把a(bǔ)的任一個(gè)元素與b的任意一個(gè)元素搭配起來作為f的自變量計(jì)算得到新的元素值，外積是兩個(gè)元素相乘的情況。函數(shù)當(dāng)然也可以是加、減、除，或其它一般函數(shù)。當(dāng)函數(shù)為乘積時(shí)可以省略不寫。

例如，我們希望計(jì)算函數(shù) 在一個(gè)x和y的網(wǎng)格上的值用來繪制三維曲面圖，可以用如下方法生成網(wǎng)格及函數(shù)值：

> x <- seq(-2, 2, length=20)
> y <- seq(-pi, pi, length=20)
> f <- function(x, y) cos(y)/(1+x^2)
> z <- outer(x, y, f)

下面考慮一個(gè)有意思的例子。我們考慮簡(jiǎn)單的2×2矩陣 ，其元素均在0,1,...,9中取值。假設(shè)四個(gè)元素 a, b, c, d都是相互獨(dú)立的服從離散均勻分布的隨機(jī)變量，我們?cè)O(shè)法求矩陣行列式 ad-bc的分布。首先，隨機(jī)變量 ad和 bc同分布，它的取值由以下外積矩陣給出，每一個(gè)取值的概率均為1/100:

> d <- outer(0:9, 0:9)

這個(gè)語句產(chǎn)生一個(gè) 的外積矩陣。為了計(jì)算 ad的100個(gè)值（有重復(fù)）與 bc的100個(gè)值相減得到的10000個(gè)結(jié)果，可以使用如下外積函數(shù)：

> d2 <- outer(d, d, "-")

這樣得到一個(gè)維數(shù)向量為c(2,2,2,2)的四維數(shù)組，每一個(gè)元素為行列式的一個(gè)可能取值，概率為萬分之一。因?yàn)檫@些取值中有很多重復(fù)，我們可以用一個(gè)table()函數(shù)來計(jì)算每一個(gè)值的出現(xiàn)次數(shù)（頻數(shù)）：

> fr <- table(d2)

得到的結(jié)果是一個(gè)帶有元素名的向量fr，fr的元素名為d2的一個(gè)取值，fr的元素值為d2該取值出現(xiàn)的頻數(shù)，比如fr[1]的元素名為－81，值為19，表示值－81在數(shù)組d2中出現(xiàn)了19次。通過計(jì)算length(fr)可以知道共有163個(gè)不同值。還可以把這些值繪制一個(gè)頻數(shù)分布圖（除以10000 則為實(shí)際概率）：

> plot(as.numeric(names(fr)), fr, type="h",
+  xlab="行列式", ylab="頻數(shù)")

其中as.numeric()把向量fr中的元素名又轉(zhuǎn)換成了數(shù)值型，用來作為作圖的橫軸坐標(biāo)，fr 中的元素值即頻數(shù)作為縱軸，type="h"表示是畫垂線型圖。

數(shù)組的廣義轉(zhuǎn)置

> a <- array(1:24, dim=c(2,3,4))
> b <- aperm(a, c(2, 3, 1))

結(jié)果a的第2維變成了b的第1維，a的第3維變成了b的第2維，a的第1維變成了b的第3維。這時(shí)有a[i1,i2,i3]≡b[i2,i3,i1]。注意c(i1,i2,i3)[2,3,1]=c(i2,i3,i1)。一般地，若b <- aperm(a, p)，i是數(shù)組a的一個(gè)下標(biāo)向量，則a[rbind(i)]≡b[rbind(i[p])]，即a的一個(gè)元素下標(biāo)經(jīng)過p變換后成為b的對(duì)應(yīng)元素的下標(biāo)。

對(duì)于矩陣a，aperm(a, c(2,1))恰好是矩陣轉(zhuǎn)置。對(duì)于矩陣轉(zhuǎn)置可以簡(jiǎn)單地用t(a)表示。

apply函數(shù)

對(duì)于向量，我們有sum、mean等函數(shù)對(duì)其進(jìn)行計(jì)算。對(duì)于數(shù)組，如果我們想對(duì)其中一維（或若干維）進(jìn)行某種計(jì)算，可以用apply函數(shù)。其一般形式為：

apply(X, MARGIN, FUN, ...)

其中X為一個(gè)數(shù)組，MARGIN是固定哪些維不變，F(xiàn)UN是用來計(jì)算的函數(shù)。例如，設(shè)a是 矩陣，則apply(a, 1, sum)的意義是對(duì)a的各行求和（保留第一維即第一個(gè)下標(biāo)不變），結(jié)果是一個(gè)長(zhǎng)度為3的向量（與第一維長(zhǎng)度相同），而apply(a, 2, sum)意義是對(duì)a的各列求和，結(jié)果是一個(gè)長(zhǎng)度為4的向量（與第二維長(zhǎng)度相同）。

如果函數(shù)FUN的結(jié)果是一個(gè)標(biāo)量，MARGIN只有一個(gè)元素，則apply的結(jié)果是一個(gè)向量，其長(zhǎng)度等于MARGIN指定維的長(zhǎng)度，相當(dāng)于固定MARGIN指定的那一維的每一個(gè)值而把其它維取出作為子數(shù)組或向量送入FUN中進(jìn)行運(yùn)算。如果MARGIN指定了多個(gè)維，則結(jié)果是一個(gè)維數(shù)向量等于dim(X)[MARGIN]的數(shù)組。如果函數(shù)FUN的結(jié)果是一個(gè)長(zhǎng)度為N的向量，則結(jié)果是一個(gè)維數(shù)向量等于c(N, dim(X)[MARGIN])的數(shù)組，注意這時(shí)不論是對(duì)哪一維計(jì)算，結(jié)果都放在了第一維。所以，比如我們要把4×3矩陣a的3列分別排序，只要用apply(a, 2, sort)，這樣對(duì)每一列排序得到一個(gè)長(zhǎng)度為4的向量，用第一維來引用，結(jié)果的維向量為c(N, dim(a)[2])=c(4,3) ，保留了列維，恰好得到所需結(jié)果，運(yùn)行如下例：

> a <- cbind(c(4,9,1), c(3,7,2))
> a
     [,1] [,2]
[1,]    4    3
[2,]    9    7
[3,]    1    2
> apply(a, 2, sort)
     [,1] [,2]
[1,]    1    2
[2,]    4    3
[3,]    9    7
> 

但是，如果要對(duì)行排序，則apply(a, 1, sort)把a(bǔ)的每一行3個(gè)元素排序后的結(jié)果用第一維來引用，結(jié)果的維向量為c(N, dim(a)[1])=c(3, 4)，把原來的列變成了行，所以t(apply(a,1,sort)) 才是對(duì)a的每一行排序的結(jié)果。如：

> apply(a, 1, sort)
     [,1] [,2] [,3]
[1,]    3    7    1
[2,]    4    9    2
> t(apply(a,1,sort))
     [,1] [,2]
[1,]    3    4
[2,]    7    9
[3,]    1    2

上面我們只用了矩陣（二維數(shù)組）作為例子講解apply的用法。實(shí)際上，apply可以用于任意維數(shù)的數(shù)組，函數(shù)FUN也可以是任意可以接收一個(gè)向量或數(shù)組作為第一自變量的函數(shù)。比如，設(shè)x是一個(gè)維數(shù)向量為c(2,3,4,5)的數(shù)組，則apply(x, c(1,3), sum)可以產(chǎn)生一個(gè)2行4 列的矩陣，其每一元素是x中固定第1維和第3維下標(biāo)取出子數(shù)組求和的結(jié)果。