分别提取R中矩阵的上下三角形的元素答案

【问题标题】：Extract respective upper and lower triangle's elements of matrix in R分别提取R中矩阵的上下三角形的元素
【发布时间】：2016-07-21 12:50:46
【问题描述】：

给定方阵m如下（nxn）：

m <- matrix(1:5,ncol = 5,nrow = 5,byrow = F)

     [,1] [,2] [,3] [,4] [,5]
[1,]    1    1    1    1    1
[2,]    2    2    2    2    2
[3,]    3    3    3    3    3
[4,]    4    4    4    4    4
[5,]    5    5    5    5    5

我想分别提取上下三角形的元素并取相对频率。

我们可以通过这样的循环天真地做到这一点（这里是n=5）：

for (i in 1:(n-1))
    for (j in (i+1):n){
        x <- m[i,j]
        y <- m[j,i]
        m[i,j] <- x/(x+y)
        m[j,i] <- y/(x+y)
    }

这是所需的输出：

          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 1.0000000 0.3333333 0.2500000 0.2000000 0.1666667
[2,] 0.6666667 2.0000000 0.4000000 0.3333333 0.2857143
[3,] 0.7500000 0.6000000 3.0000000 0.4285714 0.3750000
[4,] 0.8000000 0.6666667 0.5714286 4.0000000 0.4444444
[5,] 0.8333333 0.7142857 0.6250000 0.5555556 5.0000000

我们能否更有效地生成此输出？

附言

我知道m[upper.tri(m)] 和m[lower.tri(m)]，但这不是诀窍，因为提取元素的顺序不同。例如，m[upper.tri(m)] 会给我：

[1] 1 1 2 1 2 3 1 2 3 4

而我想要的上三角形是：

[1] 1 1 1 1 2 2 2 3 3 4

【问题讨论】：

@AlexIoannides 删除的答案提供了线索 - ?upper.tri 和 ?lower.tri 构成了解决方案的基础......
@BenBolker 我知道，但这不是诀窍。
@m0h3n 看看这个帖子github.com/gastonstat/matrixkit/blob/master/R/…
你可以通过t()和upper.tri/lower.tri的适当组合来完成它，我认为......
如果你提供一些测试数据，我可以看看。我会尝试upper.tri(m)/(m+t(m))

标签： r matrix

【解决方案1】：

更容易：

f <- m/(m+t(m))
#> f
#          [,1]      [,2]      [,3]      [,4]      [,5]
#[1,] 0.5000000 0.3333333 0.2500000 0.2000000 0.1666667
#[2,] 0.6666667 0.5000000 0.4000000 0.3333333 0.2857143
#[3,] 0.7500000 0.6000000 0.5000000 0.4285714 0.3750000
#[4,] 0.8000000 0.6666667 0.5714286 0.5000000 0.4444444
#[5,] 0.8333333 0.7142857 0.6250000 0.5555556 0.5000000

即使是对角线也是这样计算的，但是没有任何信息，所以使用diag(f) = diag(m)得到想要的输出。

#> f = m/(m+t(m))
#> diag(f) = diag(m)
#> f
#          [,1]      [,2]      [,3]      [,4]      [,5]
#[1,] 1.0000000 0.3333333 0.2500000 0.2000000 0.1666667
#[2,] 0.6666667 2.0000000 0.4000000 0.3333333 0.2857143
#[3,] 0.7500000 0.6000000 3.0000000 0.4285714 0.3750000
#[4,] 0.8000000 0.6666667 0.5714286 4.0000000 0.4444444
#[5,] 0.8333333 0.7142857 0.6250000 0.5555556 5.0000000

【讨论】：

【解决方案2】：

将lower.tri 与转置t() 结合起来，我们可以按照您想要的顺序得到upper.tri

t(m)[lower.tri(t(m))]
#[1] 1 1 1 1 2 2 2 3 3 4

为了获得所需的输出，我们使用相同的策略

