Calyx.biz
  • Home | 主页
  • Study | 学习
    • CPA | 注册会计师
    • Deep Learning | 深度学习
    • R/Python/Vim | 编程 >
      • Java Learning | Java 学习
      • R Tips | R技巧
      • Python Skill | Python技术
      • Masterpiece | 作品
    • Big Data | 大数据 >
      • Bioinformatics | 生物信息
      • Data Science | 数据科学
    • Language | 语言 >
      • Cantonese | 粤语
      • English | 英语
      • 日本語 | 日语
    • Citation | 他山之石
    • Application | 运用
    • Dota2Win | 刀塔制胜
    • Others | 别的 >
      • CIpipe
      • How to LOVE
      • FMSofChangsha.2005.
  • Fun | 娱乐
    • Billboard | 公告版
    • Essay | 随想
    • Album Image | 专辑封面
    • Movie Poster | 电影海报
    • Exercise | 锻炼
    • Funny | 有趣
    • Gastronome | 美食家
    • East Era | 东方时代 >
      • Masterpieces | 美味
      • Gastronomy | 美食主义
    • Downloading | 下载
  • Record | 记录
    • Milestone | 里程碑
    • Discussion | 讨论
    • About Us | 关于我们
    • LinkedPaper | 原文

R Matrix Operation | R 矩阵运算 

7/22/2013

0 Comments

 
原文:  http://blog.csdn.net/hawksoft/article/details/7773323

创建矩阵向量;矩阵加减,乘积;矩阵的逆;行列式的值;特征值与特征向量;QR分解;奇异值分解;广义逆;backsolve与fowardsolve函数;取矩阵的上下三角元素;向量化算子等.

1   创建一个向量
在R中可以用函数c()来创建一个向量,例如:
> x=c(1,2,3,4)
> x
[1] 1 2 3 4 

2   创建一个矩阵
在R中可以用函数matrix()来创建一个矩阵,应用该函数时需要输入必要的参数值。
> args(matrix)
function (data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL) 
data项为必要的矩阵元素,nrow为行数,ncol为列数,注意nrow与ncol的乘积应为矩阵元素个数,byrow项控制排列元素时是否按行进行,dimnames给定行和列的名称。例如:
> 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
> 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
> matrix(1:12,nrow=4,ncol=3,byrow=T)
    [,1] [,2] [,3]
[1,]   1   2   3
[2,]   4   5   6
[3,]   7   8   9
[4,]   10   11   12 
> rowname
[1] "r1" "r2" "r3"
> colname=c("c1","c2","c3","c4")
> colname
[1] "c1" "c2" "c3" "c4"
> matrix(1:12,nrow=3,ncol=4,dimnames=list(rowname,colname))
  c1 c2 c3 c4
r1 1 4 7 10
r2 2 5 8 11
3   矩阵转置
A为m×n矩阵,求A'在R中可用函数t(),例如:
> A=matrix(1:12,nrow=3,ncol=4)
> A
   [,1] [,2] [,3] [,4]
[1,]   1   4   7   10
[2,]   2   5   8   11
[3,]   3   6   9   12
> t(A)
   [,1] [,2] [,3]
[1,]   1   2   3
[2,]   4   5   6
[3,]   7   8   9
[4,]   10   11   12
若将函数t()作用于一个向量x,则R默认x为列向量,返回结果为一个行向量,例如:
> x
[1] 1 2 3 4 5 6 7 8 9 10
> t(x)
  [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]   1   2   3   4   5   6   7   8   9   10
> class(x)
[1] "integer"
> class(t(x))
[1] "matrix"
若想得到一个列向量,可用t(t(x)),例如:
> x
[1] 1 2 3 4 5 6 7 8 9 10
> t(t(x))
    [,1]
[1,]   1
[2,]   2
[3,]   3
[4,]   4
[5,]   5
[6,]   6
[7,]   7
[8,]   8
[9,]   9
[10,]  10
> y=t(t(x))
> t(t(y))
    [,1]
[1,]   1
[2,]   2
[3,]   3
[4,]   4
[5,]   5
[6,]   6
[7,]   7
[8,]   8
[9,]   9
[10,]   10
4   矩阵相加减
在R中对同行同列矩阵相加减,可用符号:“+”、“-”,例如:
> 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
5   数与矩阵相乘
A为m×n矩阵,c>0,在R中求cA可用符号:“*”,例如:
> c=2
> c*A
    [,1] [,2] [,3] [,4]
[1,]   2   8   14   20
[2,]   4   10  16   22
[3,]   6   12  18   24
6   矩阵相乘
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
若A为n×m矩阵,要得到A'B,可用函数crossprod(),该函数计算结果与t(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)
    [,1] [,2] [,3]
