R官方入门教程（5）数组和矩阵

数计加油站

已于 2023-04-22 21:44:25 修改

阅读量558

点赞数

分类专栏： R语言文章标签：矩阵 r语言算法

于 2023-04-19 22:22:46 首次发布

本文链接：https://blog.csdn.net/2301_76324977/article/details/130254888

版权

R语言专栏收录该内容

14 篇文章 5 订阅

订阅专栏

R是一门著名的可用于数据和统计分析的程序语言，本文翻译自R软件官方文档教程《An Introduction to R》，仅供学习和参考。

5 数组和矩阵

5.1 数组

数组array可以看成是多维向量。R支持使用比较简单的方式来创建和操作数组。特别地，矩阵matrix就是二维数组。

维数向量dimension vector是一个由正整数组成的向量。它的长度称为它对应的数组的维数。比如,c(4,3,2)对应着一个大小为 $4\times 3 \times 2$ 的3维数组。

一个向量被赋予了维数dim属性后才会被R当成数组。假设z是一个长度为1500的向量，

> dim(z) <- c(3,5,1000)

上面这行代码将赋予z一个dim属性，z也就可以被当成一个 $3\times 5 \times 1000$ 的数组。dim(array)用于查看或设置数组的维数向量。

在R中,与FORTRAN类似,数组是按列优先储存的。也就是，前面的下标变化快·，后面的下标变化慢。

举个例子，假设a的维数向量为c(3,2,2),那么在内部,R是以a[1,1,1]，a[2,1,1]，a[3,1,1]，a[1,2,1]，a[2,2,1],a[3,2,1],a[1,1,2]，a[2,1,2]，a[3,1,2]，a[1,2,2]，a[2,2,2],a[3,2,2]

> a<-1:12
> dim(a)<-c(3,2,2)
> a
, , 1

     [,1] [,2]
[1,]    1    4
[2,]    2    5
[3,]    3    6

, , 2

     [,1] [,2]
[1,]    7   10
[2,]    8   11
[3,]    9   12

数组也可以是1维的。这种情况下和向量没什么区别(打印方式都一样)。但一旦出现区别，就让人很迷惑了。

5.2 数组索引和子数组

数组的元素可以通过数组名[下标1,下标2,······,下标n]的方式进行索引。

如果下标都是一个正整数，就可以索引某个元素。

如果有的下标是一个整数向量，就可以索引对应的多个元素。

如果有的下标不写，代表这个维度全都索引。

继续之前的例子，a的维数向量为c(3,2,2),那么

a[2,1,1]代表(2,1,1)这个位置的元素
a[c(1,3),1,1]代表(1,1,1)和(3,1,1)位置元素组成的子数组。
a[2,,]代表(2,1,1),(2,2,1),(2,1,2),(2,2,2)位置元素组成的子数组

> a
, , 1

     [,1] [,2]
[1,]    1    4
[2,]    2    5
[3,]    3    6

, , 2

     [,1] [,2]
[1,]    7   10
[2,]    8   11
[3,]    9   12
> a[2,1,1]
[1] 2
> a[c(1,3),1,1]
[1] 1 3
> a[2,,]
     [,1] [,2]
[1,]    2    8
[2,]    5   11

5.3 索引矩阵

也可以用一个矩阵来索引数组的多个元素。索引矩阵的每一行代表一个位置。索引矩阵的列数应该与被索引数组的维数相同。

特别地，我们以矩阵为例。矩阵是2维数组。所以我们的索引矩阵必须有两列。

例如，对于一个矩阵X，我们现在有这样两个任务:

提取X[1,3],X[2,2],X[3,1]
设置1中对应位置为0

为了完成以上任务，我们需要一个 $3\times 2$ 的索引矩阵

> 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
> i <- array(c(1:3,3:1), dim=c(3,2))	# 创建索引矩阵
> i
     [,1] [,2]
[1,]    1    3
[2,]    2    2
[3,]    3    1
> X[i]	# 提取X[1,3],X[2,2],X[3,1]
[1] 9 6 3
> X[i] <- 0	# 设置1对应位置为0
> X
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    5    0   13   17
[2,]    2    0   10   14   18
[3,]    0    7   11   15   19
[4,]    4    8   12   16   20

索引矩阵中不允许出现负整数。但NA和0被允许:有NA的行在结果中也是NA,有0的行将被自动忽略。

