http://blog.csdn.net/u011239443/article/details/77363466

协方差矩阵

统计学的基本概念

协方差

上面几个统计量看似已经描述的差不多了，但我们应该注意到，标准差和方差一般是用来描述一维数据的，但现实生活我们常常遇到含有多维数据的数据集，最简单的大家上学时免不了要统计多个学科的考试成绩。面对这样的数据集，我们当然可以按照每一维独立的计算其方差，但是通常我们还想了解更多，协方差就是这样一种用来度量两个随机变量关系的统计量，我们可以仿照方差的定义：

来度量各个维度偏离其均值的程度，标准差可以这么来定义：

从协方差的定义上我们也可以看出一些显而易见的性质，如：

协方差矩阵

举一个简单的三维的例子，假设数据集有{x,y,z}{x,y,z}三个维度，则协方差矩阵为：

求解协方差矩阵的步骤

举个例子：

PCA

算法步骤

1. 形成样本矩阵，样本中心化
2. 计算样本矩阵的协方差矩阵
3. 对协方差矩阵进行特征值分解，选取最大的 p 个特征值对应的特征向量组成投影矩阵
4. 对原始样本矩阵进行投影，得到降维后的新样本矩阵

推导

为什么PCA和协方差扯上关系呢？接下来我们来进行推导：

对于每个点：

会有一个对应的编码向量：

我们希望找到一个编码函数，根据输入返回编码，f(x) = c；我们也希望找到一个解码函数，给定编码重构输入：

为了简化解码器，我们使用矩阵乘法将编码映射回$R^n$ ，即 g(c) = Dc

推导到这里，我们可以看到我们的最有解和协方差矩阵的联系。其实协方差矩阵度量的是维度与维度之间的关系，而非样本与样本之间。协方差矩阵的主对角线上的元素是各个维度上的方差(即能量)，其他元素是两两维度间的协方差(即相关性)。

实现

from numpy import *


    
    
   
    

# 加载数据
def loadDataSet(filename,delim = '\t'):
    fr = open(filename)
    stringArr = [line.strip().split(delim) for line in fr.readlines()]
    dataArr = [map(float,line) for line in stringArr]
    return mat(dataArr)

def pca(dataMat,topN=999999):
    # 形成样本矩阵，样本中心化
    meanVals= mean(dataMat,axis=0)
    meanRemoved = dataMat - meanVals
    # 计算样本矩阵的协方差矩阵
    covMat = cov(meanRemoved,rowvar=0)
    #  对协方差矩阵进行特征值分解，选取最大的 p 个特征值对应的特征向量组成投影矩阵
    eigVals,eigVects =  linalg.eig(mat(covMat))
    eigValInd = argsort(eigVals)
    eigValInd = eigValInd[:-(topN+1):-1]
    redEigVects = eigVects[:,eigValInd]
    # 对原始样本矩阵进行投影，得到降维后的新样本矩阵
    lowDDataMat = meanRemoved * redEigVects
    reconMat = (lowDDataMat * redEigVects.T)+meanVals
    return lowDDataMat,reconMat

import myPca
import matplotlib
import matplotlib.pyplot as plt
dataMat = myPca.loadDataSet('testSet.txt')
lowMat,reconMat = myPca.pca(dataMat,1)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(dataMat[:,0].flatten().A[0],dataMat[:,1].flatten().A[0],marker='^',s=90)
ax.scatter(reconMat[:,0].flatten().A[0],reconMat[:,1].flatten().A[0],marker='o',s=50,c='red')
plt.show()