R中的多维数组和矩阵

1.生成数组或矩阵

数组可以看成是带多个下标的类型相同的元素的集合，常用的数值型的数组如矩阵，也可以有其他类型（字符型、逻辑型、复数型）。R可以很容易地生成和处理数组，特别是矩阵（二维数组）。

数组有一个特殊的属性叫做维数向量，维数向量是一个元素取正整数值得向量，其长度是数组的维数，比如维数向量有两个元素时数组为二维数组（矩阵）。维数向量的每一个元素指定了该下标的上界，下标的下界总为1.

（（）将向量定义成数组

向量只有定义了维数向量后才能被看作为数组，比如：

> z<-1:12

> dim(z)<-c(3,4)

> z

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

[1,]   1    4    7  10

[2,]   2    5    8  11

[3,]   3    6    9  12

注意：矩阵的元素时按列存放的。也可以吧向量定义为一维数组，例如：

> dim(z)<-12

> z

 [1] 1  2  3 4  5  6 7  8  9 10 11 12

 

（2）用array()函数构造多维数组

array(data,dim,dimnames)

其中data是一个向量数据，dim()是数组各维的长度，缺省时为原向量的长度。dimnames是数组维的名字，缺省为空。如：

> x<-array(1:20,dim=c(4,5))

> x

    [,1] [,2] [,3] [,4] [,5]

[1,]   1    5    9  13   17

[2,]   2    6   10  14   18

[3,]   3    7   11  15   19

[4,]   4    8   12  16   20

另一种方式：

> z<-array(0,dim=c(3,4,2))

> z

, , 1

 

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

[1,]   0    0    0   0

[2,]   0    0    0   0

[3,]   0    0    0   0

 

, , 2

 

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

[1,]   0    0    0   0

[2,]   0    0    0   0

[3,]   0    0    0   0

如此定义了一个3*4*2的三维数组，元素均为0。这种方法常用来对数组作初始化。

 

（3）用matrix()函数构造矩阵

matrix(data=NA,nrow=1,ncol=1,byrow=FALSE,dimnames=NULL)

其中data是一个向量数据，nrow是矩阵的行数，ncol为列数，nrow与ncol的乘积是矩阵元素个数，byrow项控制排列元素是否按行进行，当为T时按行放置，当为F时按列放置，dimnames给定行和列的名称，缺省时为空。如：

> A<-matrix(1:15,nrow=3,ncol=5,byrow=TRUE)

> A

    [,1] [,2] [,3] [,4] [,5]

[1,]   1    2    3   4    5

[2,]   6    7    8   9   10

[3,]  11   12   13  14   15

以下两种格式与前面的格式是等价的：

> A<-matrix(1:15,nrow=3,byrow=TRUE)

> A<-matrix(1:15, ncol=5,byrow=TRUE)

如果把语句中的byrow=TRUE去掉或者改为byrow=FALSE，则数据按列放置：

> A<-matrix(1:15,nrow=3,ncol=5)

> A

    [,1] [,2] [,3] [,4] [,5]

[1,]   1    4    7  10   13

[2,]   2    5    8  11   14

[3,]   3    6    9  12   15

 

（4）用合并命令构造矩阵

用rbind()、cbind()将两个或两个以上的向量或矩阵合并起来，rbind()表示按行合并，cbind()表示按列合并。例：

> x<-c(1:5)

> y<-c(6:10)

> rbind(x,y)

 [,1] [,2] [,3] [,4] [,5]

x   1    2    3   4    5

y   6    7    8   9   10

> cbind(x,y)

    x  y

[1,] 1 6

[2,] 2 7

[3,] 3 8

[4,] 4 9

[5,] 5 10

 

2.数组下标

数组和向量一样，可以对数组中的某些元素进行访问，或进行运算。

（1）数组下标

要访问数组的某个元素，只要写出数组名和方括号内的用逗号分开的下标即可。如：

> a<-1:24

> dim(a)<-c(2,3,4)

> a[2,1,2]

[1] 8

更进一步还可以在每一个下标位置写一个下标向量，表示这一维取出的所有指定下标的元素，如：

> a[1,2:3,2:3]

    [,1] [,2]