再举一个比较实用的例子。假设我们现在有性别因子chinese,数学成绩等级因子math,我们可以将这两列组成的长表拓展为宽表。

# 创建性别因子
> gender <- factor(c("男","女","女","男","男","女","男","女"))
> gender
[1] 男 女 女 男 男 女 男 女
Levels: 男 女
> math <- factor(c("及格","优秀","良好","优秀","及格","不及格","良好","不及格"))	# 创建数学成绩等级因子
> math
[1] 及格   优秀   良好   优秀   及格   不及格 良好   不及格
Levels: 不及格 及格 良好 优秀

# 因子类别个数
> g <- length(levels(gender))	
> m <- length(levels(math))
> n <- length(gender)

# 创建全0矩阵
> Xg <- matrix(0, n, g)	
> Xg
     [,1] [,2]
[1,]    0    0
[2,]    0    0
[3,]    0    0
[4,]    0    0
[5,]    0    0
[6,]    0    0
[7,]    0    0
[8,]    0    0
> Xm <- matrix(0, n, m)
> Xm
     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0
[4,]    0    0    0    0
[5,]    0    0    0    0
[6,]    0    0    0    0
[7,]    0    0    0    0
[8,]    0    0    0    0

# 创建索引矩阵，行代表第几个样本，列是其类别
> ig <- cbind(1:n, gender)	
> ig
       gender
[1,] 1      1
[2,] 2      2
[3,] 3      2
[4,] 4      1
[5,] 5      1
[6,] 6      2
[7,] 7      1
[8,] 8      2
> im <- cbind(1:n, math)
> im
       math
[1,] 1    2
[2,] 2    4
[3,] 3    3
[4,] 4    4
[5,] 5    2
[6,] 6    1
[7,] 7    3
[8,] 8    1

# 设置对应位置为1
> Xg[ig] <- 1
> Xm[im] <- 1

# 分开的宽表
> Xg
     [,1] [,2]
[1,]    1    0
[2,]    0    1
[3,]    0    1
[4,]    1    0
[5,]    1    0
[6,]    0    1
[7,]    1    0
[8,]    0    1
> Xm
     [,1] [,2] [,3] [,4]
[1,]    0    1    0    0
[2,]    0    0    0    1
[3,]    0    0    1    0
[4,]    0    0    0    1
[5,]    0    1    0    0
[6,]    1    0    0    0
[7,]    0    0    1    0
[8,]    1    0    0    0

# 合并后的宽表，前两列代表性别，后两列代表数学成绩等级
> X <- cbind(Xg, Xm)
> X
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    0    0    1    0    0
[2,]    0    1    0    0    0    1
[3,]    0    1    0    0    1    0
[4,]    1    0    0    0    0    1
[5,]    1    0    0    1    0    0
[6,]    0    1    1    0    0    0
[7,]    1    0    0    0    1    0
[8,]    0    1    1    0    0    0

我们还可以建立性别因子与数学成绩等级因子的联列表：

> N <- crossprod(Xg, Xm)
> N
     [,1] [,2] [,3] [,4]
[1,]    0    2    1    1
[2,]    2    0    1    1

当然，用table()函数更加方便:

> M <- table(gender, math)
> M
      math
gender 不及格 及格 良好 优秀
    男      0    2    1    1
    女      2    0    1    1

5.4 array()函数

除了赋予向量一个dim属性外，数组还可以通过array()函数由向量构建，其形式如下：

 Z <- array(data_vector, dim_vector)

例如，如果向量h包含24个或更少的数字，则命令

 Z <- array(h, dim=c(3,4,2))

将把h的数据重新设置到 $3\times 4\times 2$ 的数组Z中。如果h的长度正好是24，则结果与

> Z <- h
> dim(Z) <- c(3,4,2)

相同。然而，如果h的长度小于24，则会利用循环补齐机制，使其长度达到24。而直接使用dim（h）<-c（3，4，2）将报出长度不匹配的错误。下面列举一个极端但常见的例子：

