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

Question

我有1000个5x5矩阵(Xm),如下所示:

每个$(x_ij)m $是从分布中得出的点估计.我想计算cov每个$ x {ij} $ 的协方差,其中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并提出功能请求).

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

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工作来扩展到多个维度.在我的时序,它比调用一个函数快10倍以上cov(x(:,:,i))的for循环.

是的,我使用了一个for循环.可能有或没有更快的方法来执行此操作,但在这种情况下,for循环将比大多数方案更快,特别是当您的阵列的大小不是先验已知时.