[1,]   9   15

[2,]  11   17

取出了第一下标为1，第二下标为2或3.第三下标为2或3的元素。

注意，因为第一维只有一个下标，所以退化成一个2*2的数组。

另外，如果略写某一维的下标，则表示该维全选，如：

> a[1,,]

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

[1,]   1    7   13  19

[2,]   3    9   15  21

[3,]   5   11   17  23

取出所有第一下标为1的元素，得到一个3*4的数组。

> a[,2,]

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

[1,]   3    9   15  21

[2,]   4   10   16  22

取出所有第二下标为2的元素得到一个2*4的数组。

> a[1,1,]

[1] 1  7 13 19

得到一个长度为4的向量，不再是数组。

a[,,]或a[]都表示整个数组，如：

> a[]<-0

可以在不改变数组维数的条件下把元素都赋为0.

 

（2）不规则的数组下标

在R中可以把数组中的任意位置的元素作为数组访问，方法是用一个二维数组作为数组的下标，二维数组的每一行是一个元素的下标，列数为数组的维数。

> z<-1:24

> dim(z)<-c(2,3,4)

> 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

> z[b]

[1] 1 16 23 20

表示把2*3*4的数组z的第[1,1,1],[2,2,3],[1,3,4],[2,1,4]这四个元素取出。

注意：取出的是一个向量，我们还可以对这几个元素赋值，如：

> z[b]<-c(101,102,103,104)

或

> z[b]<-0

 

3.数组的四则运算

可以对数组之间进行四则运算，这时进行的是数组对应元素的四则运算，参加运算的数组一般应该是相同形状的（dim属性完全相同）。如：

> A=B=matrix(1:12,nrow=3,ncol=4)

> A+B

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

[1,]   2    8   14  20

[2,]   4   10   16  22

[3,]   6   12   18  24

> A-B

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

[1,]   0    0    0   0

[2,]   0    0    0   0

[3,]   0    0    0   0

> C<-A*B;D=A/B;C;D

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

[1,]   1   16   49 100

[2,]   4   25   64 121

[3,]   9   36   81 144

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

[1,]   1    1   1    1

[2,]   1    1    1   1

[3,]   1    1    1   1

数组的加、减运算和数乘运算满足原矩阵运算的性质，但数组的乘、除运算实际上是数组中对应位置的元素作运算。

 

形状不一致的向量（或数组）中的数据进行四则运算，一般的规则是将向量（或数组）中的数据与对应向量（或数组）中的数据进行运算，把段向量（或数组）的数据循环使用，从而可以与长向量（或数组）数据进行匹配，并尽可能保留共同的数组属性。如：

> x1<-c(100,200)

> x2<-1:6

> x1+x2

[1] 101 202 103 204 105 206

> x1+x3

    [,1] [,2]

[1,] 101  204

[2,] 202  105

[3,] 103  206

可以看到，当向量与数组共同运算时，向量按列匹配。当两个数组不匹配时，R会提出警告。如：

> x2<-1:5

> x1+x2

[1] 101 202 103 204 105

警告信息：

In x1 + x2 : 长的对象长度不是短的对象长度的整倍

 

4.矩阵的运算

（1）转置

求矩阵A的转置用函数t()。例如：

> A=matrix(1:12,nrow=3,ncol=4);A

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

[1,]   1    4    7  10

[2,]   2    5    8  11

[3,]   3    6    9  12

> t(A)

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

[1,]   1    2    3

[2,]   4    5    6

[3,]   7    8    9

[4,]  10   11   12

（2）方阵的行列式

det()函数可求方阵行列式的值，例如：

> det(matrix(1:4,ncol=2))

[1] -2

（3）向量的内积

对于n维向量x，可以看成n*1阶矩阵或1*n阶矩阵。若x与y是相同维数的向量，则x%*%y表示x与y作内积，例如：

> x<-1:5;y<-2*1:5

> x%*%y

    [,1]

[1,] 110

函数crossprod()是内积计算函数（表示交叉乘积），crossprod(x,y)加us那向量x与y的内积，即‘t(x)%*%y‘。crossprod(x)表示x与x的内积。

