PCA簡介
??PCA（PrincipalComponents Analysis）即主成分分析，是一個非監(jiān)督的機(jī)器學(xué)習(xí)算法，它是最常用的降維方法之一，通過正交變換將一組可能存在相關(guān)性的變量數(shù)據(jù)轉(zhuǎn)換為一組線性不相關(guān)的變量，轉(zhuǎn)換后的變量被稱為主成分。
??例如下圖所示，樣本有2個特征，現(xiàn)在對該樣本進(jìn)行降維處理。首先考慮直接選擇特征1或者特征2降維，經(jīng)過降維后的樣本由2維降到1維，如圖所示。
??可以看出剔除特征2降維比剔除特征1降維的樣本間的間距更大，即樣本可區(qū)分度更大。那是否還有其他的映射方式，使得映射后樣本的間距更大，事實(shí)上還可以選擇某個軸線，例如下圖所示，樣本映射到該軸線之后，有更大的間距。
??PCA降維的思想就是尋找某個軸線，使得樣本映射到該軸線上后使得樣本區(qū)分度更大。衡量可區(qū)分度的指標(biāo)即方差（事實(shí)上，方差背后的意義就是數(shù)據(jù)的稀疏程度），現(xiàn)在的問題就是如何求得該軸線，使得方差的值最大。
求解思路
??用方差來定義樣本的間距，方差越大表示樣本分布越稀疏，方差越小表示樣本分布越密集，方差的公式如下。
??在求解最大方差前，為了方便計(jì)算，可以先對樣本進(jìn)行demean（去均值）處理，即減去每個特征的均值，這種處理方式不會改變樣本的相對分布（效果就像坐標(biāo)軸進(jìn)行了移動）。去均值后，樣本x每個特征維度上的均值都是0，方差的公式轉(zhuǎn)換成圖示的公式。
??對于2個維度的樣本，我們想要求一個軸的方向 w = (w1, w2) 使得我們所有的樣本，映射到w以后方差最大，有：
??如下圖所示，原始樣本和映射后的樣本有如下的空間幾何關(guān)系。w為單位向量，所以w的模為1。
??樣本x為已知，現(xiàn)在的目標(biāo)轉(zhuǎn)換成如下:
梯度上升發(fā)求解
公式推導(dǎo)
代碼實(shí)現(xiàn)
初始樣本準(zhǔn)備
import numpy as np
import matplotlib.pyplot as plt
X = np.empty((200, 2))
X[:,0] = np.random.uniform(0., 100., size=200)
X[:,1] = 0.75 * X[:,0] + 1.5 + np.random.normal(0, 10., size=200)
plt.scatter(X[:,0], X[:,1])
plt.show()
開始對該樣本進(jìn)行降維處理
def demean(X):
# axis=0 是 矩陣X - X每一列的均值（即axis=0） 矩陣運(yùn)算
return X - np.mean(X, axis=0)
# pca 目標(biāo)函數(shù)推導(dǎo)公式
def f(w, X):
return np.sum((X.dot(w)**2)) / len(X)
# 求目標(biāo)函數(shù)對應(yīng)的梯度
def df_math(w, X):
#即通過數(shù)學(xué)推導(dǎo)得出的求梯度公式
return X.T.dot(X.dot(w)) * 2. / len(X)
def direction(w):
# 由于w不一定是一個單位向量，把它轉(zhuǎn)化成一個單位向量， w/w.模
return w / np.linalg.norm(w)
def gradient_ascent(df, X, initial_w, eta, n_iters = 1e4, epsilon=1e-8):
w = direction(initial_w)
cur_iter = 0
while cur_iter < n_iters:
gradient = df(w, X)
last_w = w
w = w + eta * gradient
w = direction(w) # 注意1：每次求一個單位方向，方便能得到一個合適的值方便計(jì)算
if(abs(f(w, X) - f(last_w, X)) < epsilon):
break
cur_iter += 1
return w
對樣本X進(jìn)行降維
X_mean = demean(X) #去均值處理
initial_w = np.random.random(X.shape[1]) # 注意：不能用0向量開始
eta = 0.001 # 初始化步長
w = gradient_ascent(df_math, X_mean, initial_w, eta) #降維并得出第一個方向，即第一個主成分，使得在該方向上進(jìn)行映射可以有最大的間距（方差）
plt.scatter(X_mean[:,0], X_mean[:,1])
plt.plot([0, w[0]*30], [0, w[1]*30], color='r') # 將w表示的方向進(jìn)行繪制，為了更清楚地展示，這里適當(dāng)擴(kuò)大w的表示范圍（這里擴(kuò)大了30倍），但是w0和w1比例要保持一致。
plt.show()
這里求的是第一主成分，在更高維度的空間，可能還需要繼續(xù)求出其他的主成分，即降更多的維度。
求前n主成分
??對于只有2維特征的樣本X，我們求解出了第一主成分，即w=（w1,w2），為了求解第二主成分，這里需要將樣本X在第一主成分上的分量去除掉，這里使用的方法即空間幾何的向量減法，得到的結(jié)果即下圖中的綠線部分。
??在得到綠線部分的分量后，再對該分量重新求第一主成分，以此類推，可以得到前n個主成分。
求前n主成分代碼實(shí)現(xiàn)
def f(w, X):
return np.sum((X.dot(w)**2)) / len(X)
def df(w, X):
return X.T.dot(X.dot(w)) * 2. / len(X)
def direction(w):
return w / np.linalg.norm(w)
def first_component(X, initial_w, eta, n_iters = 1e4, epsilon=1e-8):
w = direction(initial_w)
cur_iter = 0
while cur_iter < n_iters:
gradient = df(w, X)
last_w = w
w = w + eta * gradient
w = direction(w)
if(abs(f(w, X) - f(last_w, X)) < epsilon):
break
cur_iter += 1
return w
求解第2個主成分，首先求去除掉第一個主成分后樣本的值
X = demean(X)
initial_w = np.random.random(X.shape[1])
eta = 0.01
w = first_component(X, initial_w, eta) #這里求解出第一個主成分
X2 = np.empty(X.shape) # 初始1個矩陣，和X的維度保持一致，用來存儲去除掉第一個主成分分量后的樣本
for i in range(len(X)):
X2[i] = X[i] - X[i].dot(w) * w # 得到每個樣本在去除掉第一個主成分分量后的值
#將該分量繪制出來
plt.scatter(X2[:,0], X2[:,1])
plt.show()
對剩下的分量求第一主成分
w2 = first_component(X2, initial_w, eta)
# 第一主成分和第二主成分結(jié)果相乘接近0，因?yàn)榈谝恢鞒煞趾偷诙鞒煞执怪?，所以矩陣相乘?
w.dot(w2)
代碼封裝
#Author:Sunshine丶天
import numpy as np
class PCA:
def __init__(self, n_components):
"""初始化PCA"""
assert n_components >= 1, "n_components must be valid"
self.n_components = n_components
self.components_ = None
def fit(self, X, eta=0.01, n_iters=1e4):
"""獲得數(shù)據(jù)集X的前n個主成分"""
assert self.n_components <= X.shape[1], "n_components must not be greater than the feature number of X"
def demean(X):
return X - np.mean(X, axis=0)
def f(w, X):
return np.sum((X.dot(w) ** 2)) / len(X)
def df(w, X):
return X.T.dot(X.dot(w)) * 2. / len(X)
def direction(w):
return w / np.linalg.norm(w)
def first_component(X, initial_w, eta=0.01, n_iters=1e4, epsilon=1e-8):
w = direction(initial_w)
cur_iter = 0
while cur_iter < n_iters:
gradient = df(w, X)
last_w = w
w = w + eta * gradient
w = direction(w)
if (abs(f(w, X) - f(last_w, X)) < epsilon):
break
cur_iter += 1
return w
X_pca = demean(X)
self.components_ = np.empty(shape=(self.n_components, X.shape[1]))
for i in range(self.n_components):
initial_w = np.random.random(X_pca.shape[1])
w = first_component(X_pca, initial_w, eta, n_iters)
self.components_[i,:] = w
X_pca = X_pca - X_pca.dot(w).reshape(-1, 1) * w
return self
def transform(self, X):
"""將給定的X，映射到各個主成分分量中"""
assert X.shape[1] == self.components_.shape[1]
return X.dot(self.components_.T)
def inverse_transform(self, X):
"""將給定的X，反向映射回原來的特征空間"""
assert X.shape[1] == self.components_.shape[0]
return X.dot(self.components_)
def __repr__(self):
return "PCA(n_components=%d)" % self.n_components
sklearn中的PCA
#Author:Sunshine丶天
from sklearn.decomposition import PCA
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsClassifier
digits = datasets.load_digits() # 加載波士頓房價數(shù)據(jù)
x = digits.data
y = digits.target
# 切分訓(xùn)練集和測試集
x_train,x_test,y_train,y_test = train_test_split(x,y,random_state = 666)
# 使用knn進(jìn)行預(yù)測。
knn = KNeighborsClassifier()
knn.fit(x_train,y_train)
# 預(yù)測準(zhǔn)確率
print(knn.score(x_test,y_test))
print('未降維之前的維數(shù)：', x.shape[1]) # 結(jié)果為64
pca = PCA(0.95) # 傳入的參數(shù)為0.95，代表保留的95%的信息
pca.fit(x_train)
print(pca.n_components_) # 降維之后結(jié)果為28，由原來的64位降到28緯，保留95%的信息
x_train_reduction = pca.transform(x_train) #訓(xùn)練數(shù)據(jù)集降維
x_test_reduction = pca .transform(x_test) #測試數(shù)據(jù)集也必須降維
knn = KNeighborsClassifier()
knn.fit(x_train_reduction,y_train)
print(knn.score(x_test_reduction,y_test))
|
