EOF與PCA算法的python實作和比較前言EOF&PCA 前提準備

前言

自然正交分解（EOF），主成分分析（PCA）原來是一種方法，或者說是比較相似的方法。因為專業的原因，氣象上用到EOF的時候比較多，而大學期間參加數模又接觸到了主成分分析，可是那是的我從來沒想到這是一種方法，哎，看算法原理沒有看仔細。直到後來參加研究所學生賽的時候，學姐告訴我真相，這原來是一種方法啊。最近氣象統計也主要講這塊，是以我就親自實踐了一下來懲罰我的不認真。

EOF&PCA 前提準備

X是經過減去時間距平處理後的資料，即

X是由m個格點的場形成t個時刻的時間序列

X可以被m個兩兩正交（模為特征值）的時間序列和m個兩兩正交的機關空間序列完全表達。此時

為0

其中

（1）

這裡的V相當于坐标旋轉系數，比如，對于一個二維空間，原來的基坐标是i=[1,0]和j=[0,1].有一個場向量是v=5i+6j, 下一時刻變為v2=7i+3j, 那麼x可以表示為

，V就是EOF裡面的模态，也是PCA中的載荷向量或者通俗的叫法是各個原成分分量的權值，對于v1=[3/sqrt(14),2/sqrt(14)]來說，基向量就變成了一個k=3/sqrt(14)i+2/sqrt(14)j, 他對應的時間系數變為=[27/sqrt(14),27/sqrt(14)]，在PCA中是X的一個主成分，也是EOF的時間系數。

還有一點我覺得比較重要的是時間系數和各個模态分量的性質：

空間場：正交，機關向量（模為1）；

（2）

時間上，正交，模為特征值，（如果除以總特征值的和，代表的是對于原場時間序列解釋的方差貢獻率）。

（3）

求證部分的思路：

1.任意向量的m維投影

得到

1.場的總誤差：

可以看到場的誤差是求一個場各個格點的時間的平均誤差

2.第一特征向量的下的E1的表達式：

第一項為場的原始方差，後一項為模态V1代表誤差，即第一時間序列的方差，算出來就是V1特征向量對應的特征值。可以了解為原來場的方差減去第一模态能解釋的方差。我們當然希望能解釋的越多越好，是以E1越小越好。

3.怎麼求v1

上面的拉格朗日求極值就是把限制條件加入到函數中來去掉限制條件。（不太明白為什麼可以這麼做，這是變成柔性邊界了），但可以看定的是lamda肯定不為0，之後我們可以得到

做多模态隻能是m和t中最小的一個，因為X的秩決定這特征值的數量，也就是lamda的數量。

根據上面求導可以知道，lamda就是X協方差的一個特征值。

為了使E1最小，就取特征值最大的了，特征值也是模态方差貢獻率。

附向量求導的方法：

附拉格朗日乘數法則：

程式實作：

機器學習實戰裡面的程式：

'''
Created on Jun 1, 2011

@author: Peter Harrington
'''
from numpy import *

def loadDataSet(fileName, delim='\t'):
    fr = open(fileName)
    stringArr = [line.strip().split(delim) for line in fr.readlines()]
    datArr = [map(float,line) for line in stringArr]
    return mat(datArr)

def pca(dataMat, topNfeat=):
    meanVals = mean(dataMat, axis=)
    meanRemoved = dataMat - meanVals #remove mean
    covMat = cov(meanRemoved, rowvar=)
    eigVals,eigVects = linalg.eig(mat(covMat))
    # print 'V',eigVects[:,0].T*eigVects[:,1]
    eigValInd = argsort(eigVals)            #sort, sort goes smallest to largest
    eigValInd = eigValInd[:-(topNfeat+):-]  #cut off unwanted dimensions
    redEigVects = eigVects[:,eigValInd]       #reorganize eig vects largest to smallest
    # print redEigVects
    lowDDataMat = meanRemoved * redEigVects#transform data into new dimensions
    # print 'T',lowDDataMat[:,0].T*lowDDataMat[:,1]
    reconMat = (lowDDataMat * redEigVects.T) + meanVals
    return lowDDataMat, reconMat

def replaceNanWithMean(): 
    datMat = loadDataSet('secom.data', ' ')
    numFeat = shape(datMat)[]
    for i in range(numFeat):
        meanVal = mean(datMat[nonzero(~isnan(datMat[:,i].A))[],i]) #values that are not NaN (a number)
        datMat[nonzero(isnan(datMat[:,i].A))[],i] = meanVal  #set NaN values to mean
    return datMat
           

我自己編的程式：

################################
# File Name: PCA.py
# Author: Trichtu
# Email: [email protected]
# Created Time: April 20th 2016
##################################
#coding=utf-8 
import numpy as np 

def PCA(mat,value,mode):
    'rows are time or something and columns are variablea and '
    'mode = 1, value is percent,such as 0.99 0.95 0.9'
    'mode = 0, value is number ,such as 3 4 no more than len(matp[1,:])'
    meanval = np.mean(mat,axis=)
    mat_abnormal = mat-meanval# every column minus average
    covMat=np.cov(mat_abnormal,rowvar=)
    eigval,eigvec = np.linalg.eig(np.mat(covMat))# eigvalue is descending ranked
    attribution=np.cumsum(eigval)/sum(eigval)
    if mode==:
        number=[i for i in range(len(attribution)) if attribution[i]>=value][]+
    elif mode==:
        number=value
    Z=eigvec[:,:number]#each columns reprsent a eig vector
    y=mat_abnormal*Z #the component
    remat=y*Z.T+meanval
    return remat,Z,y,attribution
           

用一個二維資料1000時次的資料做實驗，發現時間系數的模是1e-13，并不是完美的0，模态的模是0。用機器學習實戰書中的代碼和python裡面sklearn包的PCA接口程式作對比如下：

第二張圖下面是sklearn裡面的whiten選項為True的時候情況，不知道是怎麼做的。正常情況下whiten=False時，兩個程式一樣。紅線是主成分和載荷向量（模态）重構之後的情況。