类似地，tcrossprod(x,y)表示‘x%*% t(y) ‘，即x与y的外积，也称为叉积，tcrossprod(x)表示x与x的外积。

（4）向量的外积（叉积）

设x，y是n维向量，则x%*o%y表示x与y作外积。如：

> x%o%y

    [,1] [,2] [,3] [,4] [,5]

[1,]   2    4    6   8   10

[2,]   4    8   12  16   20

[3,]   6   12   18  24   30

[4,]   8   16   24  32   40

[5,]  10   20   30  40   50

函数outer()是外积运算函数，outer(x,y)计算x与y的外积，等价于x%o%y。

函数outer()的一般格式为：outer(x,y,fun)

其中x，y矩阵（或向量），fun是作外积运算函数，缺省值为乘法运算。函数outer()在绘制三维曲面时非常有用。

（5）矩阵的乘法

A%*%B表示通常意义下的两个矩阵的乘积（当然要求矩阵A的列数等于矩阵B的行数），如：

> A=matrix(1:12,nrow=3,ncol=4)

> B=matrix(1:12,nrow=4,ncol=3)

> A%*%B

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

[1,]  70  158  246

[2,]  80  184  288

[3,]  90  210  330

（6）生成对角阵和矩阵取对角运算

若取一个方阵的对角元素，用函数diag()。

对一个向量用函数diag()，将产生以这个向量为对角元素、其余元素全为0的对角矩阵。

对一个正整数k应用函数diag()，将产生k维单位矩阵。

例：

> A=matrix(1:16,nrow=4,ncol=4)

> A

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

[1,]   1    5    9  13

[2,]   2    6   10  14

[3,]   3    7   11  15

[4,]   4    8   12  16

> diag(A)

[1] 1  6 11 16

> diag(3)

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

[1,]   1    0    0

[2,]   0    1    0

[3,]   0    0    1

 

（7）解线性方程组和求矩阵的逆矩阵

阵求逆用函数solve()。

用solve(A,b)运算，结果是可解线性方程组Ax=b的解x。若b为单位矩阵，则结果就是A的逆矩阵。solve()中，若b缺省，系统默认为单位矩阵，则这个过程便为矩阵求逆。

例：

rnorm()为随机生成正态

> A=matrix(rnorm(16),4,4);A

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

[1,] 0.67168810  0.77743512 -2.1351143  1.1196029

[2,] 0.70557911 -0.40344244  0.5866727 1.2809767

[3,] 0.57529587 -0.02789363 -1.0624213-1.2417952

[4,] 0.08381951  0.65116106 0.1662728 -0.9488299

> solve(A)

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

[1,] 0.02554831  0.89633961  0.5293070 0.5475198

[2,] 0.34677139 -0.02034261 -0.5381588 1.0860438

[3,] -0.22913613  0.34001906 -0.2186800  0.4748699

[4,] 0.20008471  0.12480673 -0.3608891-0.1770177

（8）矩阵的特征值和特征向量

矩阵A的谱分析为A=UDU’，其中D为由A的特征值组成的对角矩阵，U的列为A的特征值对应的特征向量，在R中可用函数eigen()得到U和D。

eigen(x,symmetric,only.value=FALSE)

x为矩阵，symmetric项指定矩阵x是否为对称矩阵，若不指定，系统将自动检测x是否为对称矩阵。only.value项如果是TRUE，只有特征值计算并返回，否则返回的特征值和特征向量。

例：

> A=diag(4)+1

> A

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

[1,]   2    1    1   1

[2,]   1    2    1   1

[3,]   1    1   2    1

[4,]   1    1    1   2

> A.e=eigen(A,symmetric=T)

> A.e

$values

[1] 5 1 1 1

 

$vectors

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

[1,] -0.5 0.000000e+00  0.0000000  0.8660254

[2,] -0.5 -6.408849e-17  0.8164966 -0.2886751

[3,] -0.5 -7.071068e-01 -0.4082483-0.2886751

[4,] -0.5 7.071068e-01 -0.4082483 -0.2886751

 

>A.e$vectors%*%diag(A.e$values)%*%t(A.e$vectors)

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

[1,]   2    1    1   1

[2,]   1    2    1   1

[3,]   1    1    2   1

