UPDATE: 2021-09-26 13:55:35

はじめに

ここではデータ分析のための線形代数の章末問題をRで回答していきます。

表1.1が表すアンケートデータ

D <- matrix(c(1,3,1,1,4,5,26,
              1,2,1,3,1,3,22,
              2,3,1,1,1,5,27,
              1,2,2,1,1,5,34,
              1,2,2,2,3,2,28,
              3,4,1,1,3,5,20,
              4,1,2,3,4,4,18,
              4,4,1,2,4,4,30,
              2,3,2,4,2,4,23,
              1,2,2,1,4,1,32),
            nrow = 10, ncol = 7, byrow = TRUE)
colnames(D) <- paste0("Q", 1:7)
rownames(D) <- paste0("ID", 1:10)

# 1.1のアンケートデータ
D
##      Q1 Q2 Q3 Q4 Q5 Q6 Q7
## ID1   1  3  1  1  4  5 26
## ID2   1  2  1  3  1  3 22
## ID3   2  3  1  1  1  5 27
## ID4   1  2  2  1  1  5 34
## ID5   1  2  2  2  3  2 28
## ID6   3  4  1  1  3  5 20
## ID7   4  1  2  3  4  4 18
## ID8   4  4  1  2  4  4 30
## ID9   2  3  2  4  2  4 23
## ID10  1  2  2  1  4  1 32

6章の章末問題

6.1 次の1から6の行列の固有値と固有ベクトルを求めよ。重根の場合には、固有ベクトルは直交していても、直交していなくてもよい。

\[ (1) \begin{bmatrix} -1 & 2 \\ 2 & -1 \end{bmatrix}, (2) \begin{bmatrix} 3 & -1 \\ -1 & 3 \end{bmatrix}, (3) \begin{bmatrix} -1 & \sqrt{2} \\ \sqrt{2} & 0 \end{bmatrix}, \\ (4) \begin{bmatrix} 0 & 2 & -1 \\ 2 & 0 & 0 \\ -1 & 0 & 0 \end{bmatrix}, (5) \begin{bmatrix} 2 & \sqrt{2} & 0 \\ \sqrt{2} & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}, (6) \begin{bmatrix} 2 & 0 \\ 0 & 2 \end{bmatrix} \]

m1 <- matrix(
  c(-1, 2,
    2,-1),
   nrow = 2, ncol = 2, byrow = TRUE)

m2 <- matrix(
  c(3,-1,
    -1, 3), 
  nrow = 2, ncol = 2, byrow = TRUE)

m3 <- matrix(
  c(-1, sqrt(2),
   sqrt(2), 0),
   nrow = 2, ncol = 2, byrow = TRUE)

m4 <- matrix(
  c(0, 2,-1,
    2, 0, 0,
    -1, 0, 0),
  nrow = 3, ncol = 3, byrow = TRUE)

m5 <- matrix(
  c(2, sqrt(2), 0,
    sqrt(2), 1, 0,
    0, 0, 1),
  nrow = 3, ncol = 3, byrow = TRUE)

m6 <- matrix(
  c(2, 0,
    0, 2),
    nrow = 2, ncol = 2, byrow = TRUE)

6.1.1

# 固有ベクトルは1/sqrt(2)であり、基準化されている
eigen_m1 <- eigen(m1)
eigen_m1_val <- eigen_m1$values
eigen_m1_vec <- eigen_m1$vectors
eigen_m1_val;eigen_m1_vec
## [1]  1 -3
##           [,1]       [,2]
## [1,] 0.7071068  0.7071068
## [2,] 0.7071068 -0.7071068

6.1.2

# 固有ベクトルは1/sqrt(2)であり、基準化されている
eigen_m2 <- eigen(m2)
eigen_m2_val <- eigen_m2$values
eigen_m2_vec <- eigen_m2$vectors
eigen_m2_val;eigen_m2_vec
## [1] 4 2
##            [,1]       [,2]
## [1,] -0.7071068 -0.7071068
## [2,]  0.7071068 -0.7071068

6.1.3

# 固有ベクトルは基準化されている
eigen_m3 <- eigen(m3)
eigen_m3_val <- eigen_m3$values
eigen_m3_vec <- eigen_m3$vectors
eigen_m3_val;eigen_m3_vec
## [1]  1 -2
##            [,1]       [,2]
## [1,] -0.5773503 -0.8164966
## [2,] -0.8164966  0.5773503

6.1.4

# 固有ベクトルは基準化されている
eigen_m4 <- eigen(m4)
eigen_m4_val <- eigen_m4$values
eigen_m4_vec <- eigen_m4$vectors
eigen_m4_val;eigen_m4_vec
## [1]  2.236068  0.000000 -2.236068
##            [,1]      [,2]       [,3]
## [1,] -0.7071068 0.0000000  0.7071068
## [2,] -0.6324555 0.4472136 -0.6324555
## [3,]  0.3162278 0.8944272  0.3162278

