本文内容来自《R 语言编程艺术》(The Art of R Programming),有部分修改
矩阵 (Matrix) 是一种特殊的向量,包含两个附加属性:行数和列数。
数组 (Array) 是更一般的对象,可以有多个维度。矩阵是二维数组。
创建矩阵
R 中下标从 1 开始,矩阵按列存储。
使用 matrix()
函数创建矩阵
matrix(c(1, 2, 3, 4),
nrow=2,
ncol=2,
)print(y)
[,1] [,2]
[1,] 1 3
[2,] 2 4
在给定所有元素的情况下,nrow
和 ncol
两个参数只需要给出一个
matrix(1:4,
nrow=2
)print(y)
[,1] [,2]
[1,] 1 3
[2,] 2 4
另一种方法是创建空的矩阵,再对每一位赋值
matrix(nrow=2, ncol=2)print(y)
[,1] [,2]
[1,] NA NA
[2,] NA NA
1]
[,1] [,2]
[1,] 1 3
[2,] 2 4
设置参数 byrow=TRUE
,可以按行提供的矩阵数据。但矩阵实际还是按列存储
matrix(1:6,
nrow=2,
byrow=TRUE,
)print(m)
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
一般矩阵运算
线性代数运算
矩阵乘法
print(y
[,1] [,2]
[1,] 7 15
[2,] 10 22
矩阵元素乘法
print(
[,1] [,2]
[1,] 3 9
[2,] 6 12
矩阵加法
print(y
[,1] [,2]
[1,] 2 6
[2,] 4 8
矩阵索引
向量索引同样适用于矩阵索引
提取矩阵列
matrix(c(1, 2, 3, 4,1, 1, 0, 0,1, 0, 1, 0
),
nrow=4,
)print(z)
[,1] [,2] [,3]
[1,] 1 1 1
[2,] 2 1 0
[3,] 3 0 1
[4,] 4 0 0
print(z[,
[,1] [,2]
[1,] 1 1
[2,] 1 0
[3,] 0 1
[4,] 0 0
提取矩阵行
matrix(c(11, 21, 31,12, 22, 32
),
ncol=2,
)print(y)
[,1] [,2]
[1,] 11 12
[2,] 21 22
[3,] 31 32
print(y[2
[,1] [,2]
[1,] 21 22
[2,] 31 32
print(y[2
[1] 22 32
对子矩阵赋值
matrix(1:6,
nrow=3
)print(y)
[,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6
[c(
[,1] [,2]
[1,] 1 8
[2,] 2 5
[3,] 1 12
另一个示例
matrix(nrow=3, ncol=3)
y matrix(c(4, 5, 2, 3),
nrow=2
)print(y)
[,1] [,2]
[1,] 4 2
[2,] 5 3
matrix(1:6,
nrow=3
)print(y)
[,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6
print(y[
[,1] [,2]
[1,] 1 4
[2,] 3 6
扩展案例:图像操作
library(pixmap)
Pixmap image
Type : pixmapGrey
Size : 194x259
Resolution : 1x1
Bounding box : 0 0 259 194
plot(mtrush1)
str(mtrush1)
Formal class 'pixmapGrey' [package "pixmap"] with 6 slots
..@ grey : num [1:194, 1:259] 0.278 0.263 0.239 0.212 0.192 ...
..@ channels: chr "grey"
..@ size : int [1:2] 194 259
..@ cellres : num [1:2] 1 1
..@ bbox : num [1:4] 0 0 259 194
..@ bbcent : logi FALSE
@grey[28,
0.796078431372549
mtrush1
mtrush2@grey[84:163, 135:177] 1plot(mtrush2)
function(img, rows, cols, q) {
lrows length(rows)
lcols length(cols)
newimg img
randomnoise matrix(nrow=lrows, ncol=lcols, runif(lrows*lcols))
newimg@grey[rows, cols] (1-q) * img@grey[rows, cols] + q * randomnoisereturn(newimg)
}
blurpart(
mtrush1,84:163,135:177,0.2
)plot(mtrush3)
矩阵元素筛选
matrix(c(1, 2, 3, 2, 3, 4),
nrow=3
)print(x)
[,1] [,2]
[1,] 1 2
[2,] 2 3
[3,] 3 4
print(x[x[,
[,1] [,2]
[1,] 2 3
[2,] 3 4
分析上面的命令
x[,2] >= 3print(j)
[1] FALSE TRUE TRUE
print(x[j,])
[,1] [,2]
[1,] 2 3
[2,] 3 4
注意:求 j 是向量化计算
另一个示例
matrix(1:6,
nrow=3
)print(m)
[,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6
c(5, 12, 13)print(m[z %% 2 == 1,])
[,1] [,2]
[1,] 1 4
[2,] 3 6
print(m[m[,
[1] 3 6
分析
m[,1] > 1print(a)
[1] FALSE TRUE TRUE
m[,2] > 5print(b)
[1] FALSE FALSE TRUE
a & bprint(c)
[1] FALSE FALSE TRUE
print(m[c, ])
[1] 3 6
注意:这里使用 &
表示向量的逻辑与运算。而 &&
是用于 if 语句的标量逻辑与运算
向量运算也适用于矩阵
matrix(c(5, 2, 9, -1, 10, 11),
nrow=3
)print(m)
[,1] [,2]
[1,] 5 -1
[2,] 2 10
[3,] 9 11
print(
[1] 1 3 5 6
扩展案例:生成协方差矩阵
row()
函数和 col()
函数
协方差矩阵是对称的,对角线元素均为 1
function(rho, n) {
m matrix(nrow=n, ncol=n)
m ifelse(row(m) == col(m), 1, rho)return(m)
}
示例
print(
[,1] [,2] [,3]
[1,] 1.0 0.2 0.2
[2,] 0.2 1.0 0.2
[3,] 0.2 0.2 1.0
解释
matrix(3:8,
nrow=3
)print(z)
[,1] [,2]
[1,] 3 6
[2,] 4 7
[3,] 5 8
print(
[,1] [,2]
[1,] 1 1
[2,] 2 2
[3,] 3 3
print(
[,1] [,2]
[1,] 1 2
[2,] 1 2
[3,] 1 2
print(
[,1] [,2]
[1,] TRUE FALSE
[2,] FALSE TRUE
[3,] FALSE FALSE
对矩阵的行和列调用函数
*apply 系列函数:apply, tapply, lapply
apply()
函数在各行或各列上调用指定函数
使用 apply()
函数
一般形式
apply(m, dimcode, f, fargs)
其中
m
:矩阵dimcode
:维度编号,1 表示对每行应用函数,2 表示对每列应用函数f
:应用在行或列上的函数,第一个参数为一行或一列fargs
:可选参数,用于传递给f
函数
matrix(1:6,
nrow=3
)print(z)
[,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6
对每一列应用函数 mean
注:可以直接使用 colmean()
函数
print(
[1] 2 5
对每一行应用自定义函数
function(x) x/c(2, 8)
y apply(z, 1, f)print(y)
[,1] [,2] [,3]
[1,] 0.5 1.000 1.50
[2,] 0.5 0.625 0.75
调用的函数返回 k 个元素的向量,则 apply
函数返回的结果就有 k 行。执行几次函数,则返回结果就有多少列。
如果函数返回一个标量,则 apply
函数返回向量。
t()
函数可以实现转置
print(
[,1] [,2]
[1,] 0.5 0.500
[2,] 1.0 0.625
[3,] 1.5 0.750
附加参数
0 和 1 构成的矩阵,对于每行,前 d 个元素中 1 较多则取 1,否则取 0
function(rw, d) {
maj sum(rw[1:d]) / dreturn(if(maj > 0.5) 1 else 0)
}
matrix(c(1, 1, 1, 0,0, 1, 0, 1,1, 1, 0, 1,1, 1, 1, 1,0, 0, 1, 0
),
nrow=4
)print(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
print(
[1] 1 1 0 1
print(
[1] 0 1 0 0
扩展案例:寻找异常值
识别偏离最远的观测值
function(x) {
findol function(xrow) {
mdn median(xrow)
devs abs(xrow-mdn)return(which.max(devs))
}return(apply(x, 1, findol))
}
增加或删除矩阵的行或列
矩阵无法改变大小,增加或删除实际上是创建新的对象
改变矩阵大小
重新赋值可以改变向量的大小
c(12, 5, 13, 16, 8)print(x)
[1] 12 5 13 16 8
c(x, 20)print(x)
[1] 12 5 13 16 8 20
c(x[1:3], 20, x[4:6])print(x)
[1] 12 5 13 20 16 8 20
x[-2:-4]print(x)
[1] 12 16 8 20
使用 rbind()
和 cbind()
为矩阵增加行或列
c(1, 1, 1, 1)
z matrix(c(1, 2, 3 ,4,1, 1, 0, 0,1, 0, 1, 0
),
nrow=4
)print(z)
[,1] [,2] [,3]
[1,] 1 1 1
[2,] 2 1 0
[3,] 3 0 1
[4,] 4 0 0
print(
one
[1,] 1 1 1 1
[2,] 1 2 1 0
[3,] 1 3 0 1
[4,] 1 4 0 0
可能会用到循环补齐 (recycling)
print(
[,1] [,2] [,3] [,4]
[1,] 1 1 1 1
[2,] 1 2 1 0
[3,] 1 3 0 1
[4,] 1 4 0 0
cbind()
和 rbind()
可以用来快速生成小矩阵
cbind(c(1, 2),c(3, 4)
)print(q)
[,1] [,2]
[1,] 1 3
[2,] 2 4
注意:创建向量和矩阵均耗费时间,请谨慎使用 cbind
和 rbind
循环添加列或行的方式不可取,应该首先生成一个大矩阵,然后再往里填值。
可以通过重新赋值删除矩阵的行或列
matrix(1:6, nrow=3)print(m)
[,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6
m[c(1,3),]print(m)
[,1] [,2]
[1,] 1 4
[2,] 3 6
扩展案例:找到图中距离最近的一对端点
矩阵元素 (i, j) 表示城市 i 和 城市 j 间的距离
imin
函数寻找当前行的最小值,因为矩阵是对称的,所以仅搜索从行号 i + 1 开始到结尾。
返回最小值的索引(即矩阵中的列数)和最小值
function(x) {
lx length(x)
i x[lx]
j which.min(x[(i+1):(lx-1)])
k i + jreturn(c(k, x[k]))
}
mind
函数为每行分别求最小值 (wmins
),再找出全局最小值,返回最小值,行号和列号
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))
}
示例
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,
byrow=T,
)print(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
print(
[1] 6 3 4
分析
增加记录当前行号的一列
nrow(q)
qq cbind(q, 1:n)print(qq)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0 12 13 8 20 1
[2,] 12 0 15 28 88 2
[3,] 13 15 0 6 9 3
[4,] 8 28 6 0 33 4
[5,] 20 88 9 33 0 5
为前 n-1 行,即前 4 行,寻找最小值,并记录最小值(第 2 行)及其对应的列数(第 1 行)。
返回结果中每列代表矩阵中的行,即第 1 列表示矩阵的第 1 行,最后一列表示矩阵的第 4 行
apply(qq[-n,], 1, imin)print(wmins)
[,1] [,2] [,3] [,4]
[1,] 4 3 4 5
[2,] 8 15 6 33
求最小值,返回列号,即矩阵中的行号 i
which.min(wmins[2,])print(i)
[1] 3
从 wmins
中取出矩阵的列号 j
wmins[1, i]print(j)
[1] 4
print(
[1] 6 3 4
如果最小值唯一,则有更简单的方法
注:有两次矩阵计算,效率不是最优的
function(d) {
smallest min(d)
ij which(d == smallest, arr.ind=TRUE)return(c(smallest, ij))
}
示例,需要将 0 值和重复值都替换为比较大的数
= q
[,1] [,2] [,3] [,4] [,5]
[1,] 99999 12 13 8 20
[2,] 99999 99999 15 28 88
[3,] 99999 99999 99999 6 9
[4,] 99999 99999 99999 99999 33
[5,] 99999 99999 99999 99999 99999
print(
[1] 6 3 4
向量与矩阵的差异
矩阵实际上就是向量,只不过多了两个属性:行数和列数
matrix(1:8, nrow=4)print(z)
[,1] [,2]
[1,] 1 5
[2,] 2 6
[3,] 3 7
[4,] 4 8
求长度
print(
[1] 8
但 z
不仅仅是向量
print(
[1] "matrix" "array"
print(
$dim
[1] 4 2
用 dim()
函数访问 dim 属性
print(
[1] 4 3
nrow()
和 ncol()
是 dim()
的简单封装
print(
[1] 4
print(
[1] 2
nrow
function (x)
dim(x)[1L]
避免意外降维
从矩阵中提取一行,会减少一维,即返回向量
matrix(1:8, nrow=4)print(z)
[,1] [,2]
[1,] 1 5
[2,] 2 6
[3,] 3 7
[4,] 4 8
z[2,]print(r)
[1] 2 6
print(
$dim
[1] 3 2
print(
NULL
str(z)
int [1:3, 1:2] 1 2 3 4 5 6
str(r)
int [1:2] 2 6
使用 drop
参数禁止自动减少维度
z[2,, drop=FALSE]print(r)
[,1] [,2]
[1,] 2 5
print(
[1] 1 2
可见中括号 [
实际上也是一个函数
print(z[3,
[1] 6
print(
[1] 6
as.matrix()
函数可以将向量转为矩阵
1:3print(u)
[1] 1 2 3
as.matrix(u)print(v)
[,1]
[1,] 1
[2,] 2
[3,] 3
print(
NULL
print(
$dim
[1] 3 1
矩阵的行和列的命名问题
矩阵的行和列也可以命名
matrix(1:4, nrow=2)print(z)
[,1] [,2]
[1,] 1 3
[2,] 2 4
print(
NULL
colnames(z)
a b
[1,] 1 3
[2,] 2 4
print(
[1] "a" "b"
print(z[,
[1] 1 2
高维数组
三个学生,两次考试,每次考试记录两个成绩
第一次考试的数据,行表示学生,列表示单次考试的两个成绩
matrix(c(46, 30,21, 25,50, 50
),
nrow=3,
byrow=T
)print(firsttest)
[,1] [,2]
[1,] 46 30
[2,] 21 25
[3,] 50 50
第二次考试的成绩
matrix(c(46, 43,41, 35,50, 50
),
nrow=3,
byrow=T
)print(secondtest)
[,1] [,2]
[1,] 46 43
[2,] 41 35
[3,] 50 50
合并为一个数组,两次考试形成第三维,叫做层 (layer)
array(
data=c(firsttest, secondtest),
dim=c(3, 2, 2)
)
print(
$dim
[1] 3 2 2
学生 3 (dim 1) 在第一次考试 (dim 3) 的第二个部分 (dim 2) 中的得分
print(tests[3,
[1] 50
print(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
参考
《学习 R 语言:快速开始》
《学习 R 语言:向量》
本文的 Jupyter Notebook 代码请访问如下项目:
https://github.com/perillaroc/the-art-of-r-programming
题图来自 Pixabay。