R로 단형질 임의회귀모형의 해 구하기

 

임의회귀모형을 본격적으로 프로그램하기 전에 R의 행렬 기능을 이용하여

 

임의회귀모형을 연습함

 

잔차 분산

R.txt

3.710

유전분산

G.txt

3.297 0.594 -1.381
0.594 0.921 -0.289
-1.381 -0.289 1.005

영구환경효과 분산

P.txt

6.872 -0.254 -1.101
-0.254 3.171 0.167
-1.101 0.167 2.457

자료

strrm_data.txt

1 4 4 17.0
2 38 4 18.6
3 72 4 24.0
4 106 4 20.0
5 140 4 20.0
6 174 4 15.6
7 208 4 16.0
8 242 4 13.0
9 276 4 8.2
10 310 4 8.0
1 4 5 23.0
2 38 5 21.0
3 72 5 18.0
4 106 5 17.0
5 140 5 16.2
6 174 5 14.0
7 208 5 14.2
8 242 5 13.4
9 276 5 11.8
10 310 5 11.4
6 4 6 10.4
7 38 6 12.3
8 72 6 13.2
9 106 6 11.6
10 140 6 8.4
4 4 7 22.8
5 38 7 22.4
6 72 7 21.4
7 106 7 18.8
8 140 7 18.3
9 174 7 16.2
10 208 7 15.0
1 4 8 22.2
2 38 8 20.0
3 72 8 21.0
4 106 8 23.0
5 140 8 16.8
6 174 8 11.0
7 208 8 13.0
8 242 8 17.0
9 276 8 13.0
10 310 8 12.6

혈통

strrm_pedi.txt1 0 0
2 0 0
3 0 0
4 1 2
5 3 2
6 1 5
7 3 4
8 1 7
R 소스# 파일 -> 디렉토리 변경에서 본 스크립트가 있는 디렉토리 설정setwd("D:/Users/bhpark/2011/job/프로그램개발/단형질임의회귀모형_R")library(MASS)# R - residual variance
R <- as.numeric(read.table("R.txt"))# R 확인
R# G - covariance matrix for the random regression coefficient for animal effect
G <- as.matrix(read.table("G.txt"))Ginv <- solve(G)# G 확인
G
Ginv
G %*% Ginv
# P - covariance matrix for the random regression coefficient for permanent environmental effect
P <- as.matrix(read.table("P.txt"))Pinv = solve(P)# P 확인
P
Pinv
P %*% Pinv###########
# 자료 입력 #
###########
data <- (read.table("strrm_data.txt"))
# y 추출
y <- data[,4]