[4,]   1    1    1    2

（9）矩阵的Choleskey分解

对于正定矩阵A，可对其进行Choleskey分解，即A=P’P，其中，P为是上三角矩阵，在R中可用函数chol()来进行Choleskey分解。

正定矩阵是指

例：

> A.c=chol(A)

> A.c

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

[1,] 1.414214 0.7071068 0.7071068 0.7071068

[2,] 0.000000 1.2247449 0.4082483 0.4082483

[3,] 0.000000 0.0000000 1.1547005 0.2886751

[4,] 0.000000 0.0000000 0.0000000 1.1180340

> t(A.c)%*%A.c

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

[1,]   2    1    1   1

[2,]   1    2    1   1

[3,]   1    1    2   1

[4,]   1    1    1   2

（10）矩阵奇异值分解

A为m*n矩阵，A的秩为n，即rank(A)=n，可分解为A=UDV’，其中U’U=V’V=I。在R中可以用函数svd()进行奇异值分解。

例：

> A=matrix(1:18,3,6)

> A

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

[1,]   1    4    7  10   13   16

[2,]   2    5    8  11   14   17

[3,]   3    6    9  12   15   18

> A.s=svd(A)

> A.s

$d

[1] 4.589453e+01 1.640705e+00 3.627301e-16

 

$u

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

[1,] -0.5290354  0.74394551 0.4082483

[2,] -0.5760715  0.03840487 -0.8164966

[3,] -0.6231077 -0.66713577  0.4082483

 

$v

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

[1,] -0.07736219 -0.71960032 -0.18918124

[2,] -0.19033085 -0.50893247  0.42405898

[3,] -0.30329950 -0.29826463 -0.45330031

[4,] -0.41626816 -0.08759679 -0.01637004

[5,] -0.52923682  0.12307105 0.64231130

[6,] -0.64220548  0.33373889 -0.40751869

> A.s$u%*%diag(A.s$d)%*%t(A.s$v)

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

[1,]   1    4    7  10   13   16

[2,]   2    5    8  11   14   17

[3,]   3    6    9  12   15   18

 

5.与矩阵（数组）运算有关的函数

（1）取矩阵的维数

函数dim()返回一个矩阵的维数，nrow()返回行数，ncol()返回列数。例：

> A=matrix(1:12,3,4)

> A

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

[1,]   1    4    7  10

[2,]   2    5    8  11

[3,]   3    6    9  12

> dim(A)

[1] 3 4

> nrow(A)

[1] 3

> ncol(A)

[1] 4

（2）矩阵的拉直

as.vector()函数可将矩阵转化为向量，如：

> A<-matrix(1:6,nrow=2);A

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

[1,]   1    3    5

[2,]   2    4    6

> as.vector(A)

[1] 1 2 3 4 5 6

（3）数组的维名字

数组可以有一个属性dimnames保存各维的各个下标的名字，缺省为NULL。如：

>x<-matrix(1:6,ncol=2,dimnames=list(c("A","B","C"),c("a","b")),byrow=T);x

  a b

A 1 2

B 3 4

C 5 6

也可以先定义矩阵x，然后再为dimnames(x)赋值。如：

> x<-matrix(1:6,ncol=2,byrow=T)

>dimnames(x)<-list(c("A","B","C"),c("a","b"))

对于矩阵，还可以使用属性rownames和colnames来访问行名和列名，例如：

> x<-matrix(1:6,ncol=2,byrow=T)

>colnames(x)<-c("a","b")

>rownames(x)<-c("A","B","C")

（4）apply()函数

对于向量，可以用sum、mean等函数对其进行计算，对于数组（矩阵），如果想对其一维（或若干维）进行某种计算，可用apply()函数，其一般形式为

apply(A,margin,fun)

其中A是一个数组（矩阵），margin是固定哪些维不变，fun是用来计算的函数，如：

> A<-matrix(1:6,nrow=2);A

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

[1,]   1    3    5

[2,]   2    4    6

> apply(A,1,sum)

[1] 9 12

> apply(A,2,mean)

[1] 1.5 3.5 5.5

参考：《统计建模与R软件》 薛毅 陈立萍 编著 清华大学出版社