R的矩阵操作.docx

资源描述

R的矩阵操作.docx

《R的矩阵操作.docx》由会员分享，可在线阅读，更多相关《R的矩阵操作.docx（26页珍藏版）》请在冰豆网上搜索。

R的矩阵操作.docx

R的矩阵操作

1矩阵基本操作

1.1创建向量

R里面有多种方法来创建向量（Vector），最简单的是用函数c（）。

例如：

>X=c（1,2,3,4）

[1]1234

当然，还有别的方法。

例如：

>X=1:

[1]1234

还有seq（）函数。

例如：

>X=seq（1,4,length=4）

[1]1234

注意一点，R中的向量默认为列向量，如果要得到行向量需要对其进行转置。

1.2创建矩阵

R中创建矩阵的方法也有很多。

大致分为直接创建和由其它格式转换两种方法。

1.2.1直接创建矩阵

最简单的直接创建矩阵的方法是用matrix（）函数，matrix（）函数的使用方法如下：

>args（matrix）

function（data=NA,nrow=1,ncol=1,byrow=FALSE,dimnames=NULL）

NULL

其中，data参数输入的为矩阵的元素，不能为空；nrow参数输入的是矩阵的行数，默认为1；ncol参数输入的是矩阵的列数，默认为1；byrow参数控制矩阵元素的排列方式，TRUE表示按行排列，FALSE表示按列排列，默认为FALSE；dimnames参数输入矩阵的行名和列名，可以不输入，系统默认为NULL。

例如：

>matrix（1:

6,nrow=2,ncol=3,byrow=FALSE）

[,1] [,2] [,3]

[1,] 1 3 5

[2,] 2 4 6

改变矩阵的行数和列数：

>matrix（1:

6,nrow=3,ncol=2,byrow=FALSE）

[,1] [,2]

[1,] 1 4

[2,] 2 5

[3,] 3 6

改变byrow参数：

>matrix（1:

6,nrow=3,ncol=2,byrow=T）

[,1] [,2]

[1,] 1 2

[2,] 3 4

[3,] 5 6

设定矩阵的行名和列名：

>matrix（1:

6,nrow=3,ncol=2,byrow=T,dimnames=list（c（"A","B","C"）,c（"boy","girl"）））

boy girl

A 1 2

B 3 4

C 5 6

1.2.2由其它格式转换

也可以由其它格式的数据转换为矩阵，此时需要用到函数as.matrix（）。

1.3查看和改变矩阵的维数

矩阵有两个维数，即行维数和列维数。

在R中查看矩阵的行维数和列维数可以用函数dim（）。

例如：

>X=matrix（1:

12,ncol=3,nrow=4）

[,1][,2][,3]

[1,] 1 5 9

[2,] 2 6 10

[3,] 3 7 11

[4,] 4 8 12

>dim（X）

[1]43

只返回行维数：

>dim（X）[1]

[1]4

也可以用函数nrow（）

>nrow（X）

[1]4

只返回列维数：

>dim（X）[2]

[1]3

也可以用函数ncol（）:

>ncol（X）

[1]3

同时，函数dim（）也可以改变矩阵的维数。

例如：

>dim（X）=c（2,6）

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

[1,] 1 3 5 7 9 11

[2,] 2 4 6 8 10 12

1.4矩阵行列的名称

查看矩阵的行名和列名分别用函数rownames（）和函数colnames（）。

例如：

> X=matrix（1:

6,nrow=3,ncol=2,byrow=T,dimnames=list（c（"A","B","C"）,c（"boy","girl"）））

boygirl

A 1 2

B 3 4

C 5 6

查看矩阵的行名：

>rownames（X）

[1]"A""B""C"

查看矩阵的列名：

>colnames（X）

[1]"boy" "girl"

同时也可以改变矩阵的行名和列名，比如：

> rownames（X）=c（"E","F","G"）

boygirl

E 1 2

F 3 4

G 5 6

>colnames（X）=c（"man","woman"）

manwoman

E 1 2

F 3 4

G 5 6

1.5矩阵元素的查看及其重新赋值