# y 확인
y
# X1의 추출
X1_data <- data[,1]# X1 확인
X1_data# X1을 이용하여 incidence matrix 만들기
X1 <- contrasts( as.factor(X1_data), contrasts = FALSE)[ as.factor(X1_data), ]
X1# X1' X1 확인
t(X1) %*% X1
# X2 raw data 추출
X2_data = data[,2]# X2 확인
X2_data# X2를 이용하여 M 행렬(polynomials of the standardized DIM values) 생성
X2_1 = c(rep(1,42))
X2_2 = -1 + 2 * (X2_data - 4) / (310 - 4)
X2_3 = X2_2 * X2_2
X2_4 = X2_3 * X2_2
X2_5 = X2_4 * X2_2
M = cbind(X2_1, X2_2, X2_3, X2_4, X2_5)# Legendre polynomials 계수를 포함하는 차수 k = 5의 행렬
lambda = matrix(c(
sqrt(.5) , 0.0 , 0.0 , 0.0 , 0.0,
0.0 , sqrt(1.5) , 0.0 , 0.0 , 0.0,
-sqrt(2.5)*.5 , 0.0 , sqrt(2.5)*1.5 , 0.0 , 0.0,
0.0 , -sqrt(3.5)*1.5, 0.0 , sqrt(3.5)*2.5, 0.0,
sqrt(4.5)*3./8., 0.0 , -sqrt(4.5)*30./8., 0.0 , sqrt(4.5)*35./8.),nrow = 5)# X2
X2 = M %*% lambda# X2' X2 확인
t(X2) %*% X2# X 생성
X = cbind(X1, X2)# X 확인
X
# animal 추출
animal <- data[,3]# animal 확인
animal# Q를 만들기 위한 준비, x2의 3열을 이용
Q_data <- X2[,1:3]
Q_data# 모든 원소가 0인 42 by 24 행렬 만들기
Q <- matrix(rep(0, 42*8*3), nrow = 42)# 일종의 design matrix 만들기
for( i in 1:42 ) {
j <- (animal[i]-1)*3+1
k <- (animal[i]-1)*3+3
Q[i, j:k] <- Q_data[i,]
}# Q 확인
Q# Q' R-1 Q 확인
t(Q) %*% Q / R
# Z 생성
Z <- Q[,10:24]# Z 확인
Z# 혈통 입력
pedi = as.matrix(read.table("strrm_pedi.txt"))# 혈통 확인
pedi# inverse of numerator relationship matrix
Ainv = matrix(rep(0, 8 * 8), nrow = 8)Ainvfor ( i in 1:8) {
animal = pedi[i,1]
sire = pedi[i,2]
dam = pedi[i,3] if (sire == 0 & dam == 0) {
alpha = 1
Ainv[animal,animal] = alpha + Ainv[animal,animal]
}
else if (sire != 0 & dam == 0) {
alpha = 4/3
Ainv[animal,animal] = alpha + Ainv[animal,animal]
Ainv[animal,sire] = -alpha/2 + Ainv[animal,sire]
Ainv[sire,animal] = -alpha/2 + Ainv[sire,animal]
Ainv[sire,sire] = alpha/4 + Ainv[sire,sire]
}
else if (sire == 0 & dam != 0) {
alpha = 4/3
Ainv[animal,animal] = alpha + Ainv[animal,animal]
Ainv[animal,dam] = -alpha/2 + Ainv[animal,dam]
Ainv[dam,animal] = -alpha/2 + Ainv[dam,animal]
Ainv[dam,dam] = alpha/4 + Ainv[dam,dam]
}
else {
alpha = 2
Ainv[animal,animal] = alpha + Ainv[animal,animal]
Ainv[animal,sire] = -alpha/2 + Ainv[animal,sire]
Ainv[sire,animal] = -alpha/2 + Ainv[sire,animal]
Ainv[animal,dam] = -alpha/2 + Ainv[animal,dam]
Ainv[dam,animal] = -alpha/2 + Ainv[dam,animal]
Ainv[sire,sire] = alpha/4 + Ainv[sire,sire]
Ainv[sire,dam] = alpha/4 + Ainv[sire,dam]
Ainv[dam,sire] = alpha/4 + Ainv[dam,sire]
Ainv[dam,dam] = alpha/4 + Ainv[dam,dam]
}
}# inverse of numerator relationship matrix 확인
Ainv# LHS 생성
LHS <- rbind(
cbind(t(X) %*% X / R, t(X) %*% Q / R, t(X) %*% Z / R)
,
cbind(t(Q) %*% X / R, t(Q) %*% Q / R + Ainv %x% Ginv, t(Q) %*% Z / R)
,
cbind(t(Z) %*% X / R, t(Z) %*% Q / R, t(Z) %*% Z / R + diag(5) %x% Pinv)
)# RHS 생성
RHS <- rbind(
t(X) %*% y / R
,
t(Q) %*% y / R
,
t(Z) %*% y / R
)# 해 구하기
strrm_sol = ginv(LHS) %*% RHS# 해 확인
strrm_sol
R 실행결과> # 파일 -> 디렉토리 변경에서 본 스크립트가 있는 디렉토리 설정
>
> setwd("D:/Users/bhpark/2011/job/프로그램개발/단형질임의회귀모형_R")
>
> library(MASS)
>
> # R - residual variance
> R <- as.numeric(read.table("R.txt"))
>
> # R 확인
> R
[1] 3.71
>
> # G - covariance matrix for the random regression coefficient for animal effect
> G <- as.matrix(read.table("G.txt"))
>
> Ginv <- solve(G)
>
> # G 확인
> G
V1 V2 V3
[1,] 3.297 0.594 -1.381
[2,] 0.594 0.921 -0.289
[3,] -1.381 -0.289 1.005
> Ginv
[,1] [,2] [,3]
V1 0.7390611 -0.1736541 0.9656292
V2 -0.1736541 1.2342704 0.1163063
V3 0.9656292 0.1163063 2.3553696
> G %*% Ginv
[,1] [,2] [,3]
[1,] 1.000000e+00 -4.418124e-18 -3.382711e-17
[2,] -6.445582e-17 1.000000e+00 -8.939247e-17
[3,] 2.623227e-16 2.604796e-17 1.000000e+00
>
>
> # P - covariance matrix for the random regression coefficient for permanent environmental effect
> P <- as.matrix(read.table("P.txt"))
>
> Pinv = solve(P)
>
> # P 확인
> P
V1 V2 V3
[1,] 6.872 -0.254 -1.101
[2,] -0.254 3.171 0.167
[3,] -1.101 0.167 2.457
> Pinv
[,1] [,2] [,3]
V1 0.157023483 0.008903925 0.06975820
V2 0.008903925 0.316995726 -0.01755599
V3 0.069758201 -0.017555989 0.43945284
> P %*% Pinv
[,1] [,2] [,3]
[1,] 1.000000e+00 5.558230e-18 -3.897707e-17
[2,] -2.151464e-18 1.000000e+00 1.846532e-17
[3,] 2.340521e-17 9.893345e-19 1.000000e+00
>
> ###########
> # 자료 입력 #
> ###########
> data <- (read.table("strrm_data.txt"))
>
>
> # y 추출
> y <- data[,4]
>
> # y 확인
> y
[1] 17.0 18.6 24.0 20.0 20.0 15.6 16.0 13.0 8.2 8.0 23.0 21.0 18.0 17.0
[15] 16.2 14.0 14.2 13.4 11.8 11.4 10.4 12.3 13.2 11.6 8.4 22.8 22.4 21.4
[29] 18.8 18.3 16.2 15.0 22.2 20.0 21.0 23.0 16.8 11.0 13.0 17.0 13.0 12.6
>
>
> # X1의 추출
> X1_data <- data[,1]
>
> # X1 확인
> X1_data
[1] 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 6 7 8 9
[25] 10 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10
>
> # X1을 이용하여 incidence matrix 만들기
> X1 <- contrasts( as.factor(X1_data), contrasts = FALSE)[ as.factor(X1_data), ]
> X1
1 2 3 4 5 6 7 8 9 10
1 1 0 0 0 0 0 0 0 0 0
2 0 1 0 0 0 0 0 0 0 0
3 0 0 1 0 0 0 0 0 0 0
4 0 0 0 1 0 0 0 0 0 0
5 0 0 0 0 1 0 0 0 0 0
6 0 0 0 0 0 1 0 0 0 0
7 0 0 0 0 0 0 1 0 0 0
8 0 0 0 0 0 0 0 1 0 0
9 0 0 0 0 0 0 0 0 1 0
10 0 0 0 0 0 0 0 0 0 1
1 1 0 0 0 0 0 0 0 0 0
2 0 1 0 0 0 0 0 0 0 0
3 0 0 1 0 0 0 0 0 0 0
4 0 0 0 1 0 0 0 0 0 0
5 0 0 0 0 1 0 0 0 0 0
6 0 0 0 0 0 1 0 0 0 0
7 0 0 0 0 0 0 1 0 0 0
8 0 0 0 0 0 0 0 1 0 0
9 0 0 0 0 0 0 0 0 1 0
10 0 0 0 0 0 0 0 0 0 1
6 0 0 0 0 0 1 0 0 0 0
7 0 0 0 0 0 0 1 0 0 0
8 0 0 0 0 0 0 0 1 0 0
9 0 0 0 0 0 0 0 0 1 0
10 0 0 0 0 0 0 0 0 0 1
4 0 0 0 1 0 0 0 0 0 0
5 0 0 0 0 1 0 0 0 0 0
6 0 0 0 0 0 1 0 0 0 0
7 0 0 0 0 0 0 1 0 0 0
8 0 0 0 0 0 0 0 1 0 0
9 0 0 0 0 0 0 0 0 1 0
10 0 0 0 0 0 0 0 0 0 1
1 1 0 0 0 0 0 0 0 0 0
2 0 1 0 0 0 0 0 0 0 0
3 0 0 1 0 0 0 0 0 0 0
4 0 0 0 1 0 0 0 0 0 0
5 0 0 0 0 1 0 0 0 0 0
6 0 0 0 0 0 1 0 0 0 0
7 0 0 0 0 0 0 1 0 0 0
8 0 0 0 0 0 0 0 1 0 0
9 0 0 0 0 0 0 0 0 1 0
10 0 0 0 0 0 0 0 0 0 1
>
> # X1' X1 확인
> t(X1) %*% X1
1 2 3 4 5 6 7 8 9 10
1 3 0 0 0 0 0 0 0 0 0
2 0 3 0 0 0 0 0 0 0 0
3 0 0 3 0 0 0 0 0 0 0
4 0 0 0 4 0 0 0 0 0 0
5 0 0 0 0 4 0 0 0 0 0
6 0 0 0 0 0 5 0 0 0 0
7 0 0 0 0 0 0 5 0 0 0
8 0 0 0 0 0 0 0 5 0 0
9 0 0 0 0 0 0 0 0 5 0
10 0 0 0 0 0 0 0 0 0 5
>
>
> # X2 raw data 추출
> X2_data = data[,2]
>
> # X2 확인
> X2_data
[1] 4 38 72 106 140 174 208 242 276 310 4 38 72 106 140 174 208 242
[19] 276 310 4 38 72 106 140 4 38 72 106 140 174 208 4 38 72 106
[37] 140 174 208 242 276 310
>
> # X2를 이용하여 M 행렬(polynomials of the standardized DIM values) 생성
> X2_1 = c(rep(1,42))
> X2_2 = -1 + 2 * (X2_data - 4) / (310 - 4)
> X2_3 = X2_2 * X2_2
> X2_4 = X2_3 * X2_2
> X2_5 = X2_4 * X2_2
> M = cbind(X2_1, X2_2, X2_3, X2_4, X2_5)
>
> # Legendre polynomials 계수를 포함하는 차수 k = 5의 행렬
> lambda = matrix(c(
+ sqrt(.5) , 0.0 , 0.0 , 0.0 , 0.0,
+ 0.0 , sqrt(1.5) , 0.0 , 0.0 , 0.0,
+ -sqrt(2.5)*.5 , 0.0 , sqrt(2.5)*1.5 , 0.0 , 0.0,
+ 0.0 , -sqrt(3.5)*1.5, 0.0 , sqrt(3.5)*2.5, 0.0,
+ sqrt(4.5)*3./8., 0.0 , -sqrt(4.5)*30./8., 0.0 , sqrt(4.5)*35./8.),nrow = 5)
>
> # X2
> X2 = M %*% lambda
>
> # X2' X2 확인
> t(X2) %*% X2
[,1] [,2] [,3] [,4] [,5]
[1,] 21.0000000 -4.426352 4.058049 -0.8456242 8.715592
[2,] -4.4263521 24.629630 -4.701770 11.1718181 -3.064330
[3,] 4.0580493 -4.701770 31.063100 -6.6620895 19.087331
[4,] -0.8456242 11.171818 -6.662090 38.6554199 -8.856050
[5,] 8.7155921 -3.064330 19.087331 -8.8560497 48.292619
>
> # X 생성
> X = cbind(X1, X2)
>
> # X 확인
> X
1 2 3 4 5 6 7 8 9 10
1 1 0 0 0 0 0 0 0 0 0 0.7071068 -1.2247449 1.5811388 -1.87082869
2 0 1 0 0 0 0 0 0 0 0 0.7071068 -0.9525793 0.6441677 -0.01796406
3 0 0 1 0 0 0 0 0 0 0 0.7071068 -0.6804138 -0.0585607 0.75705688
4 0 0 0 1 0 0 0 0 0 0 0.7071068 -0.4082483 -0.5270463 0.76218947
5 0 0 0 0 1 0 0 0 0 0 0.7071068 -0.1360828 -0.7612891 0.30538905
6 0 0 0 0 0 1 0 0 0 0 0.7071068 0.1360828 -0.7612891 -0.30538905
7 0 0 0 0 0 0 1 0 0 0 0.7071068 0.4082483 -0.5270463 -0.76218947
8 0 0 0 0 0 0 0 1 0 0 0.7071068 0.6804138 -0.0585607 -0.75705688
9 0 0 0 0 0 0 0 0 1 0 0.7071068 0.9525793 0.6441677 0.01796406
10 0 0 0 0 0 0 0 0 0 1 0.7071068 1.2247449 1.5811388 1.87082869
1 1 0 0 0 0 0 0 0 0 0 0.7071068 -1.2247449 1.5811388 -1.87082869
2 0 1 0 0 0 0 0 0 0 0 0.7071068 -0.9525793 0.6441677 -0.01796406
3 0 0 1 0 0 0 0 0 0 0 0.7071068 -0.6804138 -0.0585607 0.75705688
4 0 0 0 1 0 0 0 0 0 0 0.7071068 -0.4082483 -0.5270463 0.76218947
5 0 0 0 0 1 0 0 0 0 0 0.7071068 -0.1360828 -0.7612891 0.30538905
6 0 0 0 0 0 1 0 0 0 0 0.7071068 0.1360828 -0.7612891 -0.30538905
7 0 0 0 0 0 0 1 0 0 0 0.7071068 0.4082483 -0.5270463 -0.76218947
8 0 0 0 0 0 0 0 1 0 0 0.7071068 0.6804138 -0.0585607 -0.75705688
9 0 0 0 0 0 0 0 0 1 0 0.7071068 0.9525793 0.6441677 0.01796406
10 0 0 0 0 0 0 0 0 0 1 0.7071068 1.2247449 1.5811388 1.87082869
6 0 0 0 0 0 1 0 0 0 0 0.7071068 -1.2247449 1.5811388 -1.87082869
7 0 0 0 0 0 0 1 0 0 0 0.7071068 -0.9525793 0.6441677 -0.01796406
8 0 0 0 0 0 0 0 1 0 0 0.7071068 -0.6804138 -0.0585607 0.75705688
9 0 0 0 0 0 0 0 0 1 0 0.7071068 -0.4082483 -0.5270463 0.76218947
10 0 0 0 0 0 0 0 0 0 1 0.7071068 -0.1360828 -0.7612891 0.30538905
4 0 0 0 1 0 0 0 0 0 0 0.7071068 -1.2247449 1.5811388 -1.87082869
5 0 0 0 0 1 0 0 0 0 0 0.7071068 -0.9525793 0.6441677 -0.01796406
6 0 0 0 0 0 1 0 0 0 0 0.7071068 -0.6804138 -0.0585607 0.75705688
7 0 0 0 0 0 0 1 0 0 0 0.7071068 -0.4082483 -0.5270463 0.76218947
8 0 0 0 0 0 0 0 1 0 0 0.7071068 -0.1360828 -0.7612891 0.30538905
9 0 0 0 0 0 0 0 0 1 0 0.7071068 0.1360828 -0.7612891 -0.30538905
10 0 0 0 0 0 0 0 0 0 1 0.7071068 0.4082483 -0.5270463 -0.76218947
1 1 0 0 0 0 0 0 0 0 0 0.7071068 -1.2247449 1.5811388 -1.87082869
2 0 1 0 0 0 0 0 0 0 0 0.7071068 -0.9525793 0.6441677 -0.01796406
3 0 0 1 0 0 0 0 0 0 0 0.7071068 -0.6804138 -0.0585607 0.75705688
4 0 0 0 1 0 0 0 0 0 0 0.7071068 -0.4082483 -0.5270463 0.76218947
5 0 0 0 0 1 0 0 0 0 0 0.7071068 -0.1360828 -0.7612891 0.30538905
6 0 0 0 0 0 1 0 0 0 0 0.7071068 0.1360828 -0.7612891 -0.30538905
7 0 0 0 0 0 0 1 0 0 0 0.7071068 0.4082483 -0.5270463 -0.76218947
8 0 0 0 0 0 0 0 1 0 0 0.7071068 0.6804138 -0.0585607 -0.75705688
9 0 0 0 0 0 0 0 0 1 0 0.7071068 0.9525793 0.6441677 0.01796406
10 0 0 0 0 0 0 0 0 0 1 0.7071068 1.2247449 1.5811388 1.87082869

