#?-*-?coding:?utf-8?-*-
“““
Created?-0on?Sun?Feb?28?10:04:26?2016
PCA?source?code
@author:?MM
PCA基本步驟：
對數(shù)據(jù)進行歸一化處理（代碼中并非這么做的，而是直接減去均值）
計算歸一化后的數(shù)據(jù)集的協(xié)方差矩陣
計算協(xié)方差矩陣的特征值和特征向量
保留最重要的k個特征（通常k要小于n），也可以自己制定，也可以選擇一個閾值，然后通過前k個特征值之和減去后面n-k個特征值之和大于這個閾值，則選擇這個k
找出k個特征值對應的特征向量
將m?*?n的數(shù)據(jù)集乘以k個n維的特征向量的特征向量（n?*?k）得到最后降維的數(shù)據(jù)。
“““
import?numpy?as?np
import?gdal
import?os
from?osgeo?import?gdal_array?as?ga
import?matplotlib.pyplot?as?plt

#?計算均值要求輸入數(shù)據(jù)為numpy的矩陣格式，行表示樣本數(shù)，列表示特征
def?meanX（dataX）:
????return?np.mean（dataX?axis=0）??#?axis=0表示按照列來求均值，如果輸入list則axis=1


#?計算方差傳入的是一個numpy的矩陣格式，行表示樣本數(shù)，列表示特征
def?variance（X）:
????m?n?=?np.shape（X）
????mu?=?meanX（X）
????muAll?=?np.tile（mu?（m?1））
????X1?=?X?-?muAll
????variance?=?1.?/?m?*?np.diag（X1.T?*?X1）
????return?variance


#?標準化傳入的是一個numpy的矩陣格式，行表示樣本數(shù)，列表示特征
def?normalize（X）:
????m?n?=?np.shape（X）
????mu?=?meanX（X）
????muAll?=?np.tile（mu?（m?1））
????X1?=?X?-?muAll
????X2?=?np.tile（np.diag（X.T?*?X）?（m?1））
????XNorm?=?X1?/?X2
????return?XNorm


“““
參數(shù)：
	-?XMat：傳入的是一個numpy的矩陣格式，行表示樣本數(shù)，列表示特征????
	-?k：表示取前k個特征值對應的特征向量
返回值：
	-?finalData：參數(shù)一指的是返回的低維矩陣，對應于輸入?yún)?shù)二
	-?reconData：參數(shù)二對應的是移動坐標軸后的矩陣
“““
def?pca（XMat?k）:
????average?=?meanX（XMat）?????##計算均值為進行歸一化做準備
????m?n?=?np.shape（XMat）?????#行，列（7）
????print?（mn）
????data_adjust?=?[]
????avgs?=?np.tile（average?（m?1））???#重復某個數(shù)組。比如tile（An），功能是將數(shù)組A重復n次，構(gòu)成一個新的數(shù)組
????data_adjust?=?XMat?-?avgs????????#?對數(shù)據(jù)進行歸一化處理（代碼中并非這么做的，而是直接減去均值）

????covX?=?np.cov（data_adjust.T）??#?計算協(xié)方差矩陣
????featValue?featVec?=?np.linalg.eig（covX）??#?求解協(xié)方差矩陣的特征值和特征向量
????print?（np.shape（featVec））
????index?=?np.argsort（-featValue）??#?按照featValue進行從大到小排序
????#確定k的主成分的個數(shù)
????sum=0
????max=0
????for?i?in?range（len（index））:
????????sum=sum+featValue[index[i]]
????????if?max/sum<0.99:
????????????max=max+featValue[index[i]]
????????????k=k+1
????finalData?=?[]
????print?（“k?=?“k）
????if?k?>?n:
????????print?（“k?must?lower?than?feature?number“）
????????return
????else:
????????#?注意特征向量是列向量，而numpy的二維矩陣（數(shù)組）a[m][n]中，a[1]表示第1行值
????????selectVec?=?np.matrix（featVec[index[:k]]）??#?找出k個特征值對應的特征向量
????????print?（“特征向量維數(shù)“）
????????print?（np.shape（selectVec））
????????#finalData?=?data_adjust?*?selectVec.T
????????reconData?=?（data_adjust?*?selectVec.T）?+?average[:6]?????#將m?*?n的數(shù)據(jù)集乘以k個n維的特征向量的特征向量（m??*?k）得到最后降維的數(shù)據(jù)。
????????print（?np.shape（reconData））
????????array?=?[]
????????array?=?np.array（array）??#?列表轉(zhuǎn)數(shù)組
????????for?j?in?range（k）:
????????????array?=?np.append（array?reconData[:j]）???????????#將reconData從二維（3列）數(shù)組轉(zhuǎn)化從一維數(shù)組（1維）
????????print?（np.shape（array））
????????array1?=?array.reshape（k?im_height?im_width）?????????????#將?一維數(shù)組??轉(zhuǎn)成??3維矩陣（k高（行），寬（列））
????????out?=?ga.SaveArray（array1?os.path.join（path?“after.img“）?format=“GTiff“?prototype=img）
????????#array2=array.reshape（im_heightim_widthk）
????????#print?（np.shape（array2））
????????#out1=ga.SaveArray（array2os.path.join（path“after22.img“）format=“GTiff“prototyp