importnumpy as npimportpandas as pdimportscipy.io as spioimportmatplotlib.pyplot as pltfrom scipy import misc #图片操作
deffeatureNormalize(X):
n= X.shape[1]
mu= np.mean(X,axis=0)
sigma= np.std(X,axis=0)for i inrange(n):
X[:,i]= (X[:,i] - mu[i]) /sigma[i]returnX,mu,sigmadefplot_data_2d(X,marker):
plt.plot(X[:,0],X[:,1],marker)returnpltdefdrawline(plt,p1,p2,line_type):
plt.plot(np.array([p1[0],p2[0]]),np.array([p1[1],p2[1]]),line_type)defprojectData(X_norm,U,K):
U_reduce=U[:,0:K]returnnp.dot(X_norm,U_reduce)defdisplay_imageData(imgData):
sum=0'''显示100个数(若是一个一个绘制将会非常慢,可以将要画的图片整理好,放到一个矩阵中,显示这个矩阵即可)
- 初始化一个二维数组
- 将每行的数据调整成图像的矩阵,放进二维数组
- 显示即可'''m,n=imgData.shape
width=np.int32(np.round(np.sqrt(n)))
height= np.int32(n/width);
rows_count=np.int32(np.floor(np.sqrt(m)))
cols_count= np.int32(np.ceil(m/rows_count))
pad= 1display_array= -np.ones((pad+rows_count*(height+pad),pad+cols_count*(width+pad)))for i inrange(rows_count):for j inrange(cols_count):
max_val=np.max(np.abs(imgData[sum,:]))
display_array[pad+i*(height+pad):pad+i*(height+pad)+height,pad+j*(width+pad):pad+j*(width+pad)+width] = imgData[sum,:].reshape(height,width,order="F")/max_val #order=F指定以列优先,在matlab中是这样的,python中需要指定,默认以行
sum += 1plt.figure(figsize=(10,8))
plt.imshow(display_array,cmap='gray') #显示灰度图像
plt.axis('off')
plt.show()defrecoverData(Z,U,K):
X_rec=np.zeros((Z.shape[0],U.shape[0]))
U_recude=U[:,0:K]
X_rec= np.dot(Z,np.transpose(U_recude)) #还原数据(近似)
returnX_recdefPCA_Face():print("加载图像数据")
data= spio.loadmat("data_faces.mat")
X= data['X']
m,n=X.shape#print(pd.DataFrame(X))
display_imageData(X[0:100,:])
X_norm,mu,sigma=featureNormalize(X)
Sigma= np.dot(np.transpose(X_norm),X_norm)/m #求Sigma
U,S,V = np.linalg.svd(Sigma) #奇异值分解
display_imageData(np.transpose(U[:,0:100])) #显示U的数据
K= 100 #降维100维(原先是32*32=1024维的)
Z =projectData(X_norm, U, K)print (u'投影之后Z向量的大小:%d %d' %Z.shape)print (u'显示降维之后的数据......')
X_rec= recoverData(Z, U, K) #恢复数据
display_imageData(X_rec[0:100,:])#数据降维
defPCA_2D():
data= spio.loadmat("pcadata.mat")
X= data['X']#print(X)
m,n =X.shape#画图
plt = plot_data_2d(X,'o')
plt.show()
X_copy=X.copy()
X_norm,mu,sigma=featureNormalize(X_copy)#print(X_norm,mu,sigma)
#plot_data_2d(X_norm,'o') #画归一化数据
#plt.show()
Sigma= np.dot(np.transpose(X_norm),X_norm) /m
U,S,V=np.linalg.svd(Sigma)#print(U,S,V)
plt= plot_data_2d(X,'o')
drawline(plt, mu, mu+S[0]*(U[:,0]), 'r--')
plt.show()
K= 1 #设置降维维数, 这里是降了 1 维
Z = projectData(X_norm,U,K) #投影
#print(Z) #能观察到数据以及是1维的了, 说明成功将数据降到了一维
#恢复数据
X_rec = recoverData(Z,U,K) #恢复 此处的 U 等于 V的转置
plt = plot_data_2d(X_norm,'bo')
plot_data_2d(X_rec,'ro')for i inrange(X_norm.shape[0]):
drawline(plt, X_norm[i,:], X_rec[i,:],'--k')
plt.axis('square')
plt.show()if __name__ == "__main__":
PCA_2D()
PCA_Face()