1 2.12132034
2 -0.62045629
3 -0.77565120
4 0.02618914
5 0.69870039
6 0.69870039
7 0.02618914
8 -0.77565120
9 -0.62045629
10 2.12132034
1 2.12132034
2 -0.62045629
3 -0.77565120
4 0.02618914
5 0.69870039
6 0.69870039
7 0.02618914
8 -0.77565120
9 -0.62045629
10 2.12132034
6 2.12132034
7 -0.62045629
8 -0.77565120
9 0.02618914
10 0.69870039
4 2.12132034
5 -0.62045629
6 -0.77565120
7 0.02618914
8 0.69870039
9 0.69870039
10 0.02618914
1 2.12132034
2 -0.62045629
3 -0.77565120
4 0.02618914
5 0.69870039
6 0.69870039
7 0.02618914
8 -0.77565120
9 -0.62045629
10 2.12132034
>
>
> # animal 추출
> animal <- data[,3]
>
> # animal 확인
> animal
[1] 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 6 6 6 6 6 7 7 7 7 7 7 7 8 8 8 8
[37] 8 8 8 8 8 8
>
> # Q를 만들기 위한 준비, x2의 3열을 이용
> Q_data <- X2[,1:3]
> Q_data
[,1] [,2] [,3]
[1,] 0.7071068 -1.2247449 1.5811388
[2,] 0.7071068 -0.9525793 0.6441677
[3,] 0.7071068 -0.6804138 -0.0585607
[4,] 0.7071068 -0.4082483 -0.5270463
[5,] 0.7071068 -0.1360828 -0.7612891
[6,] 0.7071068 0.1360828 -0.7612891
[7,] 0.7071068 0.4082483 -0.5270463
[8,] 0.7071068 0.6804138 -0.0585607
[9,] 0.7071068 0.9525793 0.6441677
[10,] 0.7071068 1.2247449 1.5811388
[11,] 0.7071068 -1.2247449 1.5811388
[12,] 0.7071068 -0.9525793 0.6441677
[13,] 0.7071068 -0.6804138 -0.0585607
[14,] 0.7071068 -0.4082483 -0.5270463
[15,] 0.7071068 -0.1360828 -0.7612891
[16,] 0.7071068 0.1360828 -0.7612891
[17,] 0.7071068 0.4082483 -0.5270463
[18,] 0.7071068 0.6804138 -0.0585607
[19,] 0.7071068 0.9525793 0.6441677
[20,] 0.7071068 1.2247449 1.5811388
[21,] 0.7071068 -1.2247449 1.5811388
[22,] 0.7071068 -0.9525793 0.6441677
[23,] 0.7071068 -0.6804138 -0.0585607
[24,] 0.7071068 -0.4082483 -0.5270463
[25,] 0.7071068 -0.1360828 -0.7612891
[26,] 0.7071068 -1.2247449 1.5811388
[27,] 0.7071068 -0.9525793 0.6441677
[28,] 0.7071068 -0.6804138 -0.0585607
[29,] 0.7071068 -0.4082483 -0.5270463
[30,] 0.7071068 -0.1360828 -0.7612891
[31,] 0.7071068 0.1360828 -0.7612891
[32,] 0.7071068 0.4082483 -0.5270463
[33,] 0.7071068 -1.2247449 1.5811388
[34,] 0.7071068 -0.9525793 0.6441677
[35,] 0.7071068 -0.6804138 -0.0585607
[36,] 0.7071068 -0.4082483 -0.5270463
[37,] 0.7071068 -0.1360828 -0.7612891
[38,] 0.7071068 0.1360828 -0.7612891
[39,] 0.7071068 0.4082483 -0.5270463
[40,] 0.7071068 0.6804138 -0.0585607
[41,] 0.7071068 0.9525793 0.6441677
[42,] 0.7071068 1.2247449 1.5811388
>
> # 모든 원소가 0인 42 by 24 행렬 만들기
> Q <- matrix(rep(0, 42*8*3), nrow = 42)
>
> # 일종의 design matrix 만들기
> for( i in 1:42 ) {
+ j <- (animal[i]-1)*3+1
+ k <- (animal[i]-1)*3+3
+ Q[i, j:k] <- Q_data[i,]
+ }
>
> # Q 확인
> Q
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
[1,] 0 0 0 0 0 0 0 0 0 0.7071068 -1.2247449
[2,] 0 0 0 0 0 0 0 0 0 0.7071068 -0.9525793
[3,] 0 0 0 0 0 0 0 0 0 0.7071068 -0.6804138
[4,] 0 0 0 0 0 0 0 0 0 0.7071068 -0.4082483
[5,] 0 0 0 0 0 0 0 0 0 0.7071068 -0.1360828
[6,] 0 0 0 0 0 0 0 0 0 0.7071068 0.1360828
[7,] 0 0 0 0 0 0 0 0 0 0.7071068 0.4082483
[8,] 0 0 0 0 0 0 0 0 0 0.7071068 0.6804138
[9,] 0 0 0 0 0 0 0 0 0 0.7071068 0.9525793
[10,] 0 0 0 0 0 0 0 0 0 0.7071068 1.2247449
[11,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[12,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[13,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[14,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[15,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[16,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[17,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[18,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[19,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[20,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[21,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[22,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[23,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[24,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[25,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[26,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[27,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[28,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[29,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[30,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[31,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[32,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[33,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[34,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[35,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[36,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[37,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[38,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[39,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[40,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[41,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[42,] 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000
[,12] [,13] [,14] [,15] [,16] [,17]
[1,] 1.5811388 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[2,] 0.6441677 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[3,] -0.0585607 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[4,] -0.5270463 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[5,] -0.7612891 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[6,] -0.7612891 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[7,] -0.5270463 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[8,] -0.0585607 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[9,] 0.6441677 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[10,] 1.5811388 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[11,] 0.0000000 0.7071068 -1.2247449 1.5811388 0.0000000 0.0000000
[12,] 0.0000000 0.7071068 -0.9525793 0.6441677 0.0000000 0.0000000
[13,] 0.0000000 0.7071068 -0.6804138 -0.0585607 0.0000000 0.0000000
[14,] 0.0000000 0.7071068 -0.4082483 -0.5270463 0.0000000 0.0000000
[15,] 0.0000000 0.7071068 -0.1360828 -0.7612891 0.0000000 0.0000000
[16,] 0.0000000 0.7071068 0.1360828 -0.7612891 0.0000000 0.0000000
[17,] 0.0000000 0.7071068 0.4082483 -0.5270463 0.0000000 0.0000000
[18,] 0.0000000 0.7071068 0.6804138 -0.0585607 0.0000000 0.0000000
[19,] 0.0000000 0.7071068 0.9525793 0.6441677 0.0000000 0.0000000
[20,] 0.0000000 0.7071068 1.2247449 1.5811388 0.0000000 0.0000000
[21,] 0.0000000 0.0000000 0.0000000 0.0000000 0.7071068 -1.2247449
[22,] 0.0000000 0.0000000 0.0000000 0.0000000 0.7071068 -0.9525793
[23,] 0.0000000 0.0000000 0.0000000 0.0000000 0.7071068 -0.6804138
[24,] 0.0000000 0.0000000 0.0000000 0.0000000 0.7071068 -0.4082483
[25,] 0.0000000 0.0000000 0.0000000 0.0000000 0.7071068 -0.1360828
[26,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[27,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[28,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[29,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[30,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[31,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[32,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[33,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[34,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[35,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[36,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[37,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[38,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[39,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[40,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[41,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[42,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[,18] [,19] [,20] [,21] [,22] [,23]
[1,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[2,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[3,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[4,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[5,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[6,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[7,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[8,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[9,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[10,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[11,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[12,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[13,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[14,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[15,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[16,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[17,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[18,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[19,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[20,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[21,] 1.5811388 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[22,] 0.6441677 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[23,] -0.0585607 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[24,] -0.5270463 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[25,] -0.7612891 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[26,] 0.0000000 0.7071068 -1.2247449 1.5811388 0.0000000 0.0000000
[27,] 0.0000000 0.7071068 -0.9525793 0.6441677 0.0000000 0.0000000
[28,] 0.0000000 0.7071068 -0.6804138 -0.0585607 0.0000000 0.0000000
[29,] 0.0000000 0.7071068 -0.4082483 -0.5270463 0.0000000 0.0000000
[30,] 0.0000000 0.7071068 -0.1360828 -0.7612891 0.0000000 0.0000000
[31,] 0.0000000 0.7071068 0.1360828 -0.7612891 0.0000000 0.0000000
[32,] 0.0000000 0.7071068 0.4082483 -0.5270463 0.0000000 0.0000000
[33,] 0.0000000 0.0000000 0.0000000 0.0000000 0.7071068 -1.2247449
[34,] 0.0000000 0.0000000 0.0000000 0.0000000 0.7071068 -0.9525793
[35,] 0.0000000 0.0000000 0.0000000 0.0000000 0.7071068 -0.6804138
[36,] 0.0000000 0.0000000 0.0000000 0.0000000 0.7071068 -0.4082483
[37,] 0.0000000 0.0000000 0.0000000 0.0000000 0.7071068 -0.1360828
[38,] 0.0000000 0.0000000 0.0000000 0.0000000 0.7071068 0.1360828
[39,] 0.0000000 0.0000000 0.0000000 0.0000000 0.7071068 0.4082483
[40,] 0.0000000 0.0000000 0.0000000 0.0000000 0.7071068 0.6804138
[41,] 0.0000000 0.0000000 0.0000000 0.0000000 0.7071068 0.9525793
[42,] 0.0000000 0.0000000 0.0000000 0.0000000 0.7071068 1.2247449
[,24]
[1,] 0.0000000
[2,] 0.0000000
[3,] 0.0000000
[4,] 0.0000000
[5,] 0.0000000
[6,] 0.0000000
[7,] 0.0000000
[8,] 0.0000000
[9,] 0.0000000
[10,] 0.0000000
[11,] 0.0000000
[12,] 0.0000000
[13,] 0.0000000
[14,] 0.0000000
[15,] 0.0000000
[16,] 0.0000000
[17,] 0.0000000
[18,] 0.0000000
[19,] 0.0000000
[20,] 0.0000000
[21,] 0.0000000
[22,] 0.0000000
[23,] 0.0000000
[24,] 0.0000000
[25,] 0.0000000
[26,] 0.0000000
[27,] 0.0000000
[28,] 0.0000000
[29,] 0.0000000
[30,] 0.0000000
[31,] 0.0000000
[32,] 0.0000000
[33,] 1.5811388
[34,] 0.6441677
[35,] -0.0585607
[36,] -0.5270463
[37,] -0.7612891
[38,] -0.7612891
[39,] -0.5270463
[40,] -0.0585607
[41,] 0.6441677
[42,] 1.5811388
>
> # Q' R-1 Q 확인
> t(Q) %*% Q / R
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0 0 0 0 0 0 0 0 0 0.000000e+00
[2,] 0 0 0 0 0 0 0 0 0 0.000000e+00
[3,] 0 0 0 0 0 0 0 0 0 0.000000e+00
[4,] 0 0 0 0 0 0 0 0 0 0.000000e+00
[5,] 0 0 0 0 0 0 0 0 0 0.000000e+00
[6,] 0 0 0 0 0 0 0 0 0 0.000000e+00
[7,] 0 0 0 0 0 0 0 0 0 0.000000e+00
[8,] 0 0 0 0 0 0 0 0 0 0.000000e+00
[9,] 0 0 0 0 0 0 0 0 0 0.000000e+00
[10,] 0 0 0 0 0 0 0 0 0 1.347709e+00
[11,] 0 0 0 0 0 0 0 0 0 -5.103933e-17
[12,] 0 0 0 0 0 0 0 0 0 3.348410e-01
[13,] 0 0 0 0 0 0 0 0 0 0.000000e+00
[14,] 0 0 0 0 0 0 0 0 0 0.000000e+00
[15,] 0 0 0 0 0 0 0 0 0 0.000000e+00
[16,] 0 0 0 0 0 0 0 0 0 0.000000e+00
[17,] 0 0 0 0 0 0 0 0 0 0.000000e+00
[18,] 0 0 0 0 0 0 0 0 0 0.000000e+00
[19,] 0 0 0 0 0 0 0 0 0 0.000000e+00
[20,] 0 0 0 0 0 0 0 0 0 0.000000e+00
[21,] 0 0 0 0 0 0 0 0 0 0.000000e+00
[22,] 0 0 0 0 0 0 0 0 0 0.000000e+00
[23,] 0 0 0 0 0 0 0 0 0 0.000000e+00
[24,] 0 0 0 0 0 0 0 0 0 0.000000e+00
[,11] [,12] [,13] [,14] [,15]
[1,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[2,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[3,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[4,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[5,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[6,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[7,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[8,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[9,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[10,] -5.103933e-17 3.348410e-01 0.000000e+00 0.000000e+00 0.000000e+00
[11,] 1.647200e+00 -1.063746e-16 0.000000e+00 0.000000e+00 0.000000e+00
[12,] -1.063746e-16 2.035429e+00 0.000000e+00 0.000000e+00 0.000000e+00
[13,] 0.000000e+00 0.000000e+00 1.347709e+00 -5.103933e-17 3.348410e-01
[14,] 0.000000e+00 0.000000e+00 -5.103933e-17 1.647200e+00 -1.063746e-16
[15,] 0.000000e+00 0.000000e+00 3.348410e-01 -1.063746e-16 2.035429e+00
[16,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[17,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[18,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[19,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[20,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[21,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[22,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[23,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[24,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[,16] [,17] [,18] [,19] [,20] [,21]
[1,] 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[2,] 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[3,] 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[4,] 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[5,] 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[6,] 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[7,] 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[8,] 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[9,] 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[10,] 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[11,] 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[12,] 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[13,] 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[14,] 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[15,] 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[16,] 0.6738544 -0.6484167 0.1674205 0.00000000 0.0000000 0.00000000
[17,] -0.6484167 0.8235999 -0.5907016 0.00000000 0.0000000 0.00000000
[18,] 0.1674205 -0.5907016 1.0177143 0.00000000 0.0000000 0.00000000
[19,] 0.0000000 0.0000000 0.0000000 0.94339623 -0.5446701 -0.07812956
[20,] 0.0000000 0.0000000 0.0000000 -0.54467007 0.8735150 -0.67662183
[21,] 0.0000000 0.0000000 0.0000000 -0.07812956 -0.6766218 1.24880296
[22,] 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[23,] 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[24,] 0.0000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[,22] [,23] [,24]
[1,] 0.000000e+00 0.000000e+00 0.000000e+00
[2,] 0.000000e+00 0.000000e+00 0.000000e+00
[3,] 0.000000e+00 0.000000e+00 0.000000e+00
[4,] 0.000000e+00 0.000000e+00 0.000000e+00
[5,] 0.000000e+00 0.000000e+00 0.000000e+00
[6,] 0.000000e+00 0.000000e+00 0.000000e+00
[7,] 0.000000e+00 0.000000e+00 0.000000e+00
[8,] 0.000000e+00 0.000000e+00 0.000000e+00
[9,] 0.000000e+00 0.000000e+00 0.000000e+00
[10,] 0.000000e+00 0.000000e+00 0.000000e+00
[11,] 0.000000e+00 0.000000e+00 0.000000e+00
[12,] 0.000000e+00 0.000000e+00 0.000000e+00
[13,] 0.000000e+00 0.000000e+00 0.000000e+00
[14,] 0.000000e+00 0.000000e+00 0.000000e+00
[15,] 0.000000e+00 0.000000e+00 0.000000e+00
[16,] 0.000000e+00 0.000000e+00 0.000000e+00
[17,] 0.000000e+00 0.000000e+00 0.000000e+00
[18,] 0.000000e+00 0.000000e+00 0.000000e+00
[19,] 0.000000e+00 0.000000e+00 0.000000e+00
[20,] 0.000000e+00 0.000000e+00 0.000000e+00
[21,] 0.000000e+00 0.000000e+00 0.000000e+00
[22,] 1.347709e+00 -5.103933e-17 3.348410e-01
[23,] -5.103933e-17 1.647200e+00 -1.063746e-16
[24,] 3.348410e-01 -1.063746e-16 2.035429e+00
>
>
> # Z 생성
> Z <- Q[,10:24]
>
> # Z 확인
> Z
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.7071068 -1.2247449 1.5811388 0.0000000 0.0000000 0.0000000
[2,] 0.7071068 -0.9525793 0.6441677 0.0000000 0.0000000 0.0000000
[3,] 0.7071068 -0.6804138 -0.0585607 0.0000000 0.0000000 0.0000000
[4,] 0.7071068 -0.4082483 -0.5270463 0.0000000 0.0000000 0.0000000
[5,] 0.7071068 -0.1360828 -0.7612891 0.0000000 0.0000000 0.0000000
[6,] 0.7071068 0.1360828 -0.7612891 0.0000000 0.0000000 0.0000000
[7,] 0.7071068 0.4082483 -0.5270463 0.0000000 0.0000000 0.0000000
[8,] 0.7071068 0.6804138 -0.0585607 0.0000000 0.0000000 0.0000000
[9,] 0.7071068 0.9525793 0.6441677 0.0000000 0.0000000 0.0000000
[10,] 0.7071068 1.2247449 1.5811388 0.0000000 0.0000000 0.0000000
[11,] 0.0000000 0.0000000 0.0000000 0.7071068 -1.2247449 1.5811388
[12,] 0.0000000 0.0000000 0.0000000 0.7071068 -0.9525793 0.6441677
[13,] 0.0000000 0.0000000 0.0000000 0.7071068 -0.6804138 -0.0585607
[14,] 0.0000000 0.0000000 0.0000000 0.7071068 -0.4082483 -0.5270463
[15,] 0.0000000 0.0000000 0.0000000 0.7071068 -0.1360828 -0.7612891
[16,] 0.0000000 0.0000000 0.0000000 0.7071068 0.1360828 -0.7612891
[17,] 0.0000000 0.0000000 0.0000000 0.7071068 0.4082483 -0.5270463
[18,] 0.0000000 0.0000000 0.0000000 0.7071068 0.6804138 -0.0585607
[19,] 0.0000000 0.0000000 0.0000000 0.7071068 0.9525793 0.6441677
[20,] 0.0000000 0.0000000 0.0000000 0.7071068 1.2247449 1.5811388
[21,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[22,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[23,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[24,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[25,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[26,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[27,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[28,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[29,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[30,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[31,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[32,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[33,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[34,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[35,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[36,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[37,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[38,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[39,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[40,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[41,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[42,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[,7] [,8] [,9] [,10] [,11] [,12]
[1,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[2,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[3,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[4,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[5,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[6,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[7,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[8,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[9,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[10,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[11,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[12,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[13,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[14,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[15,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[16,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[17,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[18,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[19,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[20,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[21,] 0.7071068 -1.2247449 1.5811388 0.0000000 0.0000000 0.0000000
[22,] 0.7071068 -0.9525793 0.6441677 0.0000000 0.0000000 0.0000000
[23,] 0.7071068 -0.6804138 -0.0585607 0.0000000 0.0000000 0.0000000
[24,] 0.7071068 -0.4082483 -0.5270463 0.0000000 0.0000000 0.0000000
[25,] 0.7071068 -0.1360828 -0.7612891 0.0000000 0.0000000 0.0000000
[26,] 0.0000000 0.0000000 0.0000000 0.7071068 -1.2247449 1.5811388
[27,] 0.0000000 0.0000000 0.0000000 0.7071068 -0.9525793 0.6441677
[28,] 0.0000000 0.0000000 0.0000000 0.7071068 -0.6804138 -0.0585607
[29,] 0.0000000 0.0000000 0.0000000 0.7071068 -0.4082483 -0.5270463
[30,] 0.0000000 0.0000000 0.0000000 0.7071068 -0.1360828 -0.7612891
[31,] 0.0000000 0.0000000 0.0000000 0.7071068 0.1360828 -0.7612891
[32,] 0.0000000 0.0000000 0.0000000 0.7071068 0.4082483 -0.5270463
[33,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[34,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[35,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[36,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[37,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[38,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[39,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[40,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[41,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[42,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[,13] [,14] [,15]
[1,] 0.0000000 0.0000000 0.0000000
[2,] 0.0000000 0.0000000 0.0000000
[3,] 0.0000000 0.0000000 0.0000000
[4,] 0.0000000 0.0000000 0.0000000
[5,] 0.0000000 0.0000000 0.0000000
[6,] 0.0000000 0.0000000 0.0000000
[7,] 0.0000000 0.0000000 0.0000000
[8,] 0.0000000 0.0000000 0.0000000
[9,] 0.0000000 0.0000000 0.0000000
[10,] 0.0000000 0.0000000 0.0000000
[11,] 0.0000000 0.0000000 0.0000000
[12,] 0.0000000 0.0000000 0.0000000
[13,] 0.0000000 0.0000000 0.0000000
[14,] 0.0000000 0.0000000 0.0000000
[15,] 0.0000000 0.0000000 0.0000000
[16,] 0.0000000 0.0000000 0.0000000
[17,] 0.0000000 0.0000000 0.0000000
[18,] 0.0000000 0.0000000 0.0000000
[19,] 0.0000000 0.0000000 0.0000000
[20,] 0.0000000 0.0000000 0.0000000
[21,] 0.0000000 0.0000000 0.0000000
[22,] 0.0000000 0.0000000 0.0000000
[23,] 0.0000000 0.0000000 0.0000000
[24,] 0.0000000 0.0000000 0.0000000
[25,] 0.0000000 0.0000000 0.0000000
[26,] 0.0000000 0.0000000 0.0000000
[27,] 0.0000000 0.0000000 0.0000000
[28,] 0.0000000 0.0000000 0.0000000
[29,] 0.0000000 0.0000000 0.0000000
[30,] 0.0000000 0.0000000 0.0000000
[31,] 0.0000000 0.0000000 0.0000000
[32,] 0.0000000 0.0000000 0.0000000
[33,] 0.7071068 -1.2247449 1.5811388
[34,] 0.7071068 -0.9525793 0.6441677
[35,] 0.7071068 -0.6804138 -0.0585607
[36,] 0.7071068 -0.4082483 -0.5270463
[37,] 0.7071068 -0.1360828 -0.7612891
[38,] 0.7071068 0.1360828 -0.7612891
[39,] 0.7071068 0.4082483 -0.5270463
[40,] 0.7071068 0.6804138 -0.0585607
[41,] 0.7071068 0.9525793 0.6441677
[42,] 0.7071068 1.2247449 1.5811388
>
> # 혈통 입력
> pedi = as.matrix(read.table("strrm_pedi.txt"))
>
> # 혈통 확인
> pedi
V1 V2 V3
[1,] 1 0 0
[2,] 2 0 0
[3,] 3 0 0
[4,] 4 1 2
[5,] 5 3 2
[6,] 6 1 5
[7,] 7 3 4
[8,] 8 1 7
>
> # inverse of numerator relationship matrix
> Ainv = matrix(rep(0, 8 * 8), nrow = 8)
>
> Ainv
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 0 0 0 0 0 0 0 0
[2,] 0 0 0 0 0 0 0 0
[3,] 0 0 0 0 0 0 0 0
[4,] 0 0 0 0 0 0 0 0
[5,] 0 0 0 0 0 0 0 0
[6,] 0 0 0 0 0 0 0 0
[7,] 0 0 0 0 0 0 0 0
[8,] 0 0 0 0 0 0 0 0
>
> for ( i in 1:8) {
+ animal = pedi[i,1]
+ sire = pedi[i,2]
+ dam = pedi[i,3]
+
+ if (sire == 0 & dam == 0) {
+ alpha = 1
+ Ainv[animal,animal] = alpha + Ainv[animal,animal]
+ }
+ else if (sire != 0 & dam == 0) {
+ alpha = 4/3
+ Ainv[animal,animal] = alpha + Ainv[animal,animal]
+ Ainv[animal,sire] = -alpha/2 + Ainv[animal,sire]
+ Ainv[sire,animal] = -alpha/2 + Ainv[sire,animal]
+ Ainv[sire,sire] = alpha/4 + Ainv[sire,sire]
+ }
+ else if (sire == 0 & dam != 0) {
+ alpha = 4/3
+ Ainv[animal,animal] = alpha + Ainv[animal,animal]
+ Ainv[animal,dam] = -alpha/2 + Ainv[animal,dam]
+ Ainv[dam,animal] = -alpha/2 + Ainv[dam,animal]
+ Ainv[dam,dam] = alpha/4 + Ainv[dam,dam]
+ }
+ else {
+ alpha = 2
+ Ainv[animal,animal] = alpha + Ainv[animal,animal]
+ Ainv[animal,sire] = -alpha/2 + Ainv[animal,sire]
+ Ainv[sire,animal] = -alpha/2 + Ainv[sire,animal]
+ Ainv[animal,dam] = -alpha/2 + Ainv[animal,dam]
+ Ainv[dam,animal] = -alpha/2 + Ainv[dam,animal]
+ Ainv[sire,sire] = alpha/4 + Ainv[sire,sire]
+ Ainv[sire,dam] = alpha/4 + Ainv[sire,dam]
+ Ainv[dam,sire] = alpha/4 + Ainv[dam,sire]
+ Ainv[dam,dam] = alpha/4 + Ainv[dam,dam]
+ }
+ }
>
> # inverse of numerator relationship matrix 확인
> Ainv
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 2.5 0.5 0.0 -1.0 0.5 -1 0.5 -1
[2,] 0.5 2.0 0.5 -1.0 -1.0 0 0.0 0
[3,] 0.0 0.5 2.0 0.5 -1.0 0 -1.0 0
[4,] -1.0 -1.0 0.5 2.5 0.0 0 -1.0 0
[5,] 0.5 -1.0 -1.0 0.0 2.5 -1 0.0 0
[6,] -1.0 0.0 0.0 0.0 -1.0 2 0.0 0
[7,] 0.5 0.0 -1.0 -1.0 0.0 0 2.5 -1
[8,] -1.0 0.0 0.0 0.0 0.0 0 -1.0 2
>
> # LHS 생성
> LHS <- rbind(
+ cbind(t(X) %*% X / R, t(X) %*% Q / R, t(X) %*% Z / R)
+ ,
+ cbind(t(Q) %*% X / R, t(Q) %*% Q / R + Ainv %x% Ginv, t(Q) %*% Z / R)
+ ,
+ cbind(t(Z) %*% X / R, t(Z) %*% Q / R, t(Z) %*% Z / R + diag(5) %x% Pinv)
+ )
>
> # RHS 생성
> RHS <- rbind(
+ t(X) %*% y / R
+ ,
+ t(Q) %*% y / R
+ ,
+ t(Z) %*% y / R
+ )
>
> # 해 구하기
> strrm_sol = ginv(LHS) %*% RHS
>
> # 해 확인
> strrm_sol
[,1]
[1,] 7.831057464
[2,] 5.335708913
[3,] 6.305006624
[4,] 5.987899812
[5,] 4.060992098
[6,] 0.754996778
[7,] 0.853419772
[8,] 0.916711568
[9,] -1.750627724
[10,] -2.255001603
[11,] 19.827389899
[12,] -0.625390755
[13,] -0.134573003
[14,] 0.347902534
[15,] -0.421767141
[16,] -0.058312421
[17,] 0.055195979
[18,] -0.044183716
[19,] -0.072778891
[20,] -0.030488798
[21,] -0.024395667
[22,] 0.131091311
[23,] -0.024707181
[24,] 0.068579383
[25,] 0.344564931
[26,] 0.006282660
[27,] -0.316413050
[28,] -0.453733267
[29,] -0.052015857
[30,] 0.279819550
[31,] -0.548551607
[32,] 0.073007447
[33,] 0.194574056
[34,] 0.851808910
[35,] -0.009501574
[36,] -0.313065142
[37,] 0.220854001
[38,] 0.012696725
[39,] -0.017440925
[40,] -0.648601897
[41,] -0.360049925
[42,] -1.471825517
[43,] -0.776143491
[44,] 0.137006827
[45,] 0.968819056
[46,] -1.992680958
[47,] 0.985107123
[48,] -0.069307860
[49,] 3.518760564
[50,] -1.050969041
[51,] -0.404755224
[52,] -0.101334218
[53,] 0.288905015
[54,] 0.977069546
>

 

1307604800_strrm_R.zip

다운로드

파일