加权方差-协方差矩阵和 lapply答案

【问题标题】：Weighted variance-covariance matrices and lapply加权方差-协方差矩阵和 lapply
【发布时间】：2012-12-04 16:10:55
【问题描述】：

我有一个包含 50 个元素的列表 prob。每个元素都是一个 601x3 的概率矩阵，其中每一行代表一个完整的样本空间（即每个矩阵的每一行之和为 1）。例如，这里是prob 的第一个元素的前五行：

> prob[[1]][1:5,]

           [,1]      [,2]       [,3]
 [1,] 0.6027004 0.3655563 0.03174335
 [2,] 0.6013667 0.3665756 0.03205767
 [3,] 0.6000306 0.3675946 0.03237481
 [4,] 0.5986921 0.3686131 0.03269480
 [5,] 0.5973513 0.3696311 0.03301765

现在，我要做的是为列表prob 中的每个矩阵/元素的每一行创建以下矩阵。取第一行，设 a = .603、b = .366 和 c = .032（四舍五入到小数点后三位）。那么，

> w
         [,1]       [,2]       [,3]
 [1,] a*(1-a)       -a*b       -a*c
 [2,]    -b*a    b*(1-b)       -b*c
 [3,]    -c*a       -c*b    c*(1-c)

这样：

> w
           [,1]       [,2]       [,3]
 [1,]  0.239391  -0.220698  -0.019296
 [2,] -0.220698   0.232044  -0.011712
 [3,] -0.019296  -0.011712   0.030976

我想再获得一个类似的 3x3 矩阵 600 次（对于该矩阵的其余行），然后为 prob 的其余元素再重复整个过程 49 次。我唯一能想到的就是在lapply 内调用apply，这样我就可以一次访问每个矩阵的每一行。我敢肯定这不是一种优雅的方式（更不用说我无法让它工作），但我想不出其他任何东西。谁能帮我解决这个问题？我也很想听听关于使用不同结构的建议（例如，在列表中使用矩阵是不是很糟糕？）。

【问题讨论】：

601*50 矩阵是您正在寻找的结果吗？你知道你想如何存储它们（作为一个列表，在一个矩阵中）吗？或者也许你不需要存储它们。
我不太确定我在寻找什么类型的输出——我不能完全确定在这种情况下什么是最好的。他们在下面的阵列解决方案看起来很有希望。不过，我绝对想存储输出。我需要在创建这些 3x3 矩阵后立即使用它们。

标签： r covariance lapply

【解决方案1】：

使用lapply 在类似维度的矩阵列表上运行此过程应该非常简单。如果它代表一个挑战，那么您应该发布dput(.) 输出，以获得具有相似矩阵的两个元素列表。真正的挑战是逐行进行处理，如下所示，输出为 3x3xN 数组：

w <- apply(M, 1, function(rw) diag( rw*(1-rw) ) + 
                    rbind( rw*c(0, -rw[1], -rw[1] ), 
                           rw*c(-rw[2],0, -rw[2] ),
                           rw*c(-rw[3], -rw[3], 0)
         )

 )
 w
             [,1]        [,2]        [,3]        [,4]        [,5]
 [1,]  0.23945263  0.23972479  0.23999388  0.24025987  0.24052272
 [2,] -0.22032093 -0.22044636 -0.22056801 -0.22068575 -0.22079962
 [3,] -0.01913173 -0.01927842 -0.01942588 -0.01957412 -0.01972314
 [4,] -0.22032093 -0.22044636 -0.22056801 -0.22068575 -0.22079962
 [5,]  0.23192489  0.23219793  0.23246881  0.23273748  0.23300395
 [6,] -0.01160398 -0.01175156 -0.01190081 -0.01205173 -0.01220435
 [7,] -0.01913173 -0.01927842 -0.01942588 -0.01957412 -0.01972314
 [8,] -0.01160398 -0.01175156 -0.01190081 -0.01205173 -0.01220435
 [9,]  0.03073571  0.03102998  0.03132668  0.03162585  0.03192748

 w <- array(w, c(3,3,5) )
 w
, , 1

            [,1]        [,2]        [,3]
[1,]  0.23945263 -0.22032093 -0.01913173
[2,] -0.22032093  0.23192489 -0.01160398
[3,] -0.01913173 -0.01160398  0.03073571

, , 2

            [,1]        [,2]        [,3]
[1,]  0.23972479 -0.22044636 -0.01927842
[2,] -0.22044636  0.23219793 -0.01175156
[3,] -0.01927842 -0.01175156  0.03102998

.... snipped remaining output

【讨论】：

感谢您使用数组@DWin 的建议——我想我最终会走这条路。我对您的示例感到有些困惑，但这很可能只是我对 R 的有限了解。不过，我认为这足以为我指明正确的方向，所以谢谢！
很抱歉让您感到困惑。基本上，我只是将一个对角矩阵添加到一个矩阵中，其中每行包含您的规范的非对角线元素。这给了我一个 9 行矩阵，每一列都是由你的输入的一行构成的。因为矩阵和数组在 R 中是列优先的，所以您可以将这些列中的每一列转换为数组的 3 x 3“切片”，该数组的切片数与输入的行数一样多。使用 w[ , , 1] 一次得到一个切片； w[ ,, 2] 等
哦，我明白了——这很有道理！感谢您的知识！