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软件》 薛毅 陈立萍 编著 清华大学出版社