学了PCA以后,不用SKlearn现成的包,基于numpy和scipy自己实现了一下PCA算法,原理和SKlearn里PCA是一样的,都是基于SVD方法实现的。可以按照顺序把代码贴进去,自己跑一下试试。
导入需要的包
import time #调用时间,显示算法运行时间
import os
import math
import numpy as np
import scipy as sp
import pandas as pd
# 使np矩阵不显示科学计数
np.set_printoptions(suppress=True)
主体PCA代码
class PCA(object):
"""
Principal component analysis (PCA).
Linear dimensionality reduction using Singular Value Decomposition of the
data to project it to a lower dimensional space. The input data is centered
but not scaled for each feature before applying the SVD.
Parameters
----------
n_component_ratio: select the number of components such that the amount of variance that needs to be
explained is greater than the percentage specified by n_component_ratio
----------
function1-fit_transform: Fit the model by computing full SVD on X.
Parameters
----------
X: x_train
----------
function2-transform: Fit the model with X and apply the dimensionality reduction on y.
Parameters
----------
y : x_test
----------
"""
def __init__(self,n_component_ratio):
self.n_component_ratio=n_component_ratio
def fit_transform(self,x):
n_samples, n_features = x.shape
X_center = x-np.mean(x, axis=0) # Center data
U, s, Vt=sp.linalg.svd(X_center, full_matrices=False)
self.__Vt=Vt
explained_variance = (s ** 2) / (n_samples - 1) # Get variance explained by singular values
total_var = explained_variance.sum()
explained_variance_ratio = explained_variance / total_var #caculate the compressed rate
self.n_components=self.__choose_ratio(explained_r=explained_variance_ratio) #find how many features we should keep on the compressed rate we select
x_compressed = U[:, :self.n_components].dot(np.diag(s[:self.n_components])) #return the features we choose
del X_center,explained_variance,U,s,Vt #del variable and release memory
self.explained_variance_ratio=explained_variance_ratio[:self.n_components] #make explained variance ratio the attributes in pca, so that we can print it out
return x_compressed
def __choose_ratio(self,explained_r):
for i in range(1, len(explained_r)):
if sum(explained_r[:i])/sum(explained_r) >= self.n_component_ratio:
return i
def transform(self,y):
y_centre = y-np.mean(y, axis=0)
y_compressed=y_centre.dot(np.linalg.inv(self.__Vt)[:,:self.n_components]) # compress x_test based on the Vt on x_train
del y_centre #del variable and release memory
gc.collect()
return y_compressed.values
找个数据集测试一下:
from sklearn.datasets import load_digits # 经典手写数字分类数据集
from sklearn.model_selection import train_test_split
digits= load_digits()
digits.target=pd.DataFrame(digits.target)
digits.data=pd.DataFrame(digits.data)
X_train, X_test, y_train, y_test = train_test_split(digits.data, digits.target, test_size=0.4, random_state=1)
数据集原来有64个特征值
在这里使用手写的PCA来降维:
start_time = time.time()
pca = PCA(n_component_ratio=0.95) #这里只支持输入保留的可解释方差率是多少,小数形式,不支持整数
X_train_pca = pca.fit_transform(X_train)
X_test_pca=pca.transform(X_test)
pca_n=X_train_pca.shape[1]
print("keep {} features after pca \nthe shape of X_train after PCA: {} \nthe shape of X_test after PCA: {}".format(pca_n,X_train_pca.shape,X_test_pca.shape))
print("--- %s seconds ---" % (time.time() - start_time))
在保留95%可解释方差的情况下降维,最后保留的特征值有28个
可视化一下:
var=np.cumsum(np.round(pca.explained_variance_ratio, decimals=3)*100) # caculate the cumulative explained variance ratio
plt.plot(list(range(1,len(var)+1)), var)
plt.title("PCA Variance Explained Based on Custom PCA")
plt.ylabel("% Variance Explained")
plt.xlabel("# Of Features")
plt.show()
PCA以后,前两行降维后的数据:
再用sklearn自带的PCA包对比一下:
from sklearn.decomposition import PCA
start_time = time.time()
pca2 = PCA(n_components=0.95,copy=True, whiten=False)
X_train_pca2 = pca2.fit_transform(X_train)
X_test_pca2=pca2.transform(X_test)
pca_n2=X_train_pca2.shape[1]
print("keep {} features after pca \nthe shape of X_train after PCA: {} \nthe shape of X_test after PCA: {}".format(pca_n2,X_train_pca2.shape,X_test_pca2.shape))
print("--- %s seconds ---" % (time.time() - start_time))
可视化一下可解释方差率曲线
var=np.cumsum(np.round(pca2.explained_variance_ratio_, decimals=3)*100) # caculate the cumulative explained variance ratio
plt.plot(list(range(1,len(var)+1)), var)
plt.title("PCA Variance Explained Based on sklearn PCA")
plt.ylabel("% Variance Explained")
plt.xlabel("# Of Features")
plt.show()
对比以后,前面手写的PCA和sklearn里PCA算法降维结果一样。sklearn的PCA稍微快一点。
算法主体部分代码放到GitHub上了,有兴趣欢迎关注:
https://github.com/JuneYaooo/ml-algorithms