将现有 Cov 矩阵转换为块对角线

Question

我有一个现有的协方差矩阵，我想根据各个列所属的组将其转换为块对角线（例如，第 2 行/列是第 1 组，接下来是第 2 组等）是否有一种简单的方法这样做：

下面是我所拥有的示例：

m1 <- matrix(1:16, ncol=4, byrow=TRUE)
rownames(m1) <- colnames(m1 ) <- c('a', 'b', 'c', 'd')

   a  b  c  d
a  1  2  3  4
b  5  6  7  8
c  9 10 11 12
d 13 14 15 16

我有2组：

第 1 组：“a”、“b”

第 2 组：“c”、“d”

我想要什么：

   a  b  c  d
a  1  2  0  0
b  5  6  0  0
c  0  0 11 12
d  0  0 15 16

Answer 1

使用for循环。

g <- list(c('a', 'b'), c('c', 'd'))
for (x in g) m1[!rownames(m1) %in% x, colnames(m1) %in% x] <- 0
m1
#   a b  c  d
# a 1 2  0  0
# b 5 6  0  0
# c 0 0 11 12
# d 0 0 15 16

Answer 2

我们可以使用以下which函数来使用逻辑索引：which我们使用逻辑索引将每个组之外的元素设置为 0。所以我们不需要循环。这在大型矩阵中可能更有效。

group1 <- c('a', 'b')
group2 <- c('c', 'd')

# logical indices for all groups
idx_group1 <- which(colnames(m1) %in% group1)
idx_group2 <- which(colnames(m1) %in% group2)


m1[-idx_group1, idx_group1] <- 0
m1[idx_group1, -idx_group1] <- 0
m1[-idx_group2, idx_group2] <- 0
m1[idx_group2, -idx_group2] <- 0

m1

  a b  c  d
a 1 2  0  0
b 5 6  0  0
c 0 0 11 12
d 0 0 15 16

Answer 3

这是一个方法。

fun <- function(x, groups) {
  y <- matrix(0, nrow(x), ncol(x), dimnames = dimnames(x))
  for(i in seq_along(groups)) 
    y[groups[[i]], groups[[i]]] <- 1L
  x * y
}

m1 <- matrix(1:16, ncol=4, byrow=TRUE)
rownames(m1) <- colnames(m1 ) <- c('a', 'b', 'c', 'd')
group1 <- c("a", "b")
group2 <- c("c", "d")

fun(m1, list(group1, group2))
#>   a b  c  d
#> a 1 2  0  0
#> b 5 6  0  0
#> c 0 0 11 12
#> d 0 0 15 16

^{创建于 2023-03-17，使用reprex v2.0.2}

Answer 4

我们可以使用++tcrossprod创建一个掩码矩阵，例如，tablestack

> tcrossprod(table(stack(setNames(g, seq_along(g)))))
      values
values a b c d
     a 1 1 0 0
     b 1 1 0 0
     c 0 0 1 1
     d 0 0 1 1

这样

> m1 * tcrossprod(table(stack(setNames(g, seq_along(g)))))
  a b  c  d
a 1 2  0  0
b 5 6  0  0
c 0 0 11 12
d 0 0 15 16

在哪里

g <- list(c("a", "b"), c("c", "d"))

Answer 5

您可以使用linpk带有功能的包blockdiag，其中您可以使用所需组的两个子集，如下所示：

m1 <- matrix(1:16, ncol=4, byrow=TRUE)
rownames(m1) <- colnames(m1 ) <- c('a', 'b', 'c', 'd')
library(linpk)
blockdiag(m1[c("a", "b"), c("a", "b")], m1[c("c", "d"), c("c", "d")])
#>   a b  c  d
#> a 1 2  0  0
#> b 5 6  0  0
#> c 0 0 11 12
#> d 0 0 15 16

^{创建于 2023-03-17，使用reprex v2.0.2}