[1,]  30  70 110
[2,]  70 174 278
[3,] 110 278 446
矩阵Hadamard积:若A={aij}m×n, B={bij}m×n, 则矩阵的Hadamard积定义为:
A⊙B={aij bij }m×n,R中Hadamard积可以直接运用运算符“*”例如:
> A=matrix(1:16,4,4)
> 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
> B=A
> A*B
    [,1] [,2] [,3] [,4]
[1,]   1   25   81 169
[2,]   4   36 100 196
[3,]   9   49 121 225
[4,]   16   64 144 256
R中这两个运算符的区别区加以注意。
7   矩阵对角元素相关运算
例如要取一个方阵的对角元素,
> A=matrix(1:16,nrow=4,ncol=4)
> 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
> diag(A)
[1] 1 6 11 16
对一个向量应用diag()函数将产生以这个向量为对角元素的对角矩阵,例如:
> diag(diag(A))
    [,1] [,2] [,3] [,4]
[1,]   1   0   0   0
[2,]   0   6   0   0
[3,]   0   0   11   0
[4,]   0   0   0   16
对一个正整数z应用diag()函数将产生以z维单位矩阵,例如:
> diag(3)
    [,1] [,2] [,3]
[1,]   1   0   0
[2,]   0   1   0
[3,]   0   0   1
8   矩阵求逆
矩阵求逆可用函数solve(),应用solve(a, b)运算结果是解线性方程组ax = b,若b缺省,则系统默认为单位矩阵,因此可用其进行矩阵求逆,例如:
> a=matrix(rnorm(16),4,4)
> a
            [,1]     [,2]     [,3]     [,4]
[1,] 1.6986019   0.5239738 0.2332094 0.3174184
[2,] -0.2010667 1.0913013 -1.2093734   0.8096514
[3,] -0.1797628 -0.7573283 0.2864535 1.3679963
[4,] -0.2217916 -0.3754700 0.1696771 -1.2424030
> solve(a)
              [,1]     [,2]     [,3]     [,4]
[1,] 0.9096360 0.54057479 0.7234861 1.3813059
[2,] -0.6464172 -0.91849017 -1.7546836 -2.6957775
[3,] -0.7841661 -1.78780083 -1.5795262 -3.1046207
[4,] -0.0741260 -0.06308603 0.1854137 -0.6607851
> solve (a) %*%a
                [,1]       [,2]           [,3]       [,4]
[1,] 1.000000e+00 2.748453e-17 -2.787755e-17 -8.023096e-17
[2,] 1.626303e-19 1.000000e+00 -4.960225e-18 6.977925e-16
[3,] 2.135878e-17 -4.629543e-17 1.000000e+00 6.201636e-17
[4,] 1.866183e-17 1.563962e-17 1.183813e-17 1.000000e+00
9   矩阵的特征值与特征向量
矩阵A的谱分解为A=UΛU',其中Λ是由A的特征值组成的对角矩阵,U的列为A的特征值对应的特征向量,在R中可以用函数eigen()函数得到U和Λ,
> args(eigen)
function (x, symmetric, only.values = FALSE, EISPACK = FALSE)
其中:x为矩阵,symmetric项指定矩阵x是否为对称矩阵,若不指定,系统将自动检测x是否为对称矩阵。例如:
> A=diag(4)+1
> A
  [,1] [,2] [,3] [,4]
[1,]   2   1   1   1
[2,]   1   2   1   1
[3,]   1   1   2   1
[4,]   1   1   1   2
> A.eigen=eigen(A,symmetric=T)
> A.eigen
$values
[1] 5 1 1 1

$vectors
        [,1]     [,2]       [,3]     [,4]
[1,] 0.5 0.8660254 0.000000e+00 0.0000000
[2,] 0.5 -0.2886751 -6.408849e-17 0.8164966
[3,] 0.5 -0.2886751 -7.071068e-01 -0.4082483
[4,] 0.5 -0.2886751 7.071068e-01 -0.4082483

> A.eigen$vectors%*%diag(A.eigen$values)%*%t(A.eigen$vectors)
  [,1] [,2] [,3] [,4]
[1,]   2   1   1   1
[2,]   1   2   1   1
[3,]   1   1   2   1
[4,]   1   1   1   2
> t(A.eigen$vectors)%*%A.eigen$vectors
            [,1]       [,2]         [,3]         [,4]
[1,] 1.000000e+00 4.377466e-17 1.626303e-17 -5.095750e-18
[2,] 4.377466e-17 1.000000e+00 -1.694066e-18 6.349359e-18
[3,] 1.626303e-17 -1.694066e-18 1.000000e+00 -1.088268e-16
[4,] -5.095750e-18 6.349359e-18 -1.088268e-16 1.000000e+00
10   矩阵的Choleskey分解
  对于正定矩阵A,可对其进行Choleskey分解,即:A=P'P,其中P为上三角矩阵,在R中可以用函数chol()进行Choleskey分解,例如:
