R语言编程艺术# 矩阵(matrix)和数组(array)

时间:2023-03-09 13:03:31
R语言编程艺术# 矩阵(matrix)和数组(array)

矩阵(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时,需要的矩阵如下:  

R语言编程艺术# 矩阵(matrix)和数组(array)

> 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 >