矩阵(matrix)是一种特殊的向量,包含两个附加的属性:行数和列数。所以矩阵也是和向量一样,有模式(数据类型)的概念。(但反过来,向量却不能看作是只有一列或一行的矩阵。
数组(array)是R里更一般的对象,矩阵是数组的一个特殊情形。数组可以是多维的。例如:一个三维数组可以包含行、列和层(layer),而一个矩阵只有行和列两个维度
1、创建矩阵
矩阵的行和列的下标都是从1开始,如:矩阵a左上角的元素记作a[1,1]。矩阵在R中是按列存储的,也就是说先存储第一列,再存储第二列,以此类推。
> y <- matrix(c(1,2,3,4),nrow=2,ncol=2) > y [,1] [,2] [1,] 1 3 [2,] 2 4 > y <- matrix(c(1,2,3,4),nrow=2) > y [,1] [,2] [1,] 1 3 [2,] 2 4 > #按列输出 > y[,2] #输出第二列 [1] 3 4 >
为矩阵中的元素赋值
> y <- matrix(nrow = 2,ncol = 2) > y [,1] [,2] [1,] NA NA [2,] NA NA > y[1,1] <- 1 > y[2,1] <- 2 > y[1,2] <- 3 > y[2,2] <-4 > y [,1] [,2] [1,] 1 3 [2,] 2 4 >
>#与上面的代码效果相同 > y <- matrix(c(1,2,3,4),nrow = 2) > y [,1] [,2] [1,] 1 3 [2,] 2 4 >
默认在R中矩阵是以列进行存储的,但通过byrow = T,参数可以将矩阵进行按行存储
> y <- matrix(c(1,2,3,4),nrow = 2, byrow = T) > y [,1] [,2] [1,] 1 2 [2,] 3 4 >
2、一般矩阵运算
常用的矩阵运算:线性代数运算、矩阵索引、矩阵元素筛选
#线性代数运算
线性代数运算包括:矩阵相乘、矩阵数量乘法、矩阵加法等
> y <- matrix(c(1,2,3,4),nrow = 2) > y [,1] [,2] [1,] 1 3 [2,] 2 4 > y %*% y #矩阵相乘 [,1] [,2] [1,] 7 15 [2,] 10 22 > 3*y #矩阵数量乘法 [,1] [,2] [1,] 3 9 [2,] 6 12 > y+y #矩阵加法 [,1] [,2] [1,] 2 6 [2,] 4 8 >
#矩阵索引
> z <- matrix(c(1,2,3,4,1,1,0,0,1,0,1,0),nrow = 4) > z [,1] [,2] [,3] [1,] 1 1 1 [2,] 2 1 0 [3,] 3 0 1 [4,] 4 0 0 > z[,2:3] #提取z中第2、3更 [,1] [,2] [1,] 1 1 [2,] 1 0 [3,] 0 1 [4,] 0 0 >
给矩阵赋值
> z [,1] [,2] [,3] [1,] 1 1 1 [2,] 2 1 0 [3,] 3 0 1 [4,] 4 0 0 > z[c(1,3),] <-matrix(c(1,1,8,12,16,20),nrow = 2) #给z1,3行进行赋新值 > z [,1] [,2] [,3] [1,] 1 8 16 [2,] 2 1 0 [3,] 1 12 20 [4,] 4 0 0 >
利用行号负值,移除行或列
> y <- matrix(c(1,2,3,4,5,6),nrow = 3) > y [,1] [,2] [1,] 1 4 [2,] 2 5 [3,] 3 6 > y[-2,] #移除第2行 [,1] [,2] [1,] 1 4 [2,] 3 6 > y[,-2] #移除第2列 [1] 1 2 3 >
#矩阵元素筛选
矩阵跟向量样也可以进行筛选,只是语法上不同而已
> x <-matrix(c(1,2,3,2,3,4),nrow = 3) > x [,1] [,2] [1,] 1 2 [2,] 2 3 [3,] 3 4 > x[x[,2]>=3,] #x中第2列所有大于等于3的行 [,1] [,2] [1,] 2 3 [2,] 3 4 >
矩阵筛选规则可以基于除被筛选变量这外的变量
> x [,1] [,2] [1,] 1 2 [2,] 2 3 [3,] 3 4 > z <- c(5,12,13) > x[z %% 2 == 1,] [,1] [,2] [1,] 1 2 [2,] 3 4 >
运算符:& and && 前者是向量的逻辑“与”运算,后者是用于if语句的标量逻辑“与”运算
> m <- matrix(c(1,2,3,4,5,6),nrow = 3) > m [,1] [,2] [1,] 1 4 [2,] 2 5 [3,] 3 6 > m[m[,1]>1 & m[,2]>5] #m中第1列中大于1,第2列中大于5的行 [1] 3 6 >
#扩展案例:生成协方差矩阵
n元正态分布,协方差矩阵有n行n列,要求n个随机变量方差都为1,每两个变量间的相关性都是rho,如:当n=3,rho=0.2时,需要的矩阵如下:
> makecov <- function(rho,n){ m<-matrix(nrow = n,ncol = n) m<-ifelse(row(m)==col(m),1,rho) return(m) } > makecov(0.2,3) [,1] [,2] [,3] [1,] 1.0 0.2 0.2 [2,] 0.2 1.0 0.2 [3,] 0.2 0.2 1.0 >
3、对矩阵的行和列调用函数
apply()函数,是R中最常用的函数,其中包括apply()、tapply()、lapply(),apply()函数允许用户在矩阵和各行或各列上调用指的函数。
apply()函数一般形式:apply(m,dimcode,f,fargs)
参数解释:
m:是一个矩阵
dimcode:是维度编号,若取值为1代表对一行应用函数,若取值为2代表对每一列应用函数
f:是应用在行或列上的函数
fargs:是f的可选参数
>######对z变量列进行mean()函数操作,做平均数计算 > z <- matrix(c(1,2,3,4,5,6),ncol = 2) > z [,1] [,2] [1,] 1 4 [2,] 2 5 [3,] 3 6 > apply(z,2,mean) [1] 2 5 >
>######当然上面的代码也可以有更简便的代码 > colMeans(z) [1] 2 5 > 其它语法参考如下: colSums (x, na.rm = FALSE, dims = 1) rowSums (x, na.rm = FALSE, dims = 1) colMeans(x, na.rm = FALSE, dims = 1) rowMeans(x, na.rm = FALSE, dims = 1) .colSums(x, m, n, na.rm = FALSE) .rowSums(x, m, n, na.rm = FALSE) .colMeans(x, m, n, na.rm = FALSE) .rowMeans(x, m, n, na.rm = FALSE)
#当然在R中apply()函数还可以使用自定义函数
> z [,1] [,2] [1,] 1 4 [2,] 2 5 [3,] 3 6 > f <- function(x) x/c(2,8) > y <- apply(z,1,f) #对z变量的行进行f函数操作 > y [,1] [,2] [,3] [1,] 0.5 1.000 1.50 [2,] 0.5 0.625 0.75 >
上面的代码输出的结果有两个重要的知识点:
1、如果向量x的长度大于2,那么(2,8)就会循环补齐,apply()对z的每行分别调用f(),形参x对应用的实参是(1,4)。
2、y输出的结果是一个2x3的矩阵而不是z一样的3x2的矩阵,因为R中的矩阵默认是以列进行存储的,所以当第一行输出的结果自然也是按列进行存储,如果调用f()返回有k个元素向量,那么apply()的结果就有k行。但是可以通过t()函数进行行列转置。
> t(apply(z,1,f)) [,1] [,2] [1,] 0.5 0.500 [2,] 1.0 0.625 [3,] 1.5 0.750 >
#所调用的函数只返回一个标量(即单个元素向量),那么apply()的结果就是一个向量,而非矩阵,在使用apply()函数时调用的函数至少需一个参数,在上例中的形参对应的实参就是z矩阵中的一行(或一列),有时待调用的函数需要多个参数,在调用这类函数时,调用的函数的参数写在函数名称的后面用逗号隔开。
> copymaj <- function(rw,d) { + maj <- sum(rw[1:d]) / d + return(ifelse(maj > 0.5,1,0)) + } > x <- matrix(c(1,1,1,0,0,1,0,1,1,1,0,1,1,1,1,1,0,0,1,0),nrow = 4) > x [,1] [,2] [,3] [,4] [,5] [1,] 1 0 1 1 0 [2,] 1 1 1 1 0 [3,] 1 0 0 1 1 [4,] 0 1 1 1 0 > apply(x,1,copymaj,3) [1] 1 1 0 1 > apply(x,1,copymaj,2) [1] 0 1 0 0 >
在R中使用apply()函数不能使程序运行速度加快,其优点是使代码更紧凑,便于阅读和修改,避免产生使用循环语句时可能带来的bug。此外并行运算是R目前发展的方向之一,apply()这类函数会变得越来越重要。如:在sonw包中的clusterApplay()函数能够把子矩阵的数据分配到多个网络节点上,在每个网络节点上对子矩阵调用给定的函数,达到并行计算的目的。
#扩展案例:寻找异常值
在统计学中,“异常值”(outlier)指的是哪些和大多数观测值离得很远的少数点。所以异常值要么是有问题(例如数字写错了),要么是不具有代表性(例如比尔盖茨的收入和华盛顿居民的收入相比),通常用到median()函数(中位数函数)
中位数(又称中值,英语:Median),统计学中的专有名词,代表一个样本、种群或概率分布中的一个数值,其可将数值集合划分为相等的上下两部分。对于有限的数集,可以通过把所有观察值高低排序后找出正中间的一个作为中位数。如果观察值有偶数个,通常取最中间的两个数值的平均数作为中位数。
> findols function(x){ findol <- function(xrow){ mdn <- median(xrow) devs <- abs(xrow - mdn) return(which.max(devs)) } return(apply(x,1,findol)) } > x [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 1 8 9 20 20 24 25 [2,] 2 13 13 17 18 19 26 [3,] 6 10 11 12 19 20 31 [4,] 4 5 12 13 24 24 28 [5,] 5 6 17 17 21 22 23 [6,] 10 10 14 15 16 23 24 [7,] 7 8 9 16 17 28 28 [8,] 2 9 10 21 21 25 26 [9,] 3 14 14 18 19 20 27 [10,] 7 11 12 13 20 21 32 [11,] 5 6 13 14 25 25 29 [12,] 6 7 18 18 22 23 24 [13,] 11 11 15 16 17 24 25 [14,] 8 9 10 17 18 29 29 [15,] 3 10 11 22 22 26 27 [16,] 4 15 15 19 20 21 28 [17,] 8 12 13 14 21 22 33 [18,] 6 7 14 15 26 26 30 [19,] 7 8 19 19 23 24 25 [20,] 12 12 16 17 18 25 26 [21,] 9 10 11 18 19 30 30 [22,] 4 11 12 23 23 27 28 [23,] 5 16 16 20 21 22 29 [24,] 9 13 14 15 22 23 34 [25,] 7 8 15 16 27 27 31 > findols(x) [1] 1 1 7 7 1 7 6 1 1 7 7 1 7 6 1 1 7 7 1 7 6 1 1 7 7 #输出为异常数的位置 >
4、增加或删除矩阵的行或列
严格来说,矩阵的长度和维度是固定的,因此不能增加或删除行或列,但是可以给矩阵重新赋值,这样可以得到和增加删除一样的效果
#改变矩阵的大小
>#####向量的增、插、删 > x <- c(1,2,3,4) > x [1] 1 2 3 4 > x <- c(x,99) #增加 > x [1] 1 2 3 4 99 > x <- c(x[1:4],88,x[5]) #插入 > x [1] 1 2 3 4 88 99 > x <- x[-4:-5] #删除第4:5个元素 > x [1] 1 2 3 99 >
改变矩阵常用到的函数rbind()、cbind(),可以给矩阵增加行或列
语法:
cbind(..., deparse.level = 1)
rbind(..., deparse.level = 1)
## S3 method for class 'data.frame'
rbind(..., deparse.level = 1, make.row.names = TRUE, stringsAsFactors = default.stringsAsFactors())
> x <- c(1,1,1) > x [1] 1 1 1 > z <- matrix(c(1,2,3,4,5,6,7,8,9), nrow = 3) > z [,1] [,2] [,3] [1,] 1 4 7 [2,] 2 5 8 [3,] 3 6 9 > cbind(x,z) x [1,] 1 1 4 7 [2,] 1 2 5 8 [3,] 1 3 6 9 > cbind(z,x) x [1,] 1 4 7 1 [2,] 2 5 8 1 [3,] 3 6 9 1 > cbind(9,z) [,1] [,2] [,3] [,4] [1,] 9 1 4 7 [2,] 9 2 5 8 [3,] 9 3 6 9 >
函数cbind()、rbind()还可以用来快速生成一些小的矩阵
> q <- cbind(c(1,2),c(3,4)) > q [,1] [,2] [1,] 1 3 [2,] 2 4 >
不过!不要高兴太早了,以会有了cbind,rbind对矩阵增、删就方便了,但事实你将要付出更大的资源,和创建向量一样,创建一个新的矩阵是很耗时间的(毕竟矩阵也属于向量),假如要在矩阵中插入10w+条记录,相当于将矩阵进行了10w+的增、删。
不过不要悲观,我们可以预先创建一个足够大的矩阵(按需),最开始矩阵是空的(NA)然后在循环过程中逐行或逐列进行赋值,这样做法避免了循环过程中每次进行耗时的矩阵内存分配。
> m <- matrix(nrow = 3,ncol = 2) > m [,1] [,2] [1,] NA NA [2,] NA NA [3,] NA NA > m[,] <- c(c(1:3),c(4:6)) > m [,1] [,2] [1,] 1 4 [2,] 2 5 [3,] 3 6 > m <- m[c(1,3),] > m [,1] [,2] [1,] 1 4 [2,] 3 6 >
#扩展案例:找到图中距离最近的一对端点
计算图中多个端点之间距离是计算机或统计学中常见的例子,这类问题在聚类算法和基因问题中经常出现。
我们以计算城市之间的距离为例,这比计算DNA链间距离更直观。
假设有一个距离矩阵,其第i行第j列的元素代表城市i和城市j间的距离。我们需要写一个函数,输入城市距离矩阵,输出城市间最短的距离,以及对应的两个城市。
mind <- function(d){ n <- nrow(d) dd <- cbind(d,1:n) wmins <- apply(dd[-n,],1,imin) i <- which.min(wmins[2,]) j <- wmins[1,i] return(c(d[i,j],i,j)) } imin <-function(x) { lx <- length(x) i <- x[lx] j <- which.min(x[(i+1):(lx-1)]) k <- i+j return(c(k,x[k])) } q <- matrix(c(0,12,13,8,20,12,0,15,28,88,13,15,0,6,9,8,28,6,0,33,20,88,9,33,0),nrow = 5) > q [,1] [,2] [,3] [,4] [,5] [1,] 0 12 13 8 20 [2,] 12 0 15 28 88 [3,] 13 15 0 6 9 [4,] 8 28 6 0 33 [5,] 20 88 9 33 0 > mind(q) [1] 6 3 4 #最小值是6,位于在第3行第4列 >
5、向量与矩阵的差异
矩阵就是一个向量,只是多了两个属性:行娄和列数
从面向对象编程的角度来说,矩阵类(matrix class)是实际存在的,R中的大部分类都是S3类,用$符号就可以访问其各组件。矩阵类有一个dim属性,是一个由矩阵的行数和列数组成的向量,可以用dim()函数访问dim属性。
> z <-matrix(1:8,nrow = 4) > z [,1] [,2] [1,] 1 5 [2,] 2 6 [3,] 3 7 [4,] 4 8 > length(z) [1] 8 > class(z) [1] "matrix" > attributes(z) $dim [1] 4 2 > y <-c(1:8) > y [1] 1 2 3 4 5 6 7 8 > length(y) [1] 8 > class(y) [1] "integer" > attributes(y) NULL >
> z [,1] [,2] [1,] 1 5 [2,] 2 6 [3,] 3 7 [4,] 4 8 > dim(z) [1] 4 2 > nrow(z) [1] 4 > ncol(z) [1] 2 > nrow function (x) dim(x)[1L] <bytecode: 0x07b80efc> <environment: namespace:base> > x <- c(1:12) ;dim(x)<-c(3,4) > x [,1] [,2] [,3] [,4] [1,] 1 4 7 10 [2,] 2 5 8 11 [3,] 3 6 9 12 >
6、避免意外降维
在统计学领域,“降维”(dimension reduction)是有益的,也存在很多降维的统计学方法。假设我们需要处理10个变量,如果能把变量个数降到3个,却还能保留数据的主要信息,何乐而不为呢?
在R中,降维指的是完全另外一件事情,而且通常要避免。
> z [,1] [,2] [1,] 1 5 [2,] 2 6 [3,] 3 7 [4,] 4 8 > r <- z[2,] > r [1] 2 6 > attributes(z) $dim [1] 4 2 > attributes(r) NULL > str(z) int [1:4, 1:2] 1 2 3 4 5 6 7 8 > str(r) int [1:2] 2 6 >
从上面的代码可以看出,r的结果显示的是向量格式,而非矩阵的格式,也就是说,r是一个长度为2的向量,而不是一个1*2的矩阵
在R中可以使用drop参数,禁止矩阵自动减少维度。
> r <- z[2,,drop = FALSE] > r [,1] [,2] [1,] 2 6 > dim(r) [1] 1 2 >
对原本就是向量的对象,可以使用as.matrix()函数将其转化成矩阵
> u <- c(1:12) > u [1] 1 2 3 4 5 6 7 8 9 10 11 12 > v <- as.matrix(u) > v [,1] [1,] 1 [2,] 2 [3,] 3 [4,] 4 [5,] 5 [6,] 6 [7,] 7 [8,] 8 [9,] 9 [10,] 10 [11,] 11 [12,] 12 > attributes(v) $dim [1] 12 1
7、矩阵的行和列的命名问题
访问矩阵元素最直接的方法就是通过行号和列号,但也可以使用行名与列名
> z [,1] [,2] [1,] 1 5 [2,] 2 6 [3,] 3 7 [4,] 4 8 > colnames(z) NULL > colnames(z) <- c("a","b") > z a b [1,] 1 5 [2,] 2 6 [3,] 3 7 [4,] 4 8 > colnames(z) [1] "a" "b" > z[,"a"] [1] 1 2 3 4 >
8、高维数组
在统计学领域,R语言中典型的矩阵用行表示不同的观测,比如不同的人,而用列表示不同变量,比如体重血压等。因此矩阵一般都是二维的数据结构。但是假如我们的数据采集自不同的时间,也就是每个人每个变量每个时刻记录一个数。时间就成为除了行和列之外的第三个维度,在R中,这样的数据称为数组(arrays)。
> firsttest <- matrix(c(46,21,50,30,25,50), nrow = 3) > firsttest [,1] [,2] [1,] 46 30 [2,] 21 25 [3,] 50 50 > secondtest <- matrix(c(46,41,50,43,35,50), nrow = 3) > secondtest [,1] [,2] [1,] 46 43 [2,] 41 35 [3,] 50 50 > tests <- array(data = c(firsttest,secondtest),dim = c(3,2,2)) > attributes(tests) $dim [1] 3 2 2 > tests[3,2,1] #第3行,第2列,第1个表 [1] 50 > tests[2,2,1] #第2行,第2列,第1个表 [1] 25 > tests[2,2,2] #第2行,第2列,第2个表 [1] 35 >
tests共分为两个数据层(layer),一层对应一次考试,每层都是3*2的矩阵
> tests , , 1 [,1] [,2] [1,] 46 30 [2,] 21 25 [3,] 50 50 , , 2 [,1] [,2] [1,] 46 43 [2,] 41 35 [3,] 50 50 >