如何計算MATLAB中沒有任何內置函數或循環的協方差矩陣？

Question

是否可以在MATLAB中使用任何內置函數或循環來查找矩陣的協方差？我對解決這個問題的想法完全無能爲力。

我想的是一樣的東西：

cov(x,y) = 1/(n-1) .* (x*y) 

不過，我不認爲這會工作。有任何想法嗎？

Answer 1

下面是如何數值計算協方差矩陣的一個很好的例子。 http://www.itl.nist.gov/div898/handbook/pmc/section5/pmc541.htm。但是，爲了完整起見，我們在這篇文章中提到這個。 我對你「內置」函數的含義有點困惑，因爲協方差要求你總結一個矩陣的列。如果您不能使用任何內置函數來總結這些元素，那麼我不知道如何在不使用for循環的情況下做到這一點。編輯：我想出瞭如何在不使用內置函數或循環的情況下做到這一點，但您需要使用size來確定矩陣中有多少行......除非您在函數中將此指定爲常數。

數字上，你計算的協方差矩陣，像這樣：

從本質上講，我^個行j ^{日的協方差矩陣的}列是這樣的，你拿的總和列i的產物減去列i的列的平均值，列j減去列j的平均值。現在，加上這些，然後除以n - 1。這被稱爲unbiased estimator。您還會注意到這個矩陣是對稱的，因爲即使您翻轉訂單（即查看第j列，然後查看第i列），答案仍然應該是相同的。我假設你不能使用MATLAB的mean，所以讓我們從第一原理開始。

首先，計算一個計算每列的平均值的行向量。如果不使用sum（因爲它也是內置函數），您可以執行的計算所有列的總和的方法是將此行向量與您的矩陣A相乘，輸出將是包含總和的行向量在所有列上。因此，這樣做：

one_vector(1:size(A,1)) = 1; 
mu = (one_vector * A)/size(A,1); 

與代碼的第一行中的特技是，我們動態地創建一個數組，它是在基質A相同的長度的行數的。我們填寫完整的1。請注意，你可以使用ones，但你說你不能使用任何內置函數。 mu將在所有列中包含我們的矢量。

現在，讓我們通過用平均值減去每列來預處理數據，因爲這就是我們所做的定義。要做到這一點沒有任何內置函數，你可以做的是用各自的方法減去所有列，重複mu的次數與one_vector中的1次一樣多。因此：

A_mean_subtract = A - mu(one_vector, :); 

這裏有點棘手（和很酷）。如果我們對矩陣A進行轉置，您會看到行成爲列，列成爲行。如果我們採用這個轉置並乘以原始矩陣，我們實際上將獲得矩陣A的列i和列j之間的乘積之和。這是我們協方差計算的第一部分。然後我們除以n - 1。因此，我們的協方差很簡單：

covA = (A_mean_subtract.' * A_mean_subtract)/(size(A,1) - 1); 

下面是一個快速示例，以及我在上面向您展示的網站上看到的內容。假如A是這樣的：

A = [4 2 0.5; 4.2 2.1 0.59; 3.9 2.0 0.58; 4.3 2.1 0.62; 4.1 2.2 0.63] 

A = 

    4.0000 2.0000 0.5000 
    4.2000 2.1000 0.5900 
    3.9000 2.0000 0.5800 
    4.3000 2.1000 0.6200 
    4.1000 2.2000 0.6300 

通過上面的代碼運行，這就是我們得到：

covA = 

    0.0250 0.0075 0.0042 
    0.0075 0.0070 0.0034 
    0.0042 0.0034 0.0026 

你會看到，這也與MATLAB中的cov功能也很相配：

>> cov(A) 

ans = 

0.0250 0.0075 0.0042 
0.0075 0.0070 0.0034 
0.0042 0.0034 0.0026 

位提示

如果您輸入edit cov在MATLAB命令提示符下，你可以看到他們是如何計算的協方差矩陣沒有任何for循環....這基本上是我給了同樣的答案你:)

如果你想這樣做，更有效地

假設您可以使用sum和bsxfun，我們可以使用更少（更有效率）的代碼行來完成此操作。首先，計算你的平均向量就像我們上面使用sum：

mu = sum(A)/size(A,1); 

現在，每列的相應的平均減去你的矩陣A，你可以用bsxfun來幫助你促進這一減法：

A_mean_subtract = bsxfun(@minus, A, mu); 

現在，計算你的協方差矩陣像之前一樣：

covA = (A_mean_subtract.' * A_mean_subtract)/(size(A,1) - 1); 

你應該得到確切的結果爲w相同以前看過。對穩定性

我們使用的計算使用的定義兩列之間的協方差的直線上升定義

小注。但是，已經表明，如果您提供某些類型的數據，那麼使用直接定義可能會導致數值不穩定。請參閱this Wikipedia page，其中介紹了計算兩個更穩定的長度向量之間的協方差的各種算法。