吴恩达机器学习作业-ex7（主成分分析）

news2024/9/24 15:17:35

data1

导入库，读取数据，并进行可视化数据

import numpy as np
import scipy.io as sio
import matplotlib.pyplot as plt

#读取数据
path = "./ex7data1.mat"
data = sio.loadmat(path)
# print(data.keys())
X = data.get("X")
# print(X)

#可视化数据
plt.scatter(X[:,0],X[:,1])
plt.show()

PCA的第一步就是要对数据预处理，其中最重要的就是均值归一化

均值归一化， $x_{j}=\frac{x_{j}-\mu _{j}}{s_{j}}$ ,其中 $\mu _{j}=\frac{1}{m}\sum_{i=1}^{m}x_{j}^{(i)}$ （均值）, $s_{j}$ 代表方差。

#将数据标准化再展示
def standardization(X):
    mean = np.mean(X, axis=0)
    std = np.std(X, axis=0)
    return (X - mean) / std

计算协方差矩阵， $\sum =\frac{1}{m}\sum_{i=1}^{m}(x^{(i)})(x^{(i)})^{T}$

#计算协方差矩阵
def covariance_matrix(X):
    m = len(X)
    cov = (X.T @ X)/m
    return cov

计算协方差矩阵的特征向量， $[U,S,V]=svd(sigma)$ （奇异值分解），这里主要要理解U，其中，U是一个n*n的矩阵，因为要降维至k维，所以我们去除矩阵U的前k列，作为k个方向向量，新的特征向量 $z^{(i)}=U^{T}*x^{(i)}$ ,是一个k*1的向量。

#计算协方差矩阵的特征向量
def eigenvector(sigma):
    U,S,V = np.linalg.svd(sigma)
    return U,S,V

降维： $z^{(i)}=U_{reduce}^{T}*x^{(i)}$

#降维
def get_z(X,U,k):
    U_reduce = U[:,:k]
    z = X @ U_reduce
    return z

压缩重现： $x_{approx}\approx U_{reduce}*z$

#压缩重现
def reproduces(z,U,k):
    U_reduce = U[:, :k]
    X_reproduces = z @ U_reduce.T
    return X_reproduces

调用上述函数

k = 1
X_standardization = standardization(X)
# plt.scatter(X_standardization[:,0],X_standardization[:,1])
# plt.show()
cov = covariance_matrix(X_standardization)
U,S,V = eigenvector(cov)
# print(X_standardization.shape)
# print(U.shape)
z = get_z(X_standardization,U,k)
# print(cov.shape)
# print(U.shape)
# print(z)
X_reproduces = reproduces(z,U,k)

绘制出最开始数据集的散点图和压缩重现后的散点图

fig,ax = plt.subplots()
ax.scatter(X_standardization[:,0],X_standardization[:,1],c="b")
ax.scatter(X_reproduces[:,0],X_reproduces[:,1],c="r")
plt.show()

如果想要展示出，相对应的点，可以加上下面的代码


# for i in range(len(X_standardization)):
#     plt.plot([X_standardization[i,0],X_reproduces[i,0]],[X_standardization[i,1],X_reproduces[i,1]],"k--")

faces

读取数据

#读取数据
path_faces = "ex7faces.mat"
data_faces = sio.loadmat(path_faces)
# print(data_faces.keys())
X_faces = data_faces.get("X")
# print(X.shape)

可视化

这里我定义了两个函数，分别用来选择其中一个样本进行绘制和选择100个样本。

#可视化
def one_image(X_faces):
    pick_one = np.random.randint(0,len(X))
    image = X_faces[pick_one,:]
    image = image.reshape(32,32)
    fig,ax = plt.subplots()
    ax.imshow(image.T,cmap="Greys_r")
    plt.xticks([])
    plt.yticks([])
    plt.show()

# one_image(X_faces)
def more_image(X_faces):
    pick_more = np.random.randint(0,len(X),100)
    images = X_faces[pick_more,:]
    fig,ax = plt.subplots(nrows=10,ncols=10, figsize=(8, 8), sharex=True, sharey=True)
    for i in range(10):
        for j in range(10):
            ax[i,j].imshow(images[10*i+j].reshape(32,32).T,cmap="Greys_r")
    plt.xticks([])
    plt.yticks([])
    plt.show()

# more_image(X_faces)

主成分数量k的选择

PCA所做的就是试图寻找一个投影平面对数据进行投影，使每个点到它们对应的支线上的投影点（投影误差）最小化。

投影的平均均方误差： $\frac{1}{m}\sum_{i=1}^{m}\left \| x^{(i)}-x_{approx}^{(i)} \right \|^{2}$

训练集的方差： $\frac{1}{m}\sum_{i=1}^{m}\left \| x^{(i)} \right \|^{2}$

我们希望找到最小的k使得 $\frac{\frac{1}{m}\sum_{i=1}^{m}\left \| x^{(i)}-x_{approx}^{(i)} \right \|^{2}}{\frac{1}{m}\sum_{i=1}^{m}\left \| x^{(i)} \right \|^{2}}$ （平均均方误差与训练集方差的比例）是一个很小的值。如果我们希望这个比例小于等于1%，也就是说我们保留了99%的方差。

#确定k值
def get_k(X_faces,S):
    for k in range(X_faces.shape[1]):
        error = 1 - np.sum(S[:k])/np.sum(S)
        if error <=0.01:
            return k
        else:
            k = k + 1

调用函数

展示降维前和降维后的图片

X_faces_standardization = standardization(X_faces)
cov_faces = covariance_matrix(X_faces_standardization)
U_faces,S_faces,V_faces = eigenvector(cov_faces)
k_faces = get_k(X_faces,S_faces)
z_faces = get_z(X_faces_standardization,U_faces,k_faces)
X_reproduces = reproduces(z_faces,U_faces,k_faces)
more_image(X_reproduces),more_image(X_faces)
print(k_faces)