> A
  [,1] [,2] [,3] [,4]
[1,]   2   1   1   1
[2,]   1   2   1   1
[3,]   1   1   2   1
[4,]   1   1   1   2
> chol(A)
        [,1]     [,2]     [,3]     [,4]
[1,] 1.414214 0.7071068 0.7071068 0.7071068
[2,] 0.000000 1.2247449 0.4082483 0.4082483
[3,] 0.000000 0.0000000 1.1547005 0.2886751
[4,] 0.000000 0.0000000 0.0000000 1.1180340
> t(chol(A))%*%chol(A)
  [,1] [,2] [,3] [,4]
[1,]   2   1   1   1
[2,]   1   2   1   1
[3,]   1   1   2   1
[4,]   1   1   1   2
> crossprod(chol(A),chol(A))
  [,1] [,2] [,3] [,4]
[1,]   2   1   1   1
[2,]   1   2   1   1
[3,]   1   1   2   1
[4,]   1   1   1   2
若矩阵为对称正定矩阵,可以利用Choleskey分解求行列式的值,如:
> prod(diag(chol(A))^2)
[1] 5
> det(A)
[1] 5
若矩阵为对称正定矩阵,可以利用Choleskey分解求矩阵的逆,这时用函数chol2inv(),这种用法更有效。如:
> chol2inv(chol(A))
      [,1] [,2] [,3] [,4]
[1,] 0.8 -0.2 -0.2 -0.2
[2,] -0.2 0.8 -0.2 -0.2
[3,] -0.2 -0.2 0.8 -0.2
[4,] -0.2 -0.2 -0.2 0.8
> solve(A)
  [,1] [,2] [,3] [,4]
[1,] 0.8 -0.2 -0.2 -0.2
[2,] -0.2 0.8 -0.2 -0.2
[3,] -0.2 -0.2 0.8 -0.2
[4,] -0.2 -0.2 -0.2 0.8
11   矩阵奇异值分解
  A为m×n矩阵,rank(A)= r, 可以分解为:A=UDV',其中U'U=V'V=I。在R中可以用函数scd()进行奇异值分解,例如:
> A=matrix(1:18,3,6)
> A
  [,1] [,2] [,3] [,4] [,5] [,6]
[1,]   1   4   7   10   13   16
[2,]   2   5   8   11   14   17
[3,]   3   6   9   12   15   18
> svd(A)
$d
[1] 4.589453e+01 1.640705e+00 3.627301e-16
  $u
          [,1]     [,2]     [,3]
[1,] -0.5290354 0.74394551 0.4082483
[2,] -0.5760715 0.03840487 -0.8164966
[3,] -0.6231077 -0.66713577 0.4082483
$v
          [,1]     [,2]     [,3]
[1,] -0.07736219 -0.7196003 -0.18918124
[2,] -0.19033085 -0.5089325 0.42405898
[3,] -0.30329950 -0.2982646 -0.45330031
[4,] -0.41626816 -0.0875968 -0.01637004
[5,] -0.52923682 0.1230711 0.64231130
[6,] -0.64220548 0.3337389 -0.40751869
> A.svd=svd(A)
> A.svd$u%*%diag(A.svd$d)%*%t(A.svd$v)
  [,1] [,2] [,3] [,4] [,5] [,6]
[1,]   1   4   7   10   13   16
[2,]   2   5   8   11   14   17
[3,]   3   6   9   12   15   18
> t(A.svd$u)%*%A.svd$u
            [,1]       [,2]       [,3]
[1,] 1.000000e+00 -1.169312e-16 -3.016793e-17
[2,] -1.169312e-16 1.000000e+00 -3.678156e-17
[3,] -3.016793e-17 -3.678156e-17 1.000000e+00
> t(A.svd$v)%*%A.svd$v
        [,1]       [,2]       [,3]
[1,] 1.000000e+00 8.248068e-17 -3.903128e-18
[2,] 8.248068e-17 1.000000e+00 -2.103352e-17
[3,] -3.903128e-18 -2.103352e-17 1.000000e+00
12   矩阵QR分解
A为m×n矩阵可以进行QR分解,A=QR,其中:Q'Q=I,在R中可以用函数qr()进行QR分解,例如:
> A=matrix(1:16,4,4)
> qr(A)
$qr
      [,1]     [,2]       [,3]       [,4]
