加载中…
个人资料
微博
  • 博客等级：
  • 博客积分：
  • 博客访问：
  • 关注人气：
  • 获赠金笔：0支
  • 赠出金笔：0支
  • 荣誉徽章：
正文 字体大小：大 中 小

R语言矩阵运算(转2)

(2011-10-08 13:36:18)
标签：

it

 分类： R
16   矩阵的行和、列和、行平均与列平均
  在R中很容易求得一个矩阵的各行的和、平均数与列的和、平均数，例如：
  > A
  [,1] [,2] [,3] [,4]
[1,]   1   4   7   10
[2,]   2   5   8   11
[3,]   3   6   9   12
> rowSums(A)
[1] 22 26 30
> rowMeans(A)
[1] 5.5 6.5 7.5
> colSums(A)
[1] 6 15 24 33
> colMeans(A)
[1] 2 5 8 11
上述关于矩阵行和列的操作，还可以使用apply()函数实现。
> args(apply)
function (X, MARGIN, FUN, ...)
其中：x为矩阵，MARGIN用来指定是对行运算还是对列运算，MARGIN＝1表示对行运算，MARGIN＝2表示对列运算，FUN用来指定运算函数, ...用来给定FUN中需要的其它的参数，例如：
> apply(A,1,sum)
[1] 22 26 30
> apply(A,1,mean)
[1] 5.5 6.5 7.5
> apply(A,2,sum)
[1] 6 15 24 33
> apply(A,2,mean)
[1] 2 5 8 11
apply()函数功能强大，我们可以对矩阵的行或者列进行其它运算，例如：
计算每一列的方差
> A=matrix(rnorm(100),20,5)
> apply(A,2,var)
[1] 0.4641787 1.4331070 0.3186012 1.3042711 0.5238485
> apply(A,2,function(x,a)x*a,a=2)
  [,1] [,2] [,3] [,4]
[1,]   2   8   14   20
[2,]   4   10   16   22
[3,]   6   12   18   24
注意：apply(A,2,function(x,a)x*a,a=2)与A*2效果相同，此处旨在说明如何应用alpply函数。
17   矩阵X'X的逆
  在统计计算中，我们常常需要计算这样矩阵的逆，如OLS估计中求系数矩阵。R中的包“strucchange”提供了有效的计算方法。
  > args(solveCrossprod)
function (X, method = c("qr", "chol", "solve"))
其中：method指定求逆方法，选用“qr”效率最高，选用“chol”精度最高，选用“slove”与slove(crossprod(x,x))效果相同，例如：
> A=matrix(rnorm(16),4,4)
> solveCrossprod(A,method="qr")
      [,1]     [,2]     [,3]     [,4]
[1,] 0.6132102 -0.1543924 -0.2900796 0.2054730
[2,] -0.1543924 0.4779277 0.1859490 -0.2097302
[3,] -0.2900796 0.1859490 0.6931232 -0.3162961
[4,] 0.2054730 -0.2097302 -0.3162961 0.3447627
> solveCrossprod(A,method="chol")
      [,1]     [,2]     [,3]     [,4]
[1,] 0.6132102 -0.1543924 -0.2900796 0.2054730
[2,] -0.1543924 0.4779277 0.1859490 -0.2097302
[3,] -0.2900796 0.1859490 0.6931232 -0.3162961
[4,] 0.2054730 -0.2097302 -0.3162961 0.3447627
> solveCrossprod(A,method="solve")
      [,1]     [,2]     [,3]     [,4]
[1,] 0.6132102 -0.1543924 -0.2900796 0.2054730
[2,] -0.1543924 0.4779277 0.1859490 -0.2097302
[3,] -0.2900796 0.1859490 0.6931232 -0.3162961
[4,] 0.2054730 -0.2097302 -0.3162961 0.3447627
> solve(crossprod(A,A))
      [,1]     [,2]     [,3]     [,4]
[1,] 0.6132102 -0.1543924 -0.2900796 0.2054730
[2,] -0.1543924 0.4779277 0.1859490 -0.2097302
[3,] -0.2900796 0.1859490 0.6931232 -0.3162961
[4,] 0.2054730 -0.2097302 -0.3162961 0.3447627
18   取矩阵的上、下三角部分
  在R中，我们可以很方便的取到一个矩阵的上、下三角部分的元素，函数lower.tri()和函数upper.tri()提供了有效的方法。
  > args(lower.tri)
function (x, diag = FALSE)
函数将返回一个逻辑值矩阵，其中下三角部分为真，上三角部分为假，选项diag为真时包含对角元素，为假时不包含对角元素。upper.tri()的效果与之孑然相反。例如：
> A
  [,1] [,2] [,3] [,4]
