第一部分-協方差的意義和計算公式

學過概率統計的孩子都知道，統計裡最基本的概念就是樣本的均值，方差，或者再加個標準差。首先我們給你一個含有n個樣本的集合，依次給出這些概念的公式描述，這些高中學過數學的孩子都應該知道吧，一帶而過。

很顯然，均值描述的是樣本集合的中間點，它告訴我們的資訊是很有限的，而標準差給我們描述的則是樣本集合的各個樣本點到均值的距離之平均。以這兩個集合為例，[0，8，12，20]和[8，9，11，12]，兩個集合的均值都是10，但顯然兩個集合差別是很大的，計算兩者的標準差，前者是8.3，後者是1.8，顯然後者較為集中，故其標準差小一些，標準差描述的就是這種“散佈度”。之所以除以n-1而不是除以n，是因為這樣能使我們以較小的樣本集更好的逼近總體的標準差，即統計上所謂的“無偏估計”。而方差則僅僅是標準差的平方。

為什麼需要協方差？

上面幾個統計量看似已經描述的差不多了，但我們應該注意到，標準差和方差一般是用來描述一維資料的，但現實生活我們常常遇到含有多維資料的資料集，最簡單的大家上學時免不了要統計多個學科的考試成績。面對這樣的資料集，我們當然可以按照每一維獨立的計算其方差，但是通常我們還想了解更多，比如，一個男孩子的猥瑣程度跟他受女孩子歡迎程度是否存在一些聯絡啊，嘿嘿~協方差就是這樣一種用來度量兩個隨機變數關係的統計量，我們可以仿照方差的定義：

來度量各個維度偏離其均值的程度，標準差可以這麼來定義：

協方差的結果有什麼意義呢？如果結果為正值，則說明兩者是正相關的(從協方差可以引出“相關係數”的定義)，也就是說一個人越猥瑣就越受女孩子歡迎，嘿嘿，那必須的~結果為負值就說明負相關的，越猥瑣女孩子越討厭，可能嗎？如果為0，也是就是統計上說的“相互獨立”。

從協方差的定義上我們也可以看出一些顯而易見的性質，如：

協方差多了就是協方差矩陣

上一節提到的猥瑣和受歡迎的問題是典型二維問題，而協方差也只能處理二維問題，那維數多了自然就需要計算多個協方差，比如n維的資料集就需要計算 n! / ((n-2)!*2) 個協方差，那自然而然的我們會想到使用矩陣來組織這些資料。給出協方差矩陣的定義：

這個定義還是很容易理解的，我們可以舉一個簡單的三維的例子，假設資料集有三個維度，則協方差矩陣為

可見，協方差矩陣是一個對稱的矩陣，而且對角線是各個維度上的方差。

Matlab協方差實戰

上面涉及的內容都比較容易，協方差矩陣似乎也很簡單，但實戰起來就很容易讓人迷茫了。必須要明確一點，協方差矩陣計算的是不同維度之間的協方差，而不是不同樣本之間的。

首先，隨機產生一個10*3維的整數矩陣作為樣本集，10為樣本的個數，3為樣本的維數。

mysample = fix(rand(10,3)*50)

根據公式，計算協方差需要計算均值，那是按行計算均值還是按列呢，我一開始就老是困擾這個問題。前面我們也特別強調了，協方差矩陣是計算不同維度間的協方差，要時刻牢記這一點。樣本矩陣的每行是一個樣本，每列為一個維度，所以我們要按列計算均值。為了描述方便，我們先將三個維度的資料分別賦值：

>> dim1 = mysample(:,1);
>> dim2 = mysample(:,2);
>> dim3 = mysample(:,3);

計算dim1與dim2，dim1與dim3，dim2與dim3的協方差：