> Z <- array(0, 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\times 4\times 2$ 的全0矩阵。

这时，dim(Z)代表维度向量c(3，4，2)，Z[1:24]代表h中的数据向量，带有空下标的Z[]或没有下标的Z代表整个数组。

数组可以在算术表达式中使用，结果是由数据向量上的逐元素运算形成的数组。参与运算的数组的dim属性通常需要相同，这将成为结果的维度向量。因此，如果A、B和C都是相似的数组，那么

> D <- 2*A*B + C + 1

使得D是一个类似的数组，其数据向量是相关的数组的逐元素运算的结果。然而，必须更加仔细地考虑关于数组和向量混合计算的精确规则。

5.4.1 数组和向量混合计算

影响向量和数组逐元素混合计算的精确规则有些古怪，在参考文献中很难找到。根据经验，我们发现以下内容至少是可靠：

表达式从左到右进行计算。
任何比矩阵或数组的长度还要长的向量参与运算都会产生错误。
任何短向量都会利用循环补齐机制来进行扩展，直到它们与其他进行操作的数组或向量的大小相匹配。
只要遇到数组参与运算，那么这些数组都必须具有相同的dim属性，否则将产生错误。
如果表达式中存在数组，并且没有与向量相关的错误或强制转化，则结果是和数组拥有一样dim属性的结果数组。

5.5 两个数组的外积

数组上的一个重要操作是外积。如果a和b是两个数值数组，则它们的外积也是一个数值数组，其维度向量是通过连接它们的两个维度向量来形成的（顺序很重要），并且其数据向量是通过将a的数据向量的元素与b的数据向量形成所有可能的积来获得的。外积由特殊运算符%o%形成：

> ab <- a %o% b

用函数来进行等价表达为:

> ab <- outer(a, b, "*")

乘法"*"可以用接受两个变量的任意函数来代替。例如，如果我们想在由向量x和y定义的坐标网格上计算函数 $f(x,y)=\frac{cos(y)}{1+x^2}$ 的值，我们可以这样：

> f <- function(x, y) cos(y)/(1 + x^2)	# 定义一个函数
> z <- outer(x, y, f)

特别地，两个向量的外积是一个二维数组（一个秩最多为1的矩阵）。请注意，外积运算符当然是不满足交换律的。

一个例子：2乘2的个位数矩阵的行列式

作为一个没啥用但很好玩的例子，让我们来考虑 $2\times 2$ 矩阵 $\begin{bmatrix}a & b\\c&d\end{bmatrix}$ 的行列式，其中 $a, b, c, d$ 都是 $0,1,2,\cdots,9$ 中的某个数

问题是找到这种形式的所有可能矩阵的行列式 $a d - b c$ ，并将每个行列式出现的频率表示为高密度图。如果每个值都是独立一致随机选择的，那么这相当于找到行列式的概率分布。

一种巧妙的方法是使用outer()函数两次：

> d <- outer(0:9, 0:9) # 得出个位数的所有可能乘积ad或bc
> d
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    0    0    0    0    0    0    0    0    0     0
 [2,]    0    1    2    3    4    5    6    7    8     9
 [3,]    0    2    4    6    8   10   12   14   16    18
 [4,]    0    3    6    9   12   15   18   21   24    27
 [5,]    0    4    8   12   16   20   24   28   32    36
 [6,]    0    5   10   15   20   25   30   35   40    45
 [7,]    0    6   12   18   24   30   36   42   48    54
 [8,]    0    7   14   21   28   35   42   49   56    63
 [9,]    0    8   16   24   32   40   48   56   64    72
[10,]    0    9   18   27   36   45   54   63   72    81
>  fr <- table(outer(d, d, "-"))	# 得出所有可能的ad-bc，并转化为频数表
> fr

-81 -80 -79 -78 -77 -76 -75 -74 -73 -72 -71 -70 -69 -68 -67 -66 -65 -64 -63 -62 
 19   1   2   2   3   2   4   2   4  41   4   4   8   6   6  10   7  27  49   8 
-61 -60 -59 -58 -57 -56 -55 -54 -53 -52 -51 -50 -49 -48 -47 -46 -45 -44 -43 -42 
  8  17   8  12  18  53  13  60  12  18  22  16  35  70  22  24  66  28  18  72 
-41 -40 -39 -38 -37 -36 -35 -34 -33 -32 -31 -30 -29 -28 -27 -26 -25 -24 -23 -22 
 22  75  37  34  26 111  63  36  45  84  34  94  36  93  97  50  53 156  42  60 
-21 -20 -19 -18 -17 -16 -15 -14 -13 -12 -11 -10  -9  -8  -7  -6  -5  -4  -3  -2 
103 107  50 168  51 140 112 116  59 191  65 126 156 185 115 206 117 179 153 156 
 -1   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18 
111 570 111 156 153 179 117 206 115 185 156 126  65 191  59 116 112 140  51 168 
 19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38 
 50 107 103  60  42 156  53  50  97  93  36  94  34  84  45  36  63 111  26  34 
 39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58 
 37  75  22  72  18  28  66  24  22  70  35  16  22  18  12  60  13  53  18  12 
 59  60  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78 
  8  17   8   8  49  27   7  10   6   6   8   4   4  41   4   2   4   2   3   2 
 79  80  81 
  2   1  19 
> plot(fr, xlab="行列式", ylab="频率") # 绘制频率图

在这里插入图片描述

5.6 数组的广义转置

函数aperm(a，perm)可用于置换数组a。参数perm必须是整数 $\{1,2,\cdots,k\}$ 的一个排列，其中 $k$ 是 $a$ 的维数。函数的结果是一个与a大小相同的数组，但由perm[j]给出的旧维度成为新的第j维度。将此运算视为矩阵转置的推广是最简单的理解方法。事实上，如果A是一个矩阵（也就是说，一个双下标数组），那么

> B <- aperm(A, c(2,1))

就是A的转置。对于这种特殊情况，可以使用更简单的函数t()，所以我们可以使用B<-t(A)。

5.7 矩阵相关

如上所述，矩阵就是一个二维数组。然而，这是一个相当重要的特例，需要我们单独进行讨论。R包含许多仅对矩阵可用的运算符和函数。例如，如上所述，t(X)是矩阵转置函数。

函数nrow(A)和ncol(A)分别给出矩阵A中的行数和列数。

5.7.1 矩阵运算

运算符%*%用于矩阵乘法。如果在上下文中合适的话， $n\times1$ 或 $1\times n$ 矩阵当然可以当做n维向量。相反，如果可能的话，出现在矩阵乘法表达式中的向量会自动提升为行或列向量，只要某一种是正确的（尽管这并不总是明确可能的，正如我们稍后看到的那样）。

例如，如果A和B是大小相同的方阵，则A*B是逐元素乘积得到的矩阵，而A%*%B是矩阵乘法得到的矩阵。如果x是一个向量，那么x %*% a %*% x是一个二次型。

函数crossprod()形成了“交叉乘积”，这意味着crossprod(X,y)与t(X) %*% y相同，但操作更有效。如果省略了crossprod()的第二个参数，它将被视为与第一个参数相同。

diag()函数的作用取决于它的参数。diag(v)，其中v是一个向量，就会给出一个对角矩阵，其中v的元素作为对角项。另一方面，diag(M)，其中M是矩阵，给出M的主对角线元素组成的的向量。这与Matlab中用于diag()的约定相同。此外，有点令人困惑的是，如果k是一个单一的数值，那么diag(k)就是k乘k的单位矩阵！

5.7.2 线性方程组和矩阵的逆

求解线性方程组的实质是矩阵乘法的逆运算。当b<- A %*% x之后，只需要给定A和b时，向量x就是对应线性方程组的解。

在R中，solve(A,b)将返回解x（有一定的精度损失）。

注意，在线性代数中，形式上有 $x=A^{−1}b$ ，其中 $A^{−1}$ 表示A的逆矩阵，可以通过solve(A)来计算。

不过，从数值上考虑，x<- solve(A) %*% b与slove(A,b)相比，效率低，而且更加不稳定。

与此类似，多元计算中使用的二次型 $x^{T}A^{-1}x$ 也应该通过x %*% solve(A,x)来计算，而不是直接计算A的逆。

5.7.3 特征值和特征向量

函数eigen(Sm)能够计算对称矩阵Sm的特征值和特征向量。该函数的结果是具有两个分量的列表，分别命名为values(特征值)和vectors(特征向量)。

> ev <- eigen(Sm)

上述代码将Sm的特征列表赋给ev变量。通过ev$val可以获取特征值，ev$c=vec可以获取特征向量。

如果我们只需要特征值，我们可以直接这样:

> evals <- eigen(Sm)$values

现在evals只包含特征值，而特征向量则被直接抛弃了。

如果直接用表达式

> eigen(Sm)

那么特征值和特征向量都会被直接打印，而没有储存起来。

对于大型矩阵，如果并不需要获取特征向量，则最好避免计算这些特征向量:

evals <- eigen(Sm, only.values = TRUE)$values

5.7.4 奇异值分解与行列式

函数svd(M)可以计算矩阵M的奇异值分解。它由 $U, D, V$ 三个矩阵组成，满足M＝U %*% D %*% t(V)。D实际上是以对角线上元素组成的向量形式返回的。svd(M)的结果实际上是一个名为d、u和v的三个元素组成的列表，具有很高的现实价值。

如果M是个方阵，那么不难看出absdetM <- prod(svd(M)$d)可以计算M的行列式的绝对值。

如果经常需要对各种矩阵进行这种计算，则可以将其自定义为一个函数

> absdet <- function(M) prod(svd(M)$d)

之后,我们就可以直接使用absdet()。

作为另一个琐碎但可能有用的例子，您可能会考虑编写一个函数，比如tr()，来计算方阵的迹。

R有一个内置函数det()用来计算行列式的值。

5.7.5 最小二乘拟合与QR分解

函数lsfit()返回一个列表，给出最小二乘法拟合的结果。

> ans <- lsfit(X, y)

上述代码给出了最小二乘拟合的结果，其中y是观测向量，X是设计矩阵。有关更多详细信息，请参阅帮助工具，也可以参阅后续函数ls.diag()，其中包括回归诊断。请注意，大平均项是自动包含的，不需要显式地作为X的列。此外，请注意，在之后，您几乎总是更喜欢使用lm()而不是lsfit()来处理线性模型。

另一个有用的的函数是qr()及其相关函数。

> Xplus <- qr(X)
> b <- qr.coef(Xplus, y)
> fit <- qr.fitted(Xplus, y)
> res <- qr.resid(Xplus, y)

以上代码将计算y在拟合的X范围上的正交投影，在res中的正交补码上的投影，以及在b中的投影的系数向量，也就是说，b本质上是Matlab中/运算符的结果。

不需要假设X列满秩。冗余将被发现并删除。

这种替代方法是执行最小二乘法计算的较老、较低级别的方法。尽管在某些情况下仍然有用，但现在通常会被统计模型所取代，这将在第11章中进行讨论。、

形成分块矩阵——cbind()和rbind()

正如我们已经非正式地看到的，可以通过函数cbind()将矩阵水平或按列绑定在一起，通过函数rbind()将矩阵垂直或按行绑定在一起来形成新的矩阵。

X<-cbind（arg_1,arg_2,arg_3,...）

cbind()的参数必须向量或者具有相同行数的矩阵。结果是一个从左到右包含串联参数矩阵arg 1、arg 2，…的矩阵。

如果cbind()的一些参数是向量，则它们可能比存在的任何矩阵的列长度都短，在这种情况下，它们被循环扩展以匹配矩阵列大小（或者，如果没有给定矩阵，则为最长向量的长度）。

函数rbind()对行执行类似的操作。在这种情况下，任何向量参数，可能是循环扩展的，当然都被视为行向量。

假设X1和X2具有相同的行数。要将这两个矩阵按列组合成矩阵X，再加上为1的一个初始列，我们可以：

> X<-cbind(1，X1，X2)

rbind()或cbind()的结果始终具有矩阵状态。因此，cbind(x)和rbind(x)可能是明确允许向量x分别被视为列矩阵或行矩阵的最简单方法。

5.8 c()与矩阵

需要注意的是,rbind()和cbind()处理的是对应维度的矩阵连接问题。然而，基本的c()函数却不是这样。它清除掉矩阵的所有维度属性，再将其按向量逐元素连接。这偶尔很有用。

将数组强制返回到简单向量对象的官方方法是使用as.vector()函数

> vec <- as.vector(X)

然而，只需将X传入c()，就可以获得类似的结果:

> vec <- c(X)

两者之间略有不同，但最终两者之间的选择主要是风格问题（前者更可取）。

5.9 因子频率表

回想一下，一个因子可以将一组数据划分为多个组。类似地，一对因子可以划分双向交叉组，依此类推。函数table()允许根据等长因子计算频率表。如果有 $k$ 个因子变量，结果将是一个 $k$ 向频率数组。

假设grade是一个成绩因子factor(c("优","良","差","优","良","差","优","良","差","优"))，那么

> gradefr <- table(grade)

将返回一个频数表，以grade的Levels作为列名：

> gradefr
grade
差 良 优 
 3  3  4

这相当于，但简单于:

> gradefr <- tapply(grade,grade,length)

更进一步，考虑性别因子gender <- factor(c("男","女","男","男","女","女","男","女","男","女")，形成双向频率表。

> gradefr <- table(grade,gender)
> gradefr
     gender
grade 男 女
   差  2  1
   良  0  3
   优  3  1