16 矩陣的行和、列和、行平均與列平均

在R中很容易求得一個矩陣的各行的和、平均數(shù)與列的和、平均數(shù)，例如：

> 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

17 矩陣X'X的逆

在統(tǒng)計計算中，我們常常需要計算這樣矩陣的逆，如OLS估計中求系數(shù)矩陣。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中，我們可以很方便的取到一個矩陣的上、下三角部分的元素，函數(shù)lower.tri()

和函數(shù)upper.tri()提供了有效的方法。

> args(lower.tri)

function (x, diag = FALSE)

函數(shù)將返回一個邏輯值矩陣，其中下三角部分為真，上三角部分為假，選項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函數(shù)

這兩個函數(shù)用于解特殊線性方程組，其特殊之處在于系數(shù)矩陣為上或下三角。

> args(backsolve)

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

> args(forwardsolve)

function (l, x, = ncol(l), upper.tri = FALSE, transpose = FALSE)upper.tri =F transpose= T

11

upper.tri =F transpose=F

例如：

> 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

upper.tri =F transpose=F

例如：

> 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