[1,]   1   5   9   13
[2,]   2   6   10   14
[3,]   3   7   11   15
[4,]   4   8   12   16
> lower.tri(A)
    [,1] [,2] [,3] [,4]
[1,] FALSE FALSE FALSE FALSE
[2,] TRUE FALSE FALSE FALSE
[3,] TRUE TRUE FALSE FALSE
[4,] TRUE TRUE TRUE FALSE
> lower.tri(A,diag=T)
  [,1] [,2] [,3] [,4]
[1,] TRUE FALSE FALSE FALSE
[2,] TRUE TRUE FALSE FALSE
[3,] TRUE TRUE TRUE FALSE
[4,] TRUE TRUE TRUE TRUE
> upper.tri(A)
    [,1] [,2] [,3] [,4]
[1,] FALSE TRUE TRUE TRUE
[2,] FALSE FALSE TRUE TRUE
[3,] FALSE FALSE FALSE TRUE
[4,] FALSE FALSE FALSE FALSE
> upper.tri(A,diag=T)
    [,1] [,2] [,3] [,4]
[1,] TRUE TRUE TRUE TRUE
[2,] FALSE TRUE TRUE TRUE
[3,] FALSE FALSE TRUE TRUE
[4,] FALSE FALSE FALSE TRUE
> A[lower.tri(A)]=0
> A
  [,1] [,2] [,3] [,4]
[1,]   1   5   9   13
[2,]   0   6   10   14
[3,]   0   0   11   15
[4,]   0   0   0   16
> A[upper.tri(A)]=0
> A
  [,1] [,2] [,3] [,4]
[1,]   1   0   0   0
[2,]   2   6   0   0
[3,]   3   7   11   0
[4,]   4   8   12   16
19   backsolve&fowardsolve函数
这两个函数用于解特殊线性方程组，其特殊之处在于系数矩阵为上或下三角。
> args(backsolve)
function (r, x, k = ncol(r), upper.tri = TRUE, transpose = FALSE)
> args(forwardsolve)
function (l, x, k = ncol(l), upper.tri = FALSE, transpose = FALSE)
其中：r或者l为n×n维三角矩阵，x为n×1维向量，对给定不同的upper.tri和transpose的值，方程的形式不同
对于函数backsolve()而言，
例如：
  > A=matrix(1:9,3,3)
> A
  [,1] [,2] [,3]
[1,]   1   4   7
[2,]   2   5   8
[3,]   3   6   9
> x=c(1,2,3)
> x
[1] 1 2 3
> B=A
> B[upper.tri(B)]=0
> B
  [,1] [,2] [,3]
[1,]   1   0   0
[2,]   2   5   0
[3,]   3   6   9
> C=A
> C[lower.tri(C)]=0
> C
  [,1] [,2] [,3]
[1,]   1   4   7
[2,]   0   5   8
[3,]   0   0   9
> backsolve(A,x,upper.tri=T,transpose=T)
[1] 1.00000000 -0.40000000 -0.08888889
> solve(t(C),x)
[1] 1.00000000 -0.40000000 -0.08888889
> backsolve(A,x,upper.tri=T,transpose=F)
[1] -0.8000000 -0.1333333 0.3333333
> solve(C,x)
[1] -0.8000000 -0.1333333 0.3333333
> backsolve(A,x,upper.tri=F,transpose=T)
[1] 1.111307e-17 2.220446e-17 3.333333e-01
> solve(t(B),x)
[1] 1.110223e-17 2.220446e-17 3.333333e-01
> backsolve(A,x,upper.tri=F,transpose=F)
[1] 1 0 0
> solve(B,x)
[1] 1.000000e+00 -1.540744e-33 -1.850372e-17
对于函数forwardsolve()而言，
例如：
  > A
      [,1] [,2] [,3]
[1,]   1   4   7
[2,]   2   5   8
[3,]   3   6   9
> B
  [,1] [,2] [,3]
[1,]   1   0   0
[2,]   2   5   0
[3,]   3   6   9
> C
  [,1] [,2] [,3]
[1,]   1   4   7
[2,]   0   5   8
[3,]   0   0   9
> x
[1] 1 2 3
> forwardsolve(A,x,upper.tri=T,transpose=T)
[1] 1.00000000 -0.40000000 -0.08888889
> solve(t(C),x)
[1] 1.00000000 -0.40000000 -0.08888889
> forwardsolve(A,x,upper.tri=T,transpose=F)
[1] -0.8000000 -0.1333333 0.3333333
> solve(C,x)
[1] -0.8000000 -0.1333333 0.3333333
> forwardsolve(A,x,upper.tri=F,transpose=T)
[1] 1.111307e-17 2.220446e-17 3.333333e-01
> solve(t(B),x)
[1] 1.110223e-17 2.220446e-17 3.333333e-01
> forwardsolve(A,x,upper.tri=F,transpose=F)
[1] 1 0 0
> solve(B,x)
[1] 1.000000e+00 -1.540744e-33 -1.850372e-17
20   row()与col()函数
在R中定义了的这两个函数用于取矩阵元素的行或列下标矩阵，例如矩阵A={aij}m×n，
row()函数将返回一个与矩阵A有相同维数的矩阵，该矩阵的第i行第j列元素为i，函数col()类似。例如：
> x=matrix(1:12,3,4)
> row(x)
  [,1] [,2] [,3] [,4]
