Reputation: 93
How could I build a function that extracts the diagonal blocks matrices of a larger one? The problem is as follows. The function takes a centred matrix as argument, computes the full error covariance matrix and extracts the blocks on the leading diagonal? I tried the following, but not working.
err_cov <- function(x){
m <- nrow(x)
n <- ncol(x)
#compute the full error covariance matrix as the inner product
#of vec(x) and its transpose. Note that, omega is a mnxmn matrix
vec <- as.vector(x)
omega <- vec%*%t(vec)
sigmas <- list()
for(i in 0:n-1){
#here the blocks have to be m nxn matrices along the
#leading diagonal
for (j in 1:m)
sigmas[[j]] <- omega[(n*i+1):n*(i+1), (n*i+1):n*(i+1)]
}
return(sigmas)
}
So, for instance for
A [,1] [,2] [,3] [,4] [1,] 1 4 7 10 [2,] 2 5 8 11 [3,] 3 6 9 12
> B<-as.vector(A)
> B
[1] 1 2 3 4 5 6 7 8 9 10 11 12
> C<-B%*%t(B)
> C
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,] 1 2 3 4 5 6 7 8 9 10 11 12
[2,] 2 4 6 8 10 12 14 16 18 20 22 24
[3,] 3 6 9 12 15 18 21 24 27 30 33 36
[4,] 4 8 12 16 20 24 28 32 36 40 44 48
[5,] 5 10 15 20 25 30 35 40 45 50 55 60
[6,] 6 12 18 24 30 36 42 48 54 60 66 72
[7,] 7 14 21 28 35 42 49 56 63 70 77 84
[8,] 8 16 24 32 40 48 56 64 72 80 88 96
[9,] 9 18 27 36 45 54 63 72 81 90 99 108
[10,] 10 20 30 40 50 60 70 80 90 100 110 120
[11,] 11 22 33 44 55 66 77 88 99 110 121 132
[12,] 12 24 36 48 60 72 84 96 108 120 132 144
The function should return:
> C1
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 2 4 6
[3,] 3 6 9
> C2
[,1] [,2] [,3]
[1,] 16 20 24
[2,] 20 25 30
[3,] 24 30 36
> C3
[,1] [,2] [,3]
[1,] 49 56 63
[2,] 56 64 72
[3,] 63 72 81
> C4
[,1] [,2] [,3]
[1,] 100 110 120
[2,] 110 121 132
[3,] 120 132 144
Thanks for answering.
Upvotes: 1
Views: 152
Reputation: 9687
I think a clearer solution is to reset the dimensions and then let R do the index calculations for you:
err_cov <- function(x){
m <- nrow(x)
n <- ncol(x)
#compute the full error covariance matrix as the inner product
#of vec(x) and its transpose
vec <- as.vector(x)
omega <- tcrossprod(vec)
dim(omega) <- c(n,m,n,m)
sigmas <- list()
for (j in 1:m)
sigmas[[j]] <- omega[,j,,j]
return(sigmas)
}
Here is an example:
> x
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
> tcrossprod(vec)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 2 3 4 5 6
[2,] 2 4 6 8 10 12
[3,] 3 6 9 12 15 18
[4,] 4 8 12 16 20 24
[5,] 5 10 15 20 25 30
[6,] 6 12 18 24 30 36
> err_cov(x)
[[1]]
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 2 4 6
[3,] 3 6 9
[[2]]
[,1] [,2] [,3]
[1,] 16 20 24
[2,] 20 25 30
[3,] 24 30 36
Upvotes: 1