3

我有一个像这样的 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:

在此处输入图像描述

4

2 回答 2

5

如果您查看cov（edit cov在命令窗口中），您可能会明白为什么它不支持多维数组。它执行输入矩阵的转置和矩阵乘法：xc' * xc. 这两种操作都不支持多维数组，我猜想编写函数的人决定不做推广它的工作（但是联系 Mathworks 并提出功能请求仍然可能很好）。

在您的情况下，如果我们从中获取基本代码cov并做出一些假设，我们可以编写一个协方差函数 M-file 支持 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做一些工作。在我的时间里，它比循环调用cov(x(:,:,i))的函数快 10 倍以上。for

是的，我使用了一个for循环。可能有也可能没有更快的方法来做到这一点，但在这种情况下，for循环将比大多数方案更快，特别是当你的数组大小不知道先验时。

于 2013-07-02T21:43:57.377 回答
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）。

于 2014-09-12T20:06:31.500 回答