[1,]   1   1   1   1
[2,]   2   2   2   2
[3,]   3   3   3   3
> col(x)
  [,1] [,2] [,3] [,4]
[1,]   1   2   3   4
[2,]   1   2   3   4
[3,]   1   2   3   4
这两个函数同样可以用于取一个矩阵的上下三角矩阵，例如：
> x
  [,1] [,2] [,3] [,4]
[1,]   1   4   7   10
[2,]   2   5   8   11
[3,]   3   6   9   12
> x[row(x)<col(x)]=0
> x
  [,1] [,2] [,3] [,4]
[1,]   1   0   0   0
[2,]   2   5   0   0
[3,]   3   6   9   0
> x=matrix(1:12,3,4)
> x[row(x)>col(x)]=0
> x
  [,1] [,2] [,3] [,4]
[1,]   1   4   7   10
[2,]   0   5   8   11
[3,]   0   0   9   12
21   行列式的值
在R中，函数det(x)将计算方阵x的行列式的值，例如：
> x=matrix(rnorm(16),4,4)
> x
      [,1]     [,2]     [,3]     [,4]
[1,] -1.0736375 0.2809563 -1.5796854 0.51810378
[2,] -1.6229898 -0.4175977 1.2038194 -0.06394986
[3,] -0.3989073 -0.8368334 -0.6374909 -0.23657088
[4,] 1.9413061 0.8338065 -1.5877162 -1.30568465
> det(x)
[1] 5.717667
22向量化算子
在R中可以很容易的实现向量化算子，例如：
vec<-function (x){
          t(t(as.vector(x)))
}
vech<-function (x){
          t(x[lower.tri(x,diag=T)])
}
> x=matrix(1:12,3,4)
> x
  [,1] [,2] [,3] [,4]
[1,]   1   4   7   10
[2,]   2   5   8   11
[3,]   3   6   9   12
> vec(x)
    [,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
> vech(x)
  [,1] [,2] [,3] [,4] [,5] [,6]
[1,]   1   2   3   5   6   9
23   时间序列的滞后值
  在时间序列分析中，我们常常要用到一个序列的滞后序列，R中的包“fMultivar”中的函数tslag()提供了这个功能。
  > args(tslag)
function (x, k = 1, trim = FALSE)
其中：x为一个向量，k指定滞后阶数，可以是一个自然数列，若trim为假，则返回序列与原序列长度相同，但含有NA值；若trim项为真，则返回序列中不含有NA值，例如：
> x=1:20
> tslag(x,1:4,trim=F)
    [,1] [,2] [,3] [,4]
[1,]   NA   NA   NA   NA
[2,]   1   NA   NA   NA
[3,]   2   1   NA   NA
[4,]   3   2   1   NA
[5,]   4   3   2   1
[6,]   5   4   3   2
[7,]   6   5   4   3
[8,]   7   6   5   4
[9,]   8   7   6   5
[10,]   9   8   7   6
[11,]   10   9   8   7
[12,]   11   10   9   8
[13,]   12   11   10   9
[14,]   13   12   11   10
[15,]   14   13   12   11
[16,]   15   14   13   12
[17,]   16   15   14   13
[18,]   17   16   15   14
[19,]   18   17   16   15
[20,]   19   18   17   16
> tslag(x,1:4,trim=T)
    [,1] [,2] [,3] [,4]
[1,]   4   3   2   1
[2,]   5   4   3   2
[3,]   6   5   4   3
[4,]   7   6   5   4
[5,]   8   7   6   5
[6,]   9   8   7   6
[7,]   10   9   8   7
[8,]   11   10   9   8
[9,]   12   11   10   9
[10,]   13   12   11   10
[11,]   14   13   12   11
[12,]   15   14   13   12
[13,]   16   15   14   13
[14,]   17   16   15   14
[15,]   18   17   16   15
[16,]   19   18   17   16

喜欢

0

阅读┊ 收藏 ┊ 喜欢▼ ┊打印┊举报/Report
  

新浪BLOG意见反馈留言板　欢迎批评指正

新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 产品答疑

Copyright © 1996 - 2022 SINA Corporation, All Rights Reserved

新浪公司 版权所有