>> sum((dim1 - mean(dim1)) .* (dim2 - mean(dim2))) / (size(mysample, 1) - 1)  %得到 -147.0667
>> sum((dim1 - mean(dim1)) .* (dim3 - mean(dim3))) / (size(mysample, 1) - 1)  %得到  -82.2667
>> sum((dim2 - mean(dim2)) .* (dim3 - mean(dim3))) / (size(mysample, 1) - 1)  %得到   76.5111

搞清楚了這個後面就容易多了，協方差矩陣的對角線就是各個維度上的方差，下面我們依次計算：

>> var(dim1)  %得到 227.8778
>> var(dim2)  %得到 179.8222
>> var(dim3)  %得到 156.7111

 這樣，我們就得到了計算協方差矩陣所需要的所有資料，呼叫Matlab自帶的cov函式進行驗證：

>> cov(mysample)

 把我們計算的資料對號入座，是不是一摸一樣？

Update

今天突然發現，原來協方差矩陣還可以這樣計算，先讓樣本矩陣中心化，即每一維度減去該維度的均值，使每一維度上的均值為0，然後直接用新的到的樣本矩陣乘上它的轉置，然後除以(N-1)即可。其實這種方法也是由前面的公式推導而來，只不過理解起來不是很直觀，但在抽象的公式推導時還是很常用的！同樣給出Matlab程式碼實現：

>> temp = mysample - repmat(mean(mysample), 10, 1);
>> result = temp' * temp ./ (size(mysample, 1) - 1)

總結

理解協方差矩陣的關鍵就在於牢記它計算的是不同維度之間的協方差，而不是不同樣本之間，拿到一個樣本矩陣，我們最先要明確的就是一行是一個樣本還是一個維度，心中明確這個整個計算過程就會順流而下，這麼一來就不會迷茫了~ 

第二部分-R軟體協方差cov

R軟體自帶的協方差函式

function (x, y = NULL, use = "everything", method = c("pearson", 
    "kendall", "spearman")) 
{
    na.method <- pmatch(use, c("all.obs", "complete.obs", "pairwise.complete.obs", 
        "everything", "na.or.complete"))
    if (is.na(na.method)) 
        stop("invalid 'use' argument")
    method <- match.arg(method)
    if (is.data.frame(y)) 
        y <- as.matrix(y)
    if (is.data.frame(x)) 
        x <- as.matrix(x)
    if (!is.matrix(x) && is.null(y)) 
        stop("supply both 'x' and 'y' or a matrix-like 'x'")
    stopifnot(is.numeric(x) || is.logical(x), is.atomic(x))
    if (!is.null(y)) 
        stopifnot(is.numeric(y) || is.logical(y), is.atomic(y))
    Rank <- function(u) {
        if (length(u) == 0L) 
            u
        else if (is.matrix(u)) {
            if (nrow(u) > 1L) 
                apply(u, 2L, rank, na.last = "keep")
            else row(u)
        }
        else rank(u, na.last = "keep")
    }
    if (method == "pearson") 
        .Call(C_cov, x, y, na.method, method == "kendall")
    else if (na.method %in% c(2L, 5L)) {
        if (is.null(y)) {
            .Call(C_cov, Rank(na.omit(x)), NULL, na.method, method == 
                "kendall")
        }
        else {
            nas <- attr(na.omit(cbind(x, y)), "na.action")
            dropNA <- function(x, nas) {
                if (length(nas)) {
                  if (is.matrix(x)) 
                    x[-nas, , drop = FALSE]
                  else x[-nas]
                }
                else x
            }
            .Call(C_cov, Rank(dropNA(x, nas)), Rank(dropNA(y, 
                nas)), na.method, method == "kendall")
        }
    }
    else if (na.method != 3L) {
        x <- Rank(x)
        if (!is.null(y)) 
            y <- Rank(y)
        .Call(C_cov, x, y, na.method, method == "kendall")
    }
    else stop("cannot handle 'pairwise.complete.obs'")
}
<bytecode: 0x08ab3808>
<environment: namespace:stats>
>