查看矩阵的某个特定元素，只需要知道该元素的行坐标和列坐标即可，例如：

>X=matrix（1:

12,nrow=4,ncol=3）

[,1][,2][,3]

[1,] 1 5 9

[2,] 2 6 10

[3,] 3 7 11

[4,] 4 8 12

查看位于矩阵第二行、第三列的元素：

>X[2,3]

[1]10

也可以重新对矩阵的元素进行赋值，将矩阵第二行、第三列的元素替换为0：

>X[2,3]=0

[,1][,2][,3]

[1,] 1 5 9

[2,] 2 6 0

[3,] 3 7 11

[4,] 4 8 12

R中有一个diag（）函数可以返回矩阵的全部对角元素：

>X=matrix（1:

9,ncol=3,nrow=3）

[,1] [,2] [,3]

[1,] 1 4 7

[2,] 2 5 8

[3,] 3 6 9

>diag（X）

[1]159

当然也可以对对角元素进行重新赋值：

>diag（X）=c（0,0,1）

[,1][,2][,3]

[1,] 0 4 7

[2,] 2 0 8

[3,] 3 6 1

当操作对象不是对称矩阵时，diag（）也可以进行操作。

>X=matrix（1:

12,nrow=4,ncol=3）

[,1][,2][,3]

[1,] 1 5 9

[2,] 2 6 10

[3,] 3 7 11

[4,] 4 8 12

>diag（X）

[1] 1 6 11

diag（）函数还能用来生成对角矩阵：

>diag（c（1,2,3））

[,1][,2][,3]

[1,] 1 0 0

[2,] 0 2 0

[3,] 0 0 3

也可以生成单位对角矩阵：

>diag（3）

[,1][,2][,3]

[1,] 1 0 0

[2,] 0 1 0

[3,] 0 0 1

>diag（4）

[,1][,2][,3][,4]

[1,] 1 0 0 0

[2,] 0 1 0 0

[3,] 0 0 1 0

[4,] 0 0 0 1

查看矩阵的上三角部分：

在R中查看矩阵的上三角和下三角部分很简单。

可以通过lower.tri（）和upper.tir（）来实现：

>args（lower.tri）

function（x,diag=FALSE）

NULL

>args（upper.tri）

function（x,diag=FALSE）

NULL

查看上三角：

>X=matrix（1:

12,nrow=4,ncol=3）

[,1][,2][,3]

[1,] 1 5 9

[2,] 2 6 10

[3,] 3 7 11

[4,] 4 8 12

>X[lower.tri（X）]

[1] 2 3 4 7 812

改变赋值：

>X[lower.tri（X）]=0

[,1][,2][,3]

[1,] 1 5 9

[2,] 0 6 10

[3,] 0 0 11

[4,] 0 0 0

2矩阵计算

2.1矩阵转置

R中矩阵的转置可以用t（）函数完成，例如：

>X=matrix（1:

12,nrow=4,ncol=3）

[,1][,2][,3]

[1,] 1 5 9

[2,] 2 6 10

[3,] 3 7 11

[4,] 4 8 12

>t（X）

[,1][,2][,3][,4]

[1,] 1 2 3 4

[2,] 5 6 7 8

[3,] 9 10 11 12

2.2矩阵的行和与列和及行平均值和列均值

在R中很容易计算一个矩阵的各行和和各列和以及各行的平均值和各列的平均值。

例如：

>A=matrix（1:

12,3,4）

[,1][,2][,3][,4]

[1,] 1 4 7 10

[2,] 2 5 8 11

[3,] 3 6 9 12

>rowSums（A）

[1]222630

>rowMeans（A）

[1]5.56.57.5

>colSums（A）

[1] 6152433

>colMeans（A）

[1] 2 5 811

2.3行列式的值

R中的函数det（）将计算方阵A的行列式。

例如：

>X=matrix（rnorm（9）,nrow=3,ncol=3）

[,1] [,2] [,3]

[1,] 0.05810412-1.2992698 0.5630315

[2,]-0.28070583 0.1958623-1.8202283

[3,] 0.83691209 0.4411497 1.0014306

>det（X）

[1]1.510076

2.4矩阵相加减