6.1.5

# 固有ベクトルは基準化されている
eigen_m5 <- eigen(m5)
eigen_m5_val <- eigen_m5$values
eigen_m5_vec <- eigen_m5$vectors
eigen_m5_val;eigen_m5_vec
## [1] 3.000000e+00 1.000000e+00 2.220446e-15
##           [,1] [,2]       [,3]
## [1,] 0.8164966    0  0.5773503
## [2,] 0.5773503    0 -0.8164966
## [3,] 0.0000000    1  0.0000000

6.1.6

# 固有値は重根
# ベクトルは基準化されていない
eigen_m6 <- eigen(m6)
eigen_m6_val <- eigen_m6$values
eigen_m6_vec <- eigen_m6$vectors
eigen_m6_val;eigen_m6_vec
## [1] 2 2
##      [,1] [,2]
## [1,]    0   -1
## [2,]    1    0

6.2 問題6.1の1から6の行列の固有値について、トレースと固有値の和が等しいことを示せ。

6.2.1

all.equal(sum(diag(m1)), sum(eigen_m1_val))
## [1] TRUE

6.2.2

all.equal(sum(diag(m2)), sum(eigen_m2_val))
## [1] TRUE

6.2.3

all.equal(sum(diag(m3)), sum(eigen_m3_val))
## [1] TRUE

6.2.4

all.equal(sum(diag(m4)), sum(eigen_m4_val))
## [1] TRUE

6.2.5

all.equal(sum(diag(m5)), sum(eigen_m5_val))
## [1] TRUE

6.2.6

all.equal(sum(diag(m6)), sum(eigen_m6_val))
## [1] TRUE

6.3 問題6.1の1から6の行列の固有値について、行列式と固有値の積が等しいことを示せ。

6.3.1

all.equal(det(m1), prod(eigen_m1_val))
## [1] TRUE

6.3.2

all.equal(det(m2), prod(eigen_m2_val))
## [1] TRUE

6.3.3

all.equal(det(m3), prod(eigen_m3_val))
## [1] TRUE

6.3.4

all.equal(det(m4), prod(eigen_m4_val))
## [1] TRUE

6.3.5

all.equal(det(m5), prod(eigen_m5_val))
## [1] TRUE

6.3.6

all.equal(det(m6), prod(eigen_m6_val))
## [1] TRUE

6.4 問題6.1の1から6の行列を固有ベクトルを用いて対角化せよ。重根をもつ行列の固有ベクトルについて、重複度と同じ個数の固有ベクトルを互いに直行するように求めること。

6.4.1

t(eigen_m1_vec) %*% m1  %*% eigen_m1_vec
##      [,1] [,2]
## [1,]    1    0
## [2,]    0   -3

6.4.2

t(eigen_m2_vec) %*% m2 %*% eigen_m2_vec
##      [,1] [,2]
## [1,]    4    0
## [2,]    0    2

6.4.3

t(eigen_m3_vec) %*% m3 %*% eigen_m3_vec
##      [,1] [,2]
## [1,]    1    0
## [2,]    0   -2

6.4.4

t(eigen_m4_vec) %*% m4 %*% eigen_m4_vec
##               [,1] [,2]          [,3]
## [1,]  2.236068e+00    0 -5.828671e-16
## [2,]  0.000000e+00    0  0.000000e+00
## [3,] -4.163336e-16    0 -2.236068e+00

6.4.5

t(eigen_m5_vec) %*% m5 %*% eigen_m5_vec
##              [,1] [,2]          [,3]
## [1,] 3.000000e+00    0  8.881784e-16
## [2,] 0.000000e+00    1  0.000000e+00
## [3,] 8.643898e-16    0 -6.865411e-17

6.4.6

t(eigen_m6_vec) %*% m6 %*% eigen_m6_vec
##      [,1] [,2]
## [1,]    2    0
## [2,]    0    2

6.5 次の4次行列の固有値と固有ベクトルを求めよ。重根の場合には、等しい固有値に対応する固有ベクトルが互いに直交するように求め、すべての固有ベクトルが直行することを示せ。

\[ \begin{bmatrix} 0 & 2 & 0 & 0 \\ 2 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \end{bmatrix} \]

m <- matrix(
  c(0, 2, 0, 0,
    2, 0, 0, 0,
    0, 0, 0, 1,
    0, 0, 1, 0),
  nrow = 4, ncol = 4, byrow = TRUE)
m
##      [,1] [,2] [,3] [,4]
## [1,]    0    2    0    0
## [2,]    2    0    0    0
## [3,]    0    0    0    1
## [4,]    0    0    1    0
m_eigen_vec <- eigen(m)$vectors
conbination <- t(combn(4,2))

