1. 程式人生 > >R軟體矩陣使用總結

R軟體矩陣使用總結

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={aijbij }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 256R中這兩個運算子的區別區加以注意。

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.23320940.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-186.977925e-16
[3,] 2.135878e-17 -4.629543e-17 1.000000e+006.201636e-17
[4,] 1.866183e-17 1.563962e-17 1.183813e-171.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-186.349359e-18
[3,] 1.626303e-17 -1.694066e-18 1.000000e+00-1.088268e-16
[4,] -5.095750e-18 6.349359e-18 -1.088268e-161.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-152.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+002.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-152.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-177.011739e-17
[3,] -6.760001e-17 -4.269046e-17 1.000000e+00-1.596437e-16
[4,] -7.659550e-17 7.011739e-17 -1.596437e-161.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