矩阵元素的相加减是指维数相同的矩阵，处于同行和同列的位置的元素进行加减。

实现这个功能用“＋”，“－”即可。

例如：

>A=B=matrix（1:

12,nrow=3,ncol=4）

>A+B

[,1][,2][,3][,4]

[1,] 2 8 14 20

[2,] 4 10 16 22

[3,] 6 12 18 24

>A-B

[,1][,2][,3][,4]

[1,] 0 0 0 0

[2,] 0 0 0 0

[3,] 0 0 0 0

2.5矩阵的数乘

矩阵的数乘是指一个常数与一个矩阵相乘。

设A为m×n矩阵，c≠0，在R中求cA的值，可以用符号“*”。

例如：

>c=2

>A=matrix（1:

12,nrow=3,ncol=4）

[,1][,2][,3][,4]

[1,] 1 4 7 10

[2,] 2 5 8 11

[3,] 3 6 9 12

>c*A

[,1][,2][,3][,4]

[1,] 2 8 14 20

[2,] 4 10 16 22

[3,] 6 12 18 24

结果矩阵与原矩阵的所有相应元素都差一个常数c。

2.6矩阵相乘

2.6.1矩阵的乘法

A为m×n矩阵，B为n×k矩阵，在R中求AB，可以符号“%*%”。

例如：

>A=matrix（1:

12,nrow=3,ncol=4）

>B=matrix（1:

12,nrow=4,ncol=3）

>A%*%B

[,1][,2][,3]

[1,] 70 158 246

[2,] 80 184 288

[3,] 90 210 330

注意BA无意义，因其不符合矩阵的相乘规则。

若A为n×m矩阵，B为n×k矩阵，在R中求A’B：

>A=matrix（1:

12,nrow=4,ncol=3）

>B=matrix（1:

12,nrow=4,ncol=3）

>t（A）%*%B

[,1][,2][,3]

[1,] 30 70 110

[2,] 70 174 278

[3,] 110 278 446

也可以用函数crossprod（）计算A’B：

>crossprod（A,B）

[,1][,2][,3]

[1,] 30 70 110

[2,] 70 174 278

[3,] 110 278 446

2.6.2矩阵的Kronecker积

n×m矩阵A和h×k矩阵B的Kronecker积是一个nh×mk维矩阵，公式为：

a11B …a1nB

Am×n×Bh×k= … …

am1B …amnB mh×nk

在R中Kronecker积可以用函数kronecher（）来计算。

例如：

>A=matrix（1:

4,2,2）

[,1][,2]

[1,] 1 3

[2,] 2 4

>B=matrix（rep（1,4）,2,2）

[,1][,2]

[1,] 1 1

[2,] 1 1

>kronecker（A,B）

[,1][,2][,3][,4]

[1,] 1 1 3 3

[2,] 1 1 3 3

[3,] 2 2 4 4

[4,] 2 2 4 4

2.7矩阵的伴随矩阵

求矩阵A的伴随矩阵可以用LoopAnalyst包中的函数make.adjoint（）函数。

例如：

>install.packages（"LoopAnalyst"）

>A=matrix（1:

12,nrow=3,ncol=4）

[,1][,2][,3][,4]

[1,] 1 4 7 10

[2,] 2 5 8 11

[3,] 3 6 9 12

>make.adjoint（A）

[,1][,2][,3]

[1,] -3 6 -3

[2,] 6 -12 6

[3,] -3 6 -3

2.8矩阵的逆和广义逆

2.8.1矩阵的逆

矩阵A的逆A-1可以用函数solve（），例如：

>A=matrix（rnorm（9）,nrow=3,ncol=3）

[,1] [,2] [,3]

[1,]-0.2915845 0.2831544 0.94493154

[2,]-1.6494678 0.6999185-0.06292334

[3,]-0.7224015-0.3906971 0.44799963

>solve（A）

[,1] [,2] [,3]

[1,]0.2359821-0.4050650-0.5546321

[2,]0.6405592 0.4507583-1.2877720

[3,]0.9391490-0.2600663 0.2147417

验证AA-1=1：

>A%*%solve（A）

