計算1000點5×5矩陣的MATLAB中的協方差

我有1000點5×5矩陣（XM）所示：計算1000點5×5矩陣的MATLAB中的協方差

enter image description here

每個$（x_ij）米$是從繪製的點估計分配。我想計算每個$ x {ij} $的協方差cov，其中i = 1..n，j = 1..n在紅色箭頭的方向上。

例如$ X_m $的方差是var（X，0,3），它給出了一個5×5的方差矩陣。我能以相同的方式計算協方差嗎？

嘗試在回答

到目前爲止，我已經做到了這一點：

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:

enter image description here

來源

2013-07-02 HCAI

如果你看看cov（edit cov在命令窗口）你可能會看到它爲什麼不支持多維數組。它執行輸入矩陣的轉置和矩陣乘法：xc' * xc。這兩個操作都不支持多維數組，我猜想編寫該函數的人決定不做這個工作來推廣它（然而，聯繫Mathworks可能會很好，但也可能是make a feature request）。

在你的情況，如果我們從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是沿着第三維疊加的一系列2-D矩陣。標準化標誌爲0，默認值爲cov。它可以用var等多個維度來完成，並且有一定的工作量。在我的計時中，它的速度比在for循環中調用cov(x(:,:,i))的函數快10倍。

是的，我用了一個for循環。可能有也可能不是faster ways to do this，但在這種情況下，for循環將是faster than most schemes，尤其是當您的數組的大小事先並不知道時。

來源

2013-07-02 21:43:57 horchler

啊非常好！謝謝！然而，我正在尋找基於每個變量x_ij的1000個觀察值（在第三維中彼此堆疊在一起）的5×5協方差矩陣。因此，例如Xm（1,1，:)是同一個變量的觀察結果等等。這個函數在做什麼？ – HCAI

它相當於爲一個5乘5的輸入調用'cov'函數1000次並存儲5乘5的輸出：'Xm_out（：，：1）= cov（Xm（：，：1）），Xm_out（：，：，2）= cov（Xm（：，：，2））'，...，'Xm_out ）'，'Xm_out（：，：，1000）= cov（Xm（：，：，1000））''。 – horchler

唉啊...我試圖找出每個x_ij，i \ ne = j的協方差，這樣x12是列向量x1和x2之間沿着紅色箭頭方向的協方差。我沒問過正確的問題嗎？ – HCAI

下面還答案適用於長方形矩陣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-12 20:06:31 michael

計算1000點5×5矩陣的MATLAB中的協方差

回答

相關問題