天天看點

R矩陣運算總結1 矩陣基本操作2 矩陣計算3 矩陣進階操作4 解方程組5 其它

原文出處:https://www.cnblogs.com/wentingtu/archive/2012/03/30/2425582.html

1 矩陣基本操作

1.1建立向量

R裡面有多種方法來建立向量(Vector),最簡單的是用函數c()。例如:

>X=c(1,2,3,4)

>X

[1] 1 2 3 4

當然,還有别的方法。例如:

>X=1:4

>X

[1] 1 2 3 4

還有seq()函數。例如:

> X=seq(1,4,length=4)

> X

[1] 1 2 3 4

注意一點,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)

> X

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

[1,]    1    5    9

[2,]    2    6   10

[3,]    3    7   11

[4,]    4    8   12

> dim(X)

[1] 4 3

隻傳回行維數:

> 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)

> X

     [,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”)))

> X

  boy girl

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”)

> X

  boy girl

E   1    2

F   3    4

G   5    6

> colnames(X)=c(“man”,”woman”)

> X

  man woman

E   1     2

F   3     4

G   5     6

1.5 矩陣元素的檢視及其重新指派

檢視矩陣的某個特定元素,隻需要知道該元素的行坐标和列坐标即可,例如:

> X=matrix(1:12,nrow=4,ncol=3)

