#PCA算法的底层原理实现
import numpy as np
import matplotlib.pyplot as plt
x=np.empty((100,2))
np.random.seed(666)
#噪声数据验证
x[:,0]=np.random.uniform(0.0,100.0,size=100)
x[:,1]=0.75*x[:,0]+3.0*np.random.normal(0,3,size=100)
#无噪声数据验证
#x[:,0]=np.random.uniform(0.0,100.0,size=100)
#x[:,1]=0.75*x[:,0]+3.0
plt.figure()
plt.scatter(x[:,0],x[:,1])
plt.show()
#demean操作,每一个特征进行均值归0
def demean(x):
return x-np.mean(x,axis=0)
x_demean=demean(x)
plt.scatter(x_demean[:,0],x_demean[:,1])
plt.show()
#梯度上升法的函数定义
def f(w,x):
return np.sum((x.dot(w))**2)/len(x)
def df_math(w,x):
return x.T.dot(x.dot(w))*2/len(x)
def df_debug(w,x,epsilon=0.00001):
res=np.empty(len(w))
for i in range(len(w)):
w1=w.copy()
w1[i]=w1[i]+epsilon
w2= w.copy()
w2[i] =w2[i]-epsilon
res[i]=(f(w1,x)-f(w2,x))/(2*epsilon)
return res
#向量单位化函数
def derection(w):
##求取向量的模长np.linalg.norm(w)
return w/np.linalg.norm(w)
def gradient_ascent1(df,x,eta,w_initial,erro=1e-8, n=1e6):
w=w_initial
w=derection(w)
i=0
while i<n:
gradient =df(w,x)
last_w = w
w = w + gradient * eta
w = derection(w) #注意1：每次都需要将w规定为单位向量
if (abs(f(w,x) - f(last_w,x))) < erro:
break
i+=1
return w
w0=np.random.random(x.shape[1]) #注意2：不能从0向量初始化初始点，因为0向量是极小值点，导数为0
print(w0)
eta=0.001 #注意3：不能将数据进行标准化，即不可以使用standardscaler进行数据标准化，因为标准化之后数据的方差不变了，而不是求取最大值
w1=gradient_ascent1(df_debug,x_demean,eta,w0)
w2=gradient_ascent1(df_math,x_demean,eta,w0)
print(w1,w2)
plt.scatter(x[:,0],x[:,1])
plt.plot([0,w1[0]*30],[0,w1[1]*30],color="r")
print(w1[1]/w1[0])
plt.show()
#求取数据的前n个的主成分，循环往复即可
x2=np.empty(x.shape)
for i in range(len(x)):
x2[i]=x_demean[i]-x_demean[i].dot(w1)*w1\
#计算向量化，结果是等效的的
x2=x-x.dot(w1).reshape(-1,1)*w1
plt.scatter(x2[:,0],x2[:,1],color="g")
plt.show()
w2=gradient_ascent1(df_debug,x2,eta,w0)
print(w2)
print(w1.dot(w2))
#求取n维数据的前n个主成分的封装函数
def first_n_compnent(n,x,eta=0.001,erro=1e-8, m=1e6):
x_pca=x.copy()
x_pca=demean(x_pca)
res=[]
for i in range(n):
w0=np.random.random(x_pca.shape[1])
w=gradient_ascent1(df_math,x_pca,eta,w0)
res.append(w)
x_pca=x_pca-x_pca.dot(w).reshape(-1,1)*w
return res
print(first_n_compnent(2,x))
#sklearn中调用PCA算法，基本的使用方法
#以二维数据为例
from sklearn.decomposition import PCA
pca=PCA(n_components=1) #求取第一主成分
#噪声数据验证
x[:,0]=np.random.uniform(0.0,100.0,size=100)
x[:,1]=0.75*x[:,0]+3.0*np.random.normal(0,3,size=100)
pca.fit(x)
w1=pca.components_ #求得输出前n的主成分向量方向
print(w1)
x_reduction=pca.transform(x) #进行数据的降维
x_restore=pca.inverse_transform(x_reduction)#进行数据的恢复，不过恢复之后会丢失部分信息，不完全相同
print(x_reduction.shape)
print(x.shape)
plt.scatter(x[:,0],x[:,1])
plt.scatter(x_restore[:,0],x_restore[:,1])
plt.show()
#以实际的手写字体数据为例进行数据的降维使用
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_digits
d=load_digits() #导入手写字体进行实际的使用实例
x=d.data
y=d.target
from sklearn.model_selection import train_test_split
x_train,x_test,y_train,y_test=train_test_split(x,y,random_state=666)
print(x_train.shape) #(1347,64)数据形状
from sklearn.neighbors import KNeighborsClassifier
knn=KNeighborsClassifier()
knn.fit(x_train,y_train)
print(knn.score(x_test,y_test))
from sklearn.decomposition import PCA
#指定降维的维数
pca=PCA(n_components=15)
pca.fit(x_train)
x_train_reduction=pca.transform(x_train)
x_test_reduction=pca.transform(x_test)
knn=KNeighborsClassifier()
knn.fit(x_train_reduction,y_train)
print(knn.score(x_test_reduction,y_test))
print(pca.explained_variance_ratio_) #输出各个成分轴所保持原来数据方差的大小的所占比例
#画图显示出来各个所有特征轴的主成分所占比例，然后输出累积的反映方差所占比例
pca=PCA(n_components=x_train.shape[1])
pca.fit(x_train)
print(pca.explained_variance_ratio_)
x1=[i for i in range(x_train.shape[1])]
y1=[np.sum(pca.explained_variance_ratio_[:i+1]) for i in range(x_train.shape[1])]
print(x1)
print(y1)
plt.plot([i for i in range(x_train.shape[1])],[np.sum(pca.explained_variance_ratio_[:i+1]) for i in range(x_train.shape[1])])
plt.show()
#直接使用pca算法的所占比例参数进行规定，从而求出此时满足的降低维度的信息
pca=PCA(0.95)
pca.fit(x_train)
print(pca.n_components_)#输出方差保持在0.95时需要降低到的数据维度
x_train_reduction=pca.transform(x_train)
x_test_reduction=pca.transform(x_test)
knn=KNeighborsClassifier()
knn.fit(x_train_reduction,y_train)
print(knn.score(x_test_reduction,y_test))
#降维之后，可以使用二维数据显示手写字体的区分度，可以看出二维数据也具有一定的区分度
d=load_digits() #导入手写字体进行实际的使用实例
x=d.data
y=d.target
pca1=PCA(n_components=2)
pca1.fit(x)
x_reduction=pca1.transform(x)
for i in range(10):
plt.scatter(x_reduction[y==i,0],x_reduction[y==i,1],alpha=0.8)
plt.show()
#MNIST数据集的训练与实践(7000,784)
#可以看出降维只有数据训练速度加快，准确度变化不大，也有可能提高准确度
'''import numpy as np
from sklearn.datasets import fetch_mldata
mnist=fetch_mldata("MNIST original")
print(mnist)
x=mnist["data"]
y=mnist["target"]
print(x.shape)
print(y.shape)
x_train=np.array(x[:60000],dtype=float)
y_train=np.array(y[:60000],dtype=float)
x_test=np.array(x[60000:],dtype=float)
y_test=np.array(y[60000:],dtype=float)
print(x_train.shape)
print(y_train.shape)
knn=KNeighborsClassifier()
knn.fit(x_train,y_train)
print(knn.score(x_test,y_test))
pca=PCA(0.9)
pca.fit(x_train)
x_train=pca.transform(x_train)
x_test=pca.transform(x_test)
print(x_train.shape)
knn=KNeighborsClassifier()
knn.fit(x_train,y_train)
print(knn.score(x_test,y_test))
'''
#降噪功能演示
from sklearn import datasets
digit=load_digits()
x=digit.data
y=digit.target
noisy_digit=x+np.random.normal(0,4,size=x.shape) #增加一定的噪声
example_digit=noisy_digit[y==0,:][:10]
for num in range(1,10):
x_num=noisy_digit[y==num,:][:10]
example_digit=np.vstack([example_digit,x_num])
print(example_digit.shape)
#输出手写字体数据函数
def plot_digits(data):
fig,axes=plt.subplots(10,10,figsize=(10,10),subplot_kw={"xticks":[],"yticks":[]},
gridspec_kw=dict(hspace=0.1,wspace=0.1))
for i ,ax in enumerate(axes.flat):
ax.imshow(data[i].reshape(8,8),
cmap="binary",interpolation="nearest",
clim=(0,16))
plt.show()
plot_digits(example_digit)
pca=PCA(0.5)
pca.fit(noisy_digit)
print(pca.components_)
f=pca.transform(example_digit) #降维除燥
f1=pca.inverse_transform(f) #恢复到原来高维数据
plot_digits(f1)
#人脸识别与特征脸
from sklearn.datasets import fetch_lfw_people
faces=fetch_lfw_people()
print(faces.keys())
print(faces.data.shape)
print(faces.images.shape)
index=np.random.permutation(len(faces.data))
x=faces.data[index]
example=x[:36,:]
#输出人脸数据函数
def plot_faces(data):
fig,axes=plt.subplots(6,6,figsize=(10,10),subplot_kw={"xticks":[],"yticks":[]},
gridspec_kw=dict(hspace=0.1,wspace=0.1))
for i ,ax in enumerate(axes.flat):
ax.imshow(data[i].reshape(62,47),cmap="bone")
plt.show()
plot_faces(example)
print(faces.target_names)
from sklearn.decomposition import PCA
pca=PCA(svd_solver="randomized") #定义一种新的pca算法模型
pca.fit(x)
print(pca.components_.shape)#(2914,2914)
#输出所有特征轴对应的主成分分的特征向量w(nxn方阵)
plot_faces(pca.components_[:36,:])
#输出图片数目至少为60张的样本数据
faces1=fetch_lfw_people(min_faces_per_person=60)
print(faces1.data.shape)
print(len(faces1.target_names))