[,1] [,2] [,3]

[1,] 1.000000e+008.433738e-17-1.341700e-18

[2,] 1.216339e-171.000000e+00-4.667152e-17

[3,]-2.203641e-174.283954e-17 1.000000e+00

用round函数可以更好的得到结果：

>round（A%*%solve（A））

[,1][,2][,3]

[1,] 1 0 0

[2,] 0 1 0

[3,] 0 0 1

solve（）函数也可以用来求解方程组ax=b。

2.8.2矩阵的广义逆（Moore-Penrose）

并非所有的矩阵都有逆，但是所有的矩阵都可有广义逆。

n×m矩阵A+是矩阵A的Moore-Penrose逆，如果它满足下列条件：

AA+A=A

A+AA+=A+

（AA+）T=AA+

（A+A）T=A+A

R中MASS包中的ginv（）函数可以计算矩阵的Moore-Penrose逆。

例如：

>library（MASS）

>A=matrix（1:

12,nrow=3,ncol=4）

[,1][,2][,3][,4]

[1,] 1 4 7 10

[2,] 2 5 8 11

[3,] 3 6 9 12

>solve（A）

Errorinsolve.default（A）:

onlysquarematricescanbeinverted

>ginv（A）

[,1] [,2] [,3]

[1,]-0.483333333-0.03333333 0.41666667

[2,]-0.244444444-0.01111111 0.22222222

[3,]-0.005555556 0.01111111 0.02777778

[4,] 0.233333333 0.03333333-0.16666667

验证性质①：

>A%*%ginv（A）%*%A

[,1][,2][,3][,4]

[1,] 1 4 7 10

[2,] 2 5 8 11

[3,] 3 6 9 12

[,1][,2][,3][,4]

[1,] 1 4 7 10

[2,] 2 5 8 11

[3,] 3 6 9 12

验证性质②：

>ginv（A）%*%A%*%ginv（A）

[,1] [,2] [,3]

[1,]-0.483333333-0.03333333 0.41666667

[2,]-0.244444444-0.01111111 0.22222222

[3,]-0.005555556 0.01111111 0.02777778

[4,] 0.233333333 0.03333333-0.16666667

>ginv（A）

[,1] [,2] [,3]

[1,]-0.483333333-0.03333333 0.41666667

[2,]-0.244444444-0.01111111 0.22222222

[3,]-0.005555556 0.01111111 0.02777778

[4,] 0.233333333 0.03333333-0.16666667

验证性质③：

>A%*%ginv（A）

[,1] [,2] [,3]

[1,] 0.83333330.3333333-0.1666667

[2,] 0.33333330.3333333 0.3333333

[3,]-0.16666670.3333333 0.8333333

>t（A%*%ginv（A））

[,1] [,2] [,3]

[1,] 0.83333330.3333333-0.1666667

[2,] 0.33333330.3333333 0.3333333

[3,]-0.16666670.3333333 0.8333333

验证性质④：

>ginv（A）%*%A

[,1][,2][,3][,4]

[1,] 0.7 0.4 0.1-0.2

[2,] 0.4 0.3 0.2 0.1

[3,] 0.1 0.2 0.3 0.4

[4,]-0.2 0.1 0.4 0.7

>t（ginv（A）%*%A）

[,1][,2][,3][,4]

[1,] 0.7 0.4 0.1-0.2

[2,] 0.4 0.3 0.2 0.1

[3,] 0.1 0.2 0.3 0.4

[4,]-0.2 0.1 0.4 0.7

也可以不必如此麻烦来验证性质③和④，因为③和④只是表明AA+和A+A是对称矩阵。

2.8.3X’X的逆

很多时候，我们需要计算形如X’X的逆。

这很容易实现，例如：

>x=matrix（rnorm（9）,ncol=3,nrow=3）

[,1] [,2] [,3]

[1,]-0.1806586-0.763405120.002652331

[2,]-1.8018584 0.044679431.416332187

[3,] 1.2785359-1.316535130.180653002

>solve（crossprod（x））

[,1] [,2] [,3]

[1,]1.21818370.96645761.

展开阅读全文