matlab - 在matlab中计算1000个5x5矩阵的协方差

Question

我有一个像这样的 1000 个 5x5 矩阵 (Xm):

每个 $(x_ij)m$ 都是从分布中提取的点估计。我想计算每个 $x{ij}$ 的协方差 cov，其中 i=1..n，j=1..n 在红色方向箭头。

例如，$X_m$ 的方差是 `var(X,0,3)，它给出了一个 5x5 方差矩阵。我可以用同样的方法计算协方差吗？

尝试回答

到目前为止，我已经这样做了:

for m=1:1000
Xm_new(m,:)=reshape(Xm(:,:,m)',25,1);
end

cov(Xm_new)
spy(Xm_new) gives me this unusual looking sparse matrix:

Answer 1

如果您查看cov(在命令窗口中edit cov)，您可能会明白为什么它不支持多维数组。它执行输入矩阵的转置和矩阵乘法:xc' * xc。这两种操作都不支持多维数组，我想无论是谁编写该函数，都决定不做推广它的工作(但是联系 Mathworks 和 make a feature request 可能还是不错的)。

在您的情况下，如果我们从 cov 中获取基本代码并做出一些假设，我们可以编写支持 3-D 数组的协方差函数 M 文件:

function x = cov3d(x)
% Based on Matlab's cov, version 5.16.4.10

[m,n,p] = size(x);
if m == 1
    x = zeros(n,n,p,class(x));
else
    x = bsxfun(@minus,x,sum(x,1)/m);
    for i = 1:p
        xi = x(:,:,i);
        x(:,:,i) = xi'*xi;
    end
    x = x/(m-1);
end

请注意，此简单代码假定 x 是一系列沿三维堆叠的二维矩阵。规范化标志为 0，cov 中的默认值。它可以通过一些工作扩展到多个维度，如 var。在我的计时中，它比在 for 循环中调用 cov(x(:,:,i)) 的函数快 10 倍以上。

是的，我使用了 for 循环。可能有也可能没有 faster ways to do this ，但在这种情况下，for 循环将为 faster than most schemes ，尤其是当你的数组的大小不是先验的时候。