#lower.tri and upper.tri matrices
lt <- m[lower.tri(m)]
ut <- t(m)[lower.tri(t(m))]
#defining the frequency matrix, just so it has the same dimensions of m
f <- m
#we can't assing a value to t(f), so we assing the upper.tri frequency to the lower.tri portion of f, then transpose f
f[lower.tri(f)] <- ut/(lt+ut) 
f <- t(f)
#then the lower.tri portion follows
f[lower.tri(f)] <- lt/(lt+ut)
#> f
#          [,1]      [,2]      [,3]      [,4]      [,5]
#[1,] 1.0000000 0.3333333 0.2500000 0.2000000 0.1666667
#[2,] 0.6666667 2.0000000 0.4000000 0.3333333 0.2857143
#[3,] 0.7500000 0.6000000 3.0000000 0.4285714 0.3750000
#[4,] 0.8000000 0.6666667 0.5714286 4.0000000 0.4444444
#[5,] 0.8333333 0.7142857 0.6250000 0.5555556 5.0000000

【讨论】：

【解决方案3】：

这是使用combn 的另一种解决方案（它有助于按所需顺序提取上三角形的元素。）：

out <- m
inds <- t(combn(ncol(m),2))

 # > inds
      # [,1] [,2]
 # [1,]    1    2
 # [2,]    1    3
 # [3,]    1    4
 # [4,]    1    5
 # [5,]    2    3
 # [6,]    2    4
 # [7,]    2    5
 # [8,]    3    4
 # [9,]    3    5
# [10,]    4    5

denom <- m[inds]+m[inds[,2:1]]

out[inds] <- m[inds]/denom
out[inds[,2:1]] <- m[inds[,2:1]]/denom

out
          # [,1]      [,2]      [,3]      [,4]      [,5]
# [1,] 1.0000000 0.3333333 0.2500000 0.2000000 0.1666667
# [2,] 0.6666667 2.0000000 0.4000000 0.3333333 0.2857143
# [3,] 0.7500000 0.6000000 3.0000000 0.4285714 0.3750000
# [4,] 0.8000000 0.6666667 0.5714286 4.0000000 0.4444444
# [5,] 0.8333333 0.7142857 0.6250000 0.5555556 5.0000000

基准测试

library(microbenchmark)
m <- matrix(1:5,ncol = 5,nrow = 5,byrow = F)

f_m0h3n1 <- function(m){
    n <- ncol(m)
    for (i in 1:(n-1))
        for (j in (i+1):n){x <- m[i,j];y <- m[j,i];m[i,j] <- x/(x+y);m[j,i] <- y/(x+y);}
    return(m)
}

f_catastrophic_failure1 <- function(m){f <- m/(m+t(m));diag(f) <- diag(m);return(f);}

f_m0h3n2 <- function(m){out <- m;inds <- t(combn(ncol(m),2));denom <- m[inds]+m[inds[,2:1]];out[inds] <- m[inds]/denom;out[inds[,2:1]] <- m[inds[,2:1]]/denom;return(out);}

f_catastrophic_failure2 <- function(m){
    lt <- m[lower.tri(m)];ut <- t(m)[lower.tri(t(m))];f <- m;f[lower.tri(f)] <- ut/(lt+ut);f <- t(f);f[lower.tri(f)] <- lt/(lt+ut);return(f);
}
r <- f_m0h3n1(m)
all(f_m0h3n2(m) == r)
# [1] TRUE
all(f_catastrophic_failure1(m) == r)
# [1] TRUE
all(f_catastrophic_failure2(m) == r)
# [1] TRUE

microbenchmark(f_m0h3n1(m), f_m0h3n2(m), f_catastrophic_failure1(m), f_catastrophic_failure2(m))

# Unit: microseconds
                       # expr    min      lq      mean median     uq     max neval
                # f_m0h3n1(m) 70.575 73.7825  78.95802 75.707 83.407 131.312   100
                # f_m0h3n2(m) 87.255 90.2500  96.36626 91.533 98.163 244.230   100
 # f_catastrophic_failure1(m) 33.790 35.5010  38.37556 36.785 37.640 142.432   100
 # f_catastrophic_failure2(m) 91.961 95.3825 102.27319 97.735 99.660 303.256   100

【讨论】：