R语言编程艺术#02#矩阵（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时，需要的矩阵如下：
> 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
>