> X

     [,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

> X

     [,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)

> X

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

[1,]    1    4    7

[2,]    2    5    8

[3,]    3    6    9

> diag(X)

[1] 1 5 9

當然也可以對對角元素進行重新指派:

> diag(X)=c(0,0,1)

> X

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

[1,]    0    4    7

[2,]    2    0    8

[3,]    3    6    1

當操作對象不是對稱矩陣時,diag()也可以進行操作。

> X=matrix(1:12,nrow=4,ncol=3)

> X

     [,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)

> X

     [,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  8 12

改變指派:

> X[lower.tri(X)]=0

> X

     [,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)

> X

     [,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)

> 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

2.3行列式的值

R中的函數det()将計算方陣A的行列式。例如:

> X=matrix(rnorm(9),nrow=3,ncol=3)

> X

            [,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)

> A

     [,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維矩陣,公式為:

              a­11B … a1nB

Am×n×Bh×k=     …     …

               am1B … amnB   mh×nk

在R中Kronecker積可以用函數kronecher()來計算。例如:

> A=matrix(1:4,2,2)

> A

     [,1] [,2]

[1,]    1    3

[2,]    2    4

> B=matrix(rep(1,4),2,2)

> 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

2.7矩陣的伴随矩陣

求矩陣A的伴随矩陣可以用LoopAnalyst包中的函數make.adjoint()函數。例如:

>install.packages(“LoopAnalyst”)

> 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

> 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)

> A

           [,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+00 8.433738e-17 -1.341700e-18

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

[3,] -2.203641e-17 4.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)

> A

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

[1,]    1    4    7   10

[2,]    2    5    8   11

[3,]    3    6    9   12

> solve(A)

Error in solve.default(A) : only square matrices can be inverted

> 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

> A

     [,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.8333333 0.3333333 -0.1666667

[2,]  0.3333333 0.3333333  0.3333333

[3,] -0.1666667 0.3333333  0.8333333

> t(A%*%ginv(A))

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

[1,]  0.8333333 0.3333333 -0.1666667

[2,]  0.3333333 0.3333333  0.3333333

[3,] -0.1666667 0.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.3 X’X的逆

很多時候,我們需要計算形如X’X的逆。這很容易實作,例如:

> x=matrix(rnorm(9),ncol=3,nrow=3)

> x

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

[1,] -0.1806586 -0.76340512 0.002652331

[2,] -1.8018584  0.04467943 1.416332187

[3,]  1.2785359 -1.31653513 0.180653002

> solve(crossprod(x))

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

[1,] 1.2181837 0.9664576 1.470940

[2,] 0.9664576 1.2010110 1.204599

[3,] 1.4709402 1.2045986 2.269921

R中的strucchange包中的函數solveCrossprod()也可完成:

> args(solveCrossprod)

function (X, method = c(“qr”, “chol”, “solve”))

NULL

> solveCrossprod(x,method=”qr”)

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

[1,] 1.2181837 0.9664576 1.470940

[2,] 0.9664576 1.2010110 1.204599

[3,] 1.4709402 1.2045986 2.269921

> solveCrossprod(x,method=”chol”)

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

[1,] 1.2181837 0.9664576 1.470940

[2,] 0.9664576 1.2010110 1.204599

[3,] 1.4709402 1.2045986 2.269921

> solveCrossprod(x,method=”solve”)

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

[1,] 1.2181837 0.9664576 1.470940

[2,] 0.9664576 1.2010110 1.204599

[3,] 1.4709402 1.2045986 2.269921

2.9矩陣的特征值和特征向量

可以通過對矩陣A進行譜分解來得到矩陣的特征值和特征向量。矩陣A的譜分解如下:A=UΛU’,其中U的列為A的特征值所對應的特征向量,在R中可以用eigen()函數得到U和Λ。例如:

> args(eigen)

function (x, symmetric, only.values = FALSE, EISPACK = FALSE)

NULL

其中,x參數輸入矩陣;symmetric參數判斷矩陣是否為對稱矩陣,如果參數為空,系統将自動檢測矩陣的對稱性。例如:

> A=matrix(1:9,nrow=3,ncol=3)

> A

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

[1,]    1    4    7

[2,]    2    5    8

[3,]    3    6    9

> Aeigen=eigen(A)

> Aeigen

$values

[1]  1.611684e+01 -1.116844e+00 -4.054214e-16

$vectors

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

[1,] -0.4645473 -0.8829060  0.4082483

[2,] -0.5707955 -0.2395204 -0.8164966

[3,] -0.6770438  0.4038651  0.4082483

得到矩陣A的特征值:

> Aeigen$values

[1]  1.611684e+01 -1.116844e+00 -4.054214e-16

得到矩陣A的特征向量:

> Aeigen$vectors

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

[1,] -0.4645473 -0.8829060  0.4082483

[2,] -0.5707955 -0.2395204 -0.8164966

[3,] -0.6770438  0.4038651  0.4082483

3 矩陣進階操作

3.1 Choleskey分解

對于正定矩陣A,可以對其進行Choleskey分解,A=P’P,其中P為上三角矩陣,在R中可以用函數chol()。例如:

> A=diag(3)+1

> A

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

[1,]    2    1    1

[2,]    1    2    1

[3,]    1    1    2

> chol(A)

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

[1,] 1.414214 0.7071068 0.7071068

[2,] 0.000000 1.2247449 0.4082483

[3,] 0.000000 0.0000000 1.1547005

驗證A=P’P:

> t(chol(A))%*%chol(A)

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

[1,]    2    1    1

[2,]    1    2    1

[3,]    1    1    2

也可以用crossprod()函數:

> crossprod(chol(A),chol(A))

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

[1,]    2    1    1

[2,]    1    2    1

[3,]    1    1    2

可以用Choleskey分解來計算矩陣的行列式:

> prod(diag(chol(A))^2)

[1] 4

> det(A)

[1] 4

也可以用Choleskey分解來計算矩陣的逆,這時候可以用到函數chol2inv():

> chol2inv(chol(A))

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

[1,]  0.75 -0.25 -0.25

[2,] -0.25  0.75 -0.25

[3,] -0.25 -0.25  0.75

> solve(A)

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

[1,]  0.75 -0.25 -0.25

[2,] -0.25  0.75 -0.25

[3,] -0.25 -0.25  0.75

3.2奇異值分解

A為m×n矩陣,矩陣的秩為r。A可以分解為A=UDV’,其中U’U=V’V=I。在R中可以用函數svd()。例如:

> 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 2.294505e-15

$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.67039144

[2,] -0.19033085 -0.5089325  0.55766549

[3,] -0.30329950 -0.2982646  0.28189237

[4,] -0.41626816 -0.0875968  0.07320847

[5,] -0.52923682  0.1230711  0.12920119

[6,] -0.64220548  0.3337389 -0.37157608

> A.u%*%diag(A.d)%*%t(A.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

3.3 QR分解

A為m×n矩陣可以進行QR分解:A=QR,其中Q’Q=I,在R中可以用函數qr()來完成這個過程,例如:

> A=matrix(1:12,4,3)

> qr(A)

$qr

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

[1,] -5.4772256 -12.7801930 -2.008316e+01

[2,]  0.3651484  -3.2659863 -6.531973e+00

[3,]  0.5477226  -0.3781696  7.880925e-16

[4,]  0.7302967  -0.9124744  9.277920e-01

$rank

[1] 2

$qraux

[1] 1.182574 1.156135 1.373098

$pivot

[1] 1 2 3

attr(,”class”)

[1] “qr”

Rank傳回的是矩陣的秩。Qr項包含了Q矩陣和R矩陣的資訊。要想得到Q矩陣和R矩陣,可以用qr.Q()函數和qr.R()函數:

> qr.Q(qr(A))

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

[1,] -0.1825742 -8.164966e-01 -0.4000874

[2,] -0.3651484 -4.082483e-01  0.2546329

[3,] -0.5477226  4.938541e-17  0.6909965

[4,] -0.7302967  4.082483e-01 -0.5455419

> qr.R(qr(A))

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

[1,] -5.477226 -12.780193 -2.008316e+01

[2,]  0.000000  -3.265986 -6.531973e+00

[3,]  0.000000   0.000000  7.880925e-16

4 解方程組

4.1普通方程組

解普通方程組可以用函數solve(),solve()的基本用法是solve(A,b),其中,A為方程組的系數矩陣,b為方程組的右端。例如:

已知方程組:

2x1+2x3=1

2x1+x2+2x­3=2

2x1+x­2=3

解法如下:

> A

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

[1,]    2    0    2

[2,]    2    1    2

[3,]    2    1    0

> b=1:3

>b

[1] 1 2 3

> solve(A,b)

[1]  1.0  1.0 -0.5

即x­1=1,x2=1,x3=-0.5。

4.2 特殊方程組

對于系數矩陣是上三角矩陣和下三角矩陣的方程組。R中提供了backsolve()和fowardsolve()來解決這個問題。

backsolve(r, x, k=ncol(r), upper.tri=TRUE, transpose=FALSE)

forwardsolve(l, x, k=ncol(l), upper.tri=FALSE, transpose=FALSE)

這兩個函數都是符合操作的函數,大緻可以分為三個步驟:

①通過将系數矩陣的上三角或者下三角變為0的到新的系數矩陣,這通過upper.tri參數來實作,若upper.tri=TRUR,上三角不為0。

②通過将對步驟1中得到的新系數矩陣進行轉置得到新的系數矩陣,這通過transpose參數實作,若transpose=FALSE,則步驟1中得到的系數矩陣将被轉置。

③根據步驟2得到的系數矩陣來解方程組。

X1+4X2+7X3=1

2X1+5X2+8X3=2

3X1+6X2+9X3=3

方程組的系數矩陣為:

> A

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

[1,]    1    4    7

[2,]    2    5    8

[3,]    3    6    9

> b

[1] 1 2 3

> backsolve(A,b,upper.tri=T,transpose=F)

[1] -0.8000000 -0.1333333  0.3333333

過程分解:

①upper.tri=T,說明系數矩陣的上三角不為0。

> B=A

> B[lower.tri(B)]=0

> B

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

[1,]    1    4    7

[2,]    0    5    8

[3,]    0    0    9

②transpose=F說明系數矩陣未被轉置。

③解方程:

> solve(B,b)

[1] -0.8000000 -0.1333333  0.3333333

5 其它

5.1矩陣的向量化

将矩陣向量化有時候是必要的。矩陣的向量化可以通過as.vector()來實作:

> A

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

[1,]    1    4    7   10

[2,]    2    5    8   11

[3,]    3    6    9   12

将矩陣元素向量化:

> as.vector(A)

 [1]  1  2  3  4  5  6  7  8  9 10 11 12

将矩陣的方陣部分元素向量化:

> as.vector(A[1:min(dim(A)),1:min(dim(A))])

[1] 1 2 3 4 5 6 7 8 9

5.2矩陣的合并

5.2.1矩陣的列合并

矩陣的列合并可以通過cbind()來實作。

> A

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

[1,]    1    4    7

[2,]    2    5    8

[3,]    3    6    9

> B=1:3

> cbind(A,B)

           B

[1,] 1 4 7 1

[2,] 2 5 8 2

[3,] 3 6 9 3

5.2.2矩陣的行合并

矩陣的行合并可以通過rbind()來實作。

> A

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

[1,]    1    4    7

[2,]    2    5    8

[3,]    3    6    9

> B=1:3

> rbind(A,B)

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

     1    4    7

     2    5    8

     3    6    9

B    1    2    3

5.3 時序矩陣的滞後

在時間序列中經常會用到一個序列的滞後序列,R中的包fMultivar中的函數tslag()提供了這個功能。

> library(fMultivar)

Loading required package: sn

Loading required package: mnormt

Package ‘sn’, 0.4-16 (2010-08-30). Type ‘help(SN)’ for summary information

Loading required package: timeDate

Loading required package: timeSeries

Loading required package: fBasics

Loading required package: MASS

Attaching package: ‘fBasics’

The following object(s) are masked from ‘package:base’:

    norm

> args(tslag)

function (x, k = 1, trim = FALSE)

NULL

其中:x為一個向量,k指定滞後階數,可以是一個自然數列,若trim為假,則傳回序列與原序列長度相同,但含有NA值;若trim項為真,則傳回序列中不含有NA值,例如:

> x=1:9

> x

[1] 1 2 3 4 5 6 7 8 9

> 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

> 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