使用Matlab计算1000个5x5矩阵的协方差。

Question

使用Matlab计算1000个5x5矩阵的协方差。

3

我有1000个类似于这样的5x5矩阵（Xm）：

enter image description here

每个$(x_ij)m$都是从一个分布中得出的点估计。我想计算每个$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:

- HCAI

2个回答

2

下面的答案也适用于矩形矩阵 xi=x(:,:,i)。

function xy = cov3d(x)

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

我的答案与horchler非常相似，但是horchler的代码不能用于矩形矩阵 xi（其尺寸与xi'*xi的尺寸不同）。

- michael

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- horchler · Accepted Answer

如果您查看cov（在命令窗口中编辑cov），您可能会看到为什么它不支持多维数组。它执行输入矩阵的转置和矩阵乘法：xc' * xc。这两个操作都不支持多维数组，我猜想编写该函数的人决定不去通用化它（但是联系Mathworks并提出特性请求仍然是一个好主意）。

在您的情况下，如果我们从cov中取基本代码并进行一些假设，我们可以编写一个支持三维数组的协方差函数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是沿第三维度堆叠的一系列2-D矩阵。并且归一化标志为0，这是cov中的默认值。它可以像var那样扩展到多个维度，需要进行一些工作。在我的时间测试中，它比一个使用for循环调用cov(x(:,:,i))的函数快了10倍以上。

是的，我使用了for循环。可能有更快的方法来完成这个任务，但在这种情况下，for循环将会比大多数方案更快,尤其是当你的数组大小事先未知时。