2.5 多维数组和矩阵
2.5.1 生成数组或矩阵
数组有一个特征属性叫做维数向量(dim属性)，维数向量是一个元素取正整数的向量，其长度是数组的维数，比如维数向量有两个元素时数组为2维数组（矩阵）。维数向量的每一个元素指定了该下标的上界，下标的下界总为1
1.将向量定义成数组
向量只有定义了维数向量（dim属性）后才能被看作是数组
> 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()函数构造多维数组
R软件可以用array()函数直接构造数组，其构造形式为
array(data=NA,dim=length(data),dimnames=NULL)
其中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.用matrix()函数构造矩阵
函数matrix()是构造矩阵的函数，其构造形式为
matrix(data=NA,nrow=1,ncol=1,byrow=FALSE,dimnames=NULL)
其中data是一个向量数据，nrow是矩阵的行数，ncol是矩阵的列数。当byrow=TRUE时，生成矩阵的数据按行放置，缺省时相当于byrow=FALSE，数据按列放置。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
[,1] [,2] [,3] [,4] [,5]
[1,] 1 2 3 4 5
[2,] 6 7 8 9 10
[3,] 11 12 13 14 15
2.5.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
> a
, , 1
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
, , 2
[,1] [,2] [,3]
[1,] 7 9 11
[2,] 8 10 12
, , 3
[,1] [,2] [,3]
[1,] 13 15 17
[2,] 14 16 18
, , 4
[,1] [,2] [,3]
[1,] 19 21 23
[2,] 20 22 24
> a[1,,]
[,1] [,2] [,3] [,4]
[1,] 1 7 13 19
[2,] 3 9 15 21
[3,] 5 11 17 23
> a[,2,]
[,1] [,2] [,3] [,4]
[1,] 3 9 15 21
[2,] 4 10 16 22
> a[1,1,]
[1] 1 7 13 19
> a[3:10]
[1] 3 4 5 6 7 8 9 10
2.不规则的数组下标
在R语言中，甚至可以把数组中的任意位置的元素作为数组访问，其方法是用一个二维数组作为数组的下标，二维数组的每一个行是一个元素的下标，列数为数组的维数。
例如，要把上面的形状为234的数组a的第[1,1,1],[2,2,3],[1,3,4],[2,1,4]号共四个元素作为一个整体访问，先定义一个包含这些下标作为行的二维数组：
> b<- matrix(c(1,1,1,2,2,3,1,3,4,2,1,4),ncol=3,byrow = TRUE);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[b]
[1] 101 102 103 104
> a
, , 1
[,1] [,2] [,3]
[1,] 101 3 5
[2,] 2 4 6
, , 2
[,1] [,2] [,3]
[1,] 7 9 11
[2,] 8 10 12
, , 3
[,1] [,2] [,3]
[1,] 13 15 17
[2,] 14 102 18
, , 4
[,1] [,2] [,3]
[1,] 19 21 103
[2,] 104 22 24
2.5.3 数组的四则运算
> A<- matrix(1:6,nrow = 2,byrow = TRUE);A
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
> B<- matrix(1:6,nrow = 2);B
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
> C<- matrix(c(1,2,2,3,3,4),nrow = 2);C
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 2 3 4
> D<- 2*C+A/B;D
[,1] [,2] [,3]
[1,] 3 4.666667 6.6
[2,] 6 7.250000 9.0
形状不一致的向量（或数组）也可以进行四则运算，一般的规则是将向量（或数组）中的数据与对应向量（或数组）中的数据进行运算，把短向量（或数组）的数据循环使用，从而可以与长向量（或数组）数据进行匹配，并尽可能保留共同的数组属性
> x1<-c(100,200)
> x2<-1:6
> x1+x2
[1] 101 202 103 204 105 206
> x3<- matrix(1:6,nrow = 3)
> x1+x3
[,1] [,2]
[1,] 101 204
[2,] 202 105
[3,] 103 206
> x2<- 1:5
> x1+x2
[1] 101 202 103 204 105
Warning message:
In x1 + x2 : 长的对象长度不是短的对象长度的整倍数
2.5.4 矩阵的运算
1.转置运算
> A<-matrix(1:6,nrow=2);A
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
> t(A)
[,1] [,2]
[1,] 1 2
[2,] 3 4
[3,] 5 6
> ### 2.求方阵的行列式
> det(matrix(1:4,ncol=2))
[1] -2
> ### 3.向量的内积
> x<-1:5;y<-2*1:5
> x%*%y
[,1]
[1,] 110
函数crossprod()是内积运算函数（表示交叉乘积），crossprod(x,y)计算向量x与y的内积，即't(x) %*% y'
> crossprod(x,y)
[,1]
[1,] 110
> x;t(x)
[1] 1 2 3 4 5
[,1] [,2] [,3] [,4] [,5]
[1,] 1 2 3 4 5
crossprod(x)表示x与x的内积，即(||x||^2_2)
类似地，tcrossprod(x,y)表示'x%*% t(y)',即x与y的外积，也称为叉积。tcrossprod(x)表示x与x作外积
> tcrossprod(x,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
> tcrossprod(x)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 2 3 4 5
[2,] 2 4 6 8 10
[3,] 3 6 9 12 15
[4,] 4 8 12 16 20
[5,] 5 10 15 20 25
4.向量的外积（叉积）
设x,y是n维向量，则x %o% y表示x与y作外积
> x<-1:5;y<-2*1:5
> 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()在绘制三维曲面时非常有用，它可以生成一个X和Y的网格，
> outer(x,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
5.矩阵的乘法
如果矩阵A和B具有相同的维数，则AB表示矩阵中对应的元素的乘积，A %% B表示通常意义下的两个矩阵的乘积（要求矩阵A的列数等于矩阵B的行数）
> A<-array(1:9,dim=(c(3,3)))
> B<-array(9:1,dim=(c(3,3)))
> C<-A*B;C
[,1] [,2] [,3]
[1,] 9 24 21
[2,] 16 25 16
[3,] 21 24 9
> D<- A%*%B;D
[,1] [,2] [,3]
[1,] 90 54 18
[2,] 114 69 24
[3,] 138 84 30
crossprod(A,B)表示的是t(A)%%B,函数tcrossprod(A,B)表示的是A%%t(B)
> A
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> B<-c(1:3);B
[1] 1 2 3
> crossprod(B,A)
[,1] [,2] [,3]
[1,] 14 32 50
> t(B);crossprod(t(B),A)
[,1] [,2] [,3]
[1,] 1 2 3
Error in crossprod(t(B), A) : 非整合参数
> tcrossprod(A,B)
Error in tcrossprod(A, B) : 非整合参数
> tcrossprod(A,t(B))
[,1]
[1,] 30
[2,] 36
[3,] 42
6.生成对角阵和矩阵取对角运算
函数diag()依赖于它的变量，当v是一个向量时，diag(v)表示以v的元素为对角线元素的对角阵。当M是一个矩阵时，则diag(M)表示的是取M对角线上的元素的向量
> v<-c(1,4,5)
> diag(v)
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 4 0
[3,] 0 0 5
> M<-array(1:9,dim=c(3,3));M
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> diag(M)
[1] 1 5 9
7.解线性方程组和求矩阵的逆矩阵
若求解线性方程组Ax=b，其命令形式为solve(A,b),求矩阵A的逆，其命令形式为solve(A)
> A<-t(array(c(1:8,10),dim=c(3,3)))
> b<-c(1,1,1)
> x<-solve(A,b);x
[1] -1.000000e+00 1.000000e+00 3.330669e-16
> B<-solve(A);B
[,1] [,2] [,3]
[1,] -0.6666667 -1.333333 1
[2,] -0.6666667 3.666667 -2
[3,] 1.0000000 -2.000000 1
> A %*% B
[,1] [,2] [,3]
[1,] 1 8.881784e-16 -4.440892e-16
[2,] 0 1.000000e+00 -1.776357e-15
[3,] 0 0.000000e+00 1.000000e+00
8.求矩阵的特征值与特征值向量
eigen(Sm)是求对称矩阵Sm的特征值与特征向量，其命令形式为
ev<-eigen(Sm)
ev存放着对称矩阵Sm特征值和特征向量，是由列表形式给出，其中ev(values是Sm的特征值构成的向量，ev)vectors是Sm的特征向量构成的矩阵
> A
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
[3,] 7 8 10
> Sm<-crossprod(A,A);Sm
[,1] [,2] [,3]
[1,] 66 78 97
[2,] 78 93 116
[3,] 97 116 145
> ev<-eigen(Sm);ev
eigen() decomposition
$values
[1] 303.19533618 0.76590739 0.03875643
$vectors
[,1] [,2] [,3]
[1,] -0.4646675 0.833286355 0.2995295
[2,] -0.5537546 -0.009499485 -0.8326258
[3,] -0.6909703 -0.552759994 0.4658502
9.矩阵的奇异值分解
函数svd(A)是对称A作奇异值分解，即A=UDV',其中U,V是正交阵，D为对角阵，也就是矩阵A的奇异值。svd(A)的返回值也是列表，svd(A)(d表示矩阵A的奇异值，即矩阵D的对角线上的元素，svd(A))u对应的是正交阵U，svd(A)$v对应的是正交阵V
> svdA<-svd(A);svdA
$d
[1] 17.4125052 0.8751614 0.1968665
$u
[,1] [,2] [,3]
[1,] -0.2093373 0.96438514 0.1616762
[2,] -0.5038485 0.03532145 -0.8630696
[3,] -0.8380421 -0.26213299 0.4785099
$v
[,1] [,2] [,3]
[1,] -0.4646675 -0.833286355 0.2995295
[2,] -0.5537546 0.009499485 -0.8326258
[3,] -0.6909703 0.552759994 0.4658502
> attach(svdA)
The following object is masked _by_ .GlobalEnv:
v
> u %*% diag(d) %*% t(v)
Error in u %*% diag(d) %*% t(v) : 非整合参数
10.求矩阵行列式的值
> det(A)
[1] -3
11.最小拟合与QR分解
> x<-c(0.0,0.2,0.4,0.6,0.8)
> y<-c(0.9,1.9,2.8,3.3,4.2)
> lsfit.sol<-lsfit(x,y);lsfit.sol
$coefficients
Intercept X
1.02 4.00
$residuals
[1] -0.12 0.08 0.18 -0.12 -0.02
$intercept
[1] TRUE
$qr
$qt
[1] -5.85849810 2.52982213 0.23749843 -0.02946714 0.10356728
$qr
Intercept X
[1,] -2.2360680 -0.8944272
[2,] 0.4472136 0.6324555
[3,] 0.4472136 -0.1954395
[4,] 0.4472136 -0.5116673
[5,] 0.4472136 -0.8278950
$qraux
[1] 1.447214 1.120788
$rank
[1] 2
$pivot
[1] 1 2
$tol
[1] 1e-07
attr(,"class")
[1] "qr"
(coefficients是拟合系数，)residuals是拟合残差
与lsfit()函数有密切关系的函数是ls.diag(),它给出拟合的进一步统计信息
另一个最小二乘拟合有密切关系的函数是QR分解函数qr(),qe.coef(),qr.fitted()和qr.resid()
> X<-matrix(c(rep(1,5),x),ncol = 2);X
[,1] [,2]
[1,] 1 0.0
[2,] 1 0.2
[3,] 1 0.4
[4,] 1 0.6
[5,] 1 0.8
> Xplus<-qr(X);Xplus
$qr
[,1] [,2]
[1,] -2.2360680 -0.8944272
[2,] 0.4472136 0.6324555
[3,] 0.4472136 -0.1954395
[4,] 0.4472136 -0.5116673
[5,] 0.4472136 -0.8278950
$rank
[1] 2
$qraux
[1] 1.447214 1.120788
$pivot
[1] 1 2
attr(,"class")
[1] "qr"
QR分解函数qr()输入的设计矩阵需要加以1为元素的列，其返回值是列表，其中(qr矩阵的上三角阵是QR分解中得到的R矩阵，下三角阵是QR分解得到的正交阵Q的部分信息，)qraux是Q的附加信息
可用QR分解得到的结果计算最小二乘的系数
> b<- qr.coef(Xplus,y);b
[1] 1.02 4.00
得到的系数与函数lsfit()得到的结果相同，但是为什么用这种方法计算呢？这是因为用QR分解在计算最小二乘拟合时，其计算误差比一般方法要小
类似的，可以用QR分解得到最小二乘的拟合值和残差值
> fit<-qr.fitted(Xplus,y);fit
[1] 1.02 1.82 2.62 3.42 4.22
> res<-qr.resid(Xplus,y);res
[1] -0.12 0.08 0.18 -0.12 -0.02
2.5.5 与矩阵（数组）运算相关的函数
1.取矩阵的维数
> A<-matrix(1:6,nrow=2);A
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
> dim(A)
[1] 2 3
> nrow(A)
[1] 2
> ncol(A)
[1] 3
2.矩阵的合并
函数cbind()把其自变量横向拼成一个大矩阵，rbind()把其自变量纵向拼成一个大矩阵
> 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
3.矩阵的拉直
> 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
4.数组的维名字
> 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<-matrix(1:6,ncol=2,byrow=T)
> dimnames(X)<-list(c("one","two","three"),c("First","Second"))
> colnames(X)
[1] "First" "Second"
> rownames(X)
[1] "one" "two" "three"
5.矩阵的广义转置
可以用aperm(A,perm)函数把数组A的各维按perm中指定的新次序重新排列
> A<-array(1:24,dim=c(2,3,4));A
, , 1
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
, , 2
[,1] [,2] [,3]
[1,] 7 9 11
[2,] 8 10 12
, , 3
[,1] [,2] [,3]
[1,] 13 15 17
[2,] 14 16 18
, , 4
[,1] [,2] [,3]
[1,] 19 21 23
[2,] 20 22 24
> B<-aperm(A,c(2,3,1));B
, , 1
[,1] [,2] [,3] [,4]
[1,] 1 7 13 19
[2,] 3 9 15 21
[3,] 5 11 17 23
, , 2
[,1] [,2] [,3] [,4]
[1,] 2 8 14 20
[2,] 4 10 16 22
[3,] 6 12 18 24
结果是B把A的第2维移到了第1维，A的第3维移到了第2维，A的第1维移到了第3维，这时有B[i,j,k]=A[i,k,i]
> B[1,3,2];A[3,2,1]
[1] 14
Error in A[3, 2, 1] : 下标出界
6.apply函数
对于向量，可以用sum,mean等函数对其进行计算，对于数组（矩阵），如果想对其1维（或若干维）进行某种计算，可用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