for (i in 1:nrow(conbination)){
  tmp <- conbination[i,]
  m1 <- m_eigen_vec[, tmp[1]]
  m2 <- m_eigen_vec[, tmp[2]]
  print(sprintf("vec%.0f・vec%.0f: %.0f", tmp[1], tmp[2], t(m1) %*% m2))
}
## [1] "vec1・vec2: 0"
## [1] "vec1・vec3: 0"
## [1] "vec1・vec4: 0"
## [1] "vec2・vec3: 0"
## [1] "vec2・vec4: 0"
## [1] "vec3・vec4: 0"

6.6 表1.1の回答者1から5の設問5と6に対する回答を用いて、1から5に答えよ。

D56 <- D[1:5, 5:6]
n <- nrow(D56)
D56_Q5_mean <- mean(D56[,1])
D56_Q6_mean <- mean(D56[,2])
D56_hensa <- cbind(D56[,1] - D56_Q5_mean, D56[,2] - D56_Q6_mean)
D56_hensa
##     [,1] [,2]
## ID1    2    1
## ID2   -1   -1
## ID3   -1    1
## ID4   -1    1
## ID5    1   -2

6.6.1 偏差行列を\(\sqrt(5)(\sqrt(n))\)で割った行列を求めよ。

c <- (1/sqrt(n)) * D56_hensa
c
##           [,1]       [,2]
## ID1  0.8944272  0.4472136
## ID2 -0.4472136 -0.4472136
## ID3 -0.4472136  0.4472136
## ID4 -0.4472136  0.4472136
## ID5  0.4472136 -0.8944272

6.6.2 「偏差行列を\(\sqrt(5)\)で割った行列の転置行列」と「偏差行列を\(\sqrt(5)\)で割った行列」の積を求めよ。

# 1/n * t(D56_hensa) %*% D56_hensa
cov <- t(c) %*% c
cov
##      [,1] [,2]
## [1,]  1.6 -0.2
## [2,] -0.2  1.6

6.6.3 上記の2で得られた行列の固有値と固有ベクトルを求めよ。固有ベクトルは基準化すること。

eigen_cov <- eigen(cov)
eigen_cov
## eigen() decomposition
## $values
## [1] 1.8 1.4
## 
## $vectors
##            [,1]       [,2]
## [1,] -0.7071068 -0.7071068
## [2,]  0.7071068 -0.7071068

6.6.4 上記の3で得られた固有ベクトルを列に持つ直交行列を求めよ。

eigen_cov_vec <- eigen_cov$vectors
eigen_cov_vec
##            [,1]       [,2]
## [1,] -0.7071068 -0.7071068
## [2,]  0.7071068 -0.7071068

6.6.5 上記4で得られた直交行列の転置行列を「偏差行列を\(\sqrt(5)\)で割った行列」の転置行列に左から乗じ、得られた行列から5人の回答者の散布図を描け。

# 対角化
# 1. eigen_cov$values
# 2. (1/sqrt(n) * t(eigen_cov_vec) %*% t(D56_hensa)) %*% (1/sqrt(n) * D56_hensa %*% eigen_cov_vec)
# 1と2は同じ

# (1/sqrt(n) * eigen_cov_vec %*% t(D56_hensa))
rotate <- c %*% eigen_cov$vectors
rotate
##           [,1]       [,2]
## ID1 -0.3162278 -0.9486833
## ID2  0.0000000  0.6324555
## ID3  0.6324555  0.0000000
## ID4  0.6324555  0.0000000
## ID5 -0.9486833  0.3162278
plot(rotate, xlim = c(-1, 1), ylim = c(-1, 1), asp = TRUE)
text(rotate, row.names(rotate), cex=0.6, pos=4, col="red")
arrows(0, 0, eigen_cov_vec[1,1], eigen_cov_vec[2,1])
arrows(0, 0, eigen_cov_vec[1,2], eigen_cov_vec[2,2])
text(t(eigen_cov_vec), c("Q5","Q6"), cex=0.6, pos=4, col="blue")
abline(h = 0, v = 0)

# セッション情報

sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] ja_JP.UTF-8/ja_JP.UTF-8/ja_JP.UTF-8/C/ja_JP.UTF-8/ja_JP.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] compiler_4.0.3    magrittr_2.0.1    ragg_1.1.3        tools_4.0.3      
##  [5] htmltools_0.5.1.1 yaml_2.2.1        stringi_1.5.3     rmarkdown_2.6    
##  [9] highr_0.8         knitr_1.33        stringr_1.4.0     xfun_0.24        
## [13] digest_0.6.27     textshaping_0.3.5 systemfonts_1.0.2 rlang_0.4.10     
## [17] evaluate_0.14