[1,] -5.4772256 -12.7801930 -2.008316e+01 -2.738613e+01
[2,] 0.3651484 -3.2659863 -6.531973e+00 -9.797959e+00
[3,] 0.5477226 -0.3781696 2.641083e-15 2.056562e-15
[4,] 0.7302967 -0.9124744 8.583032e-01 -2.111449e-16

$rank
[1] 2

$qraux
[1] 1.182574e+00 1.156135e+00 1.513143e+00 2.111449e-16

$pivot
[1] 1 2 3 4

attr(,"class")
[1] "qr"
rank项返回矩阵的秩,qr项包含了矩阵Q和R的信息,要得到矩阵Q和R,可以用函数qr.Q()和qr.R()作用qr()的返回结果,例如:
> qr.R(qr(A))
      [,1]     [,2]       [,3]       [,4]
[1,] -5.477226 -12.780193 -2.008316e+01 -2.738613e+01
[2,] 0.000000 -3.265986 -6.531973e+00 -9.797959e+00
[3,] 0.000000   0.000000 2.641083e-15 2.056562e-15
[4,] 0.000000   0.000000 0.000000e+00 -2.111449e-16
> qr.Q(qr(A))
      [,1]       [,2]     [,3]     [,4]
[1,] -0.1825742 -8.164966e-01 -0.4000874 -0.37407225
[2,] -0.3651484 -4.082483e-01 0.2546329 0.79697056
[3,] -0.5477226 -8.131516e-19 0.6909965 -0.47172438
[4,] -0.7302967 4.082483e-01 -0.5455419 0.04882607
> qr.Q(qr(A))%*%qr.R(qr(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
> t(qr.Q(qr(A)))%*%qr.Q(qr(A))
        [,1]       [,2]       [,3]       [,4]
[1,] 1.000000e+00 -1.457168e-16 -6.760001e-17 -7.659550e-17
[2,] -1.457168e-16 1.000000e+00 -4.269046e-17 7.011739e-17
[3,] -6.760001e-17 -4.269046e-17 1.000000e+00 -1.596437e-16
[4,] -7.659550e-17 7.011739e-17 -1.596437e-16 1.000000e+00
> qr.X(qr(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
13   矩阵广义逆(Moore-Penrose)
  n×m矩阵A+称为m×n矩阵A的Moore-Penrose逆,如果它满足下列条件:
①   A A+A=A;②A+A A+= A+;③(A A+)H=A A+;④(A+A)H= A+A
在R的MASS包中的函数ginv()可计算矩阵A的Moore-Penrose逆,例如:
library(“MASS”)
> 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
> ginv(A)
    [,1]   [,2] [,3]   [,4]
[1,] -0.285 -0.1075 0.07 0.2475
[2,] -0.145 -0.0525 0.04 0.1325
[3,] -0.005 0.0025 0.01 0.0175
[4,] 0.135 0.0575 -0.02 -0.0975
验证性质1:
> A%*%ginv(A)%*%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
验证性质2:
> ginv(A)%*%A%*%ginv(A)
    [,1]   [,2] [,3]   [,4]
[1,] -0.285 -0.1075 0.07 0.2475
[2,] -0.145 -0.0525 0.04 0.1325
[3,] -0.005 0.0025 0.01 0.0175
[4,] 0.135 0.0575 -0.02 -0.0975
验证性质3:
> t(A%*%ginv(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
> A%*%ginv(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
验证性质4:
> 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
> 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
14   矩阵Kronecker积
  n×m矩阵A与h×k矩阵B的kronecker积为一个nh×mk维矩阵,
在R中kronecker积可以用函数kronecker()来计算,例如:
> A=matrix(1:4,2,2)
> B=matrix(rep(1,4),2,2)
> A
  [,1] [,2]
[1,]   1   3
[2,]   2   4
> B
  [,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
15   矩阵的维数
  在R中很容易得到一个矩阵的维数,函数dim()将返回一个矩阵的维数,nrow()返回行数,ncol()返回列数,例如:
  > A=matrix(1:12,3,4)
> A
  [,1] [,2] [,3] [,4]
[1,]   1   4   7   10
[2,]   2   5   8   11
[3,]   3   6   9   12
> nrow(A)
[1] 3
> ncol(A)
[1] 4
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

原文: http://blog.csdn.net/hawksoft/article/details/7773323
0 Comments



Leave a Reply.

    Yamol

    Devoting to R/ Bioinformatics/ Data Mining/ Big Data/ Cooking.

    Archives

    August 2018
    July 2017
    June 2016
    March 2016
    March 2015
    May 2014
    March 2014
    January 2014
    December 2013
    November 2013
    August 2013
    July 2013

    Categories

    All
    Creation
    Fragmental Study
    Fragmental Study
    Order
    Package
    Systematic Learn

    RSS Feed

Powered by Create your own unique website with customizable templates.