2024.9.20 Python模式识别新国大EE5907,PCA主成分分析，LDA线性判别分析，GMM聚类分类，SVM支持向量机

1.PCA 主成分分析用于特征提取、可视化和分类

根据要求，我在第一个代码框中完成了从指定路径提取图像数据，将其转换为灰度图像并将其展平。在这里，我将数字 88 设置为我的照片的标签，然后将所有 10 张照片传入代码。然后我定义了 PCA 函数，计算居中数据，计算协方差矩阵，并计算协方差矩阵的特征值和特征向量。然后对特征向量进行排序，并保留最大的 n 个特征向量。
接下来，我将图像的维数降低到 2D 和 3D。图像在二维和三维空间中均有显示。在图中我使用红点来显示我的图片。

import os
import numpy as np
from PIL import Image
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
from random import sample

dataset_path = 'PIE'
my_photo_label = '88'
image_size = (32, 32)  

def load_images(path, label, num_images):
    images = []
    label_path = os.path.join(path, label)
    available_images = os.listdir(label_path)
    selected_images = sample(available_images, num_images)
    for image_name in selected_images:
        image_path = os.path.join(label_path, image_name)
        with Image.open(image_path) as img:
            img = img.resize(image_size).convert('L') 
            images.append(np.array(img).flatten()) 
    return images

def pca(X, num_components):
    X_meaned = X - np.mean(X , axis = 0)
    cov_mat = np.cov(X_meaned , rowvar = False)
    eigen_values , eigen_vectors = np.linalg.eigh(cov_mat)
    
    sorted_index = np.argsort(eigen_values)[::-1]
    sorted_eigenvalue = eigen_values[sorted_index]
    sorted_eigenvectors = eigen_vectors[:,sorted_index]
    eigenvector_subset = sorted_eigenvectors[:,0:num_components]  
    X_reduced = np.dot(eigenvector_subset.transpose() , X_meaned.transpose()).transpose()
    
    return X_reduced

data = []
for i in range(1, 26):
    data.extend(load_images(dataset_path, str(i), 490 // 25))
data.extend(load_images(dataset_path, '88', 10))
data = np.array(data)

data_2d = pca(data, 2)
data_3d = pca(data, 3)

pca_3 = PCA(n_components=3)
data_3d = pca_3.fit_transform(data)

plt.figure(figsize=(8, 6))
plt.scatter(data_2d[:-10, 0], data_2d[:-10, 1], alpha=0.5)
plt.scatter(data_2d[-10:, 0], data_2d[-10:, 1], color='red')  
plt.title('PCA to 2D')
plt.xlabel('Component 1')
plt.ylabel('Component 2')
plt.show()

from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(data_3d[:-10, 0], data_3d[:-10, 1], data_3d[:-10, 2], alpha=0.5)
ax.scatter(data_3d[-10:, 0], data_3d[-10:, 1], data_3d[-10:, 2], color='red')  
ax.set_title('PCA to 3D')
ax.set_xlabel('Component 1')
ax.set_ylabel('Component 2')
ax.set_zlabel('Component 3')
plt.show()

1.load_images函数，输入一个文件夹里的指定数量的图片进去，然后每张图像被打开、调整大小为 32x32，并转换为灰度图像（convert(‘L’)）将图像数据展平（1D数组）并存储在 images 列表中。返回展平的图像列表，每张图像以 1D 数组形式表示。label只是用来找到要找的图，在这里并没有存储
2.pca函数
输入参数:X 是数据矩阵，每一行是一个样本。num_components 是要保留的主成分数量。
中心化数据，减去每列（每个特征）的均值。
计算协方差矩阵，并对协方差矩阵求特征值和特征向量。
将特征值按降序排列，并根据指定的 num_components 选择对应的特征向量子集。
将数据投影到这些选定的特征向量空间中，获得降维后的数据。
输出: 降维后的数据矩阵。
3.加载数据并执行PCA，降维后输出图片。

这个时候的数据仅是铺开的情况，并没有辨别能力，所以还需要用knn来去辨别

from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

data_40d = pca(data, 40)
data_80d = pca(data, 80)
data_200d = pca(data, 200)

print("Data dimensionality reduced to 40 dimensions:", data_40d.shape)
print("Data dimensionality reduced to 80 dimensions:", data_80d.shape)
print("Data dimensionality reduced to 200 dimensions:", data_200d.shape)

import numpy as np
from collections import Counter

class KNearestNeighbors:
    def __init__(self, k=3):
        self.k = k

    def fit(self, X, y):
        self.X_train = X
        self.y_train = y

    def predict(self, X):
        y_pred = [self._predict(x) for x in X]
        return np.array(y_pred)

    def _predict(self, x):
        distances = [np.sqrt(np.sum((x_train - x) ** 2)) for x_train in self.X_train]
        k_indices = np.argsort(distances)[:self.k]
        k_nearest_labels = [self.y_train[i] for i in k_indices]
        most_common = Counter(k_nearest_labels).most_common(1)
        return most_common[0][0]
def train_and_evaluate_knn(X_train, X_test, y_train, y_test, n_neighbors=3):
    knn = KNearestNeighbors(k=3)
    knn.fit(X_train, y_train)
    y_pred = knn.predict(X_test)
    accuracy = accuracy_score(y_test, y_pred)
    return accuracy

def split_data(data, label_counts, test_size=0.3, random_state=42):
    X_train, X_test, y_train, y_test = [], [], [], []
    for label, count in label_counts.items():
        start_index = sum(label_counts[l] for l in label_counts if int(l) < int(label))
        end_index = start_index + count
        X_label = data[start_index:end_index]
        y_label = [label] * count
        X_label_train, X_label_test, y_label_train, y_label_test = train_test_split(
            X_label, y_label, test_size=test_size, random_state=random_state
        )
        X_train.extend(X_label_train)
        X_test.extend(X_label_test)
        y_train.extend(y_label_train)
        y_test.extend(y_label_test)
    return np.array(X_train), np.array(X_test), np.array(y_train), np.array(y_test)

label_counts = {str(i): 490 // 25 for i in range(1, 26)}
label_counts['88'] = 10 

X_train_40d, X_test_40d, y_train_40d, y_test_40d = split_data(data_40d, label_counts)
X_train_80d, X_test_80d, y_train_80d, y_test_80d = split_data(data_80d, label_counts)
X_train_200d, X_test_200d, y_train_200d, y_test_200d = split_data(data_200d, label_counts)

accuracy_40d = train_and_evaluate_knn(X_train_40d, X_test_40d, y_train_40d, y_test_40d)
accuracy_80d = train_and_evaluate_knn(X_train_80d, X_test_80d, y_train_80d, y_test_80d)
accuracy_200d = train_and_evaluate_knn(X_train_200d, X_test_200d, y_train_200d, y_test_200d)

print("Classification Accuracy with 40 dimensions:", accuracy_40d)
print("Classification Accuracy with 80 dimensions:", accuracy_80d)
print("Classification Accuracy with 200 dimensions:", accuracy_200d)

这个类定义了一个KNN分类器，它通过计算测试样本和训练样本之间的距离，来判断测试样本的类别。KNN的关键思想是基于训练集中与测试样本最相似的 k 个样本来预测测试样本的类别。
首先是class kNN类，他定义了前K个最小距离，
这里的fit其实是一个存取函数，这个函数把训练集直接存起来，也没有进行什么训练，因为之后的判断中，只需要过来找他计算就好了，predict函数是接受train的数据的，然后调用了类的内部函数进行具体的计算，_predict函数的内部是，我先计算要预测的x和所有的x的欧氏距离，然后把距离存下来，然后根据距离排序，取前k个最小值，然后再取这k个数据的原始标签，数这出现的标签最多的是谁那这个图片就应该是谁。
接下来的train and evaluation其实就是在做这个事情，split_data是用来处理原始数据集的，剩下的就没什么好说的了

2.numpy和python列表相比

numpy数组更像是一个动态矩阵，他强大，而且耗费资源少，所以很实用，具体的优点如下：
1.数据类型的一致性，能动态处理所有的数据类型
2.性能好效率高，占用资源少
3.高级运算，内建函数：数组运算更快更方便，就是矩阵的运算。
4.广播机制：np.array: numpy 提供广播机制，这意味着可以对形状不同的数组进行操作，numpy 会自动扩展较小的数组以匹配较大的数组。比如他可以直接给数组+1，那么数组中所有的元素都会+1

3.LDA线性判别分析

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

label_counts = {str(i): 490 // 25 for i in range(1, 26)}
label_counts['88'] = 10 

labels = []
for label, count in label_counts.items():
    labels.extend([label] * count)

labels = np.array(labels)

class LinearDiscriminantAnalysis:
    def __init__(self, n_components, reg_param=0.01):
        self.n_components = n_components
        self.reg_param = reg_param
        self.means_ = None
        self.scalings_ = None

    def fit(self, X, y):
        self.mean_ = np.mean(X, axis=0)
        X = (X - self.mean_)

        class_labels = np.unique(y)
        mean_vectors = [np.mean(X[y == cl], axis=0) for cl in class_labels]

        S_W = sum([(X[y == cl] - mv).T.dot(X[y == cl] - mv) for cl, mv in zip(class_labels, mean_vectors)])
        S_W += np.eye(S_W.shape[0]) * self.reg_param

        overall_mean = np.mean(X, axis=0)
        S_B = sum([len(X[y == cl]) * (mv - overall_mean).reshape(X.shape[1], 1).dot((mv - overall_mean).reshape(1, X.shape[1])) for cl, mv in zip(class_labels, mean_vectors)])

        A = np.linalg.inv(S_W).dot(S_B)
        U, _, _ = np.linalg.svd(A)

        self.scalings_ = U[:, :self.n_components]

    def transform(self, X):
        X = X - self.mean_
        return X.dot(self.scalings_)

    def fit_transform(self, X, y):
        self.fit(X, y)
        return self.transform(X)

def apply_lda(X, y, n_components):
    lda = LDA(n_components=n_components)
    return lda.fit_transform(X, y)

data_2d_lda = apply_lda(data, labels, 2)
data_3d_lda = apply_lda(data, labels, 3)
data_9d_lda = apply_lda(data, labels, 9)

label_to_int = {str(i): i for i in range(1, 26)}
label_to_int['88'] = 88  
int_labels = np.array([label_to_int[label] for label in labels])

#2D
plt.figure(figsize=(8, 6))
plt.scatter(data_2d_lda[:, 0], data_2d_lda[:, 1], c=int_labels, cmap='rainbow', alpha=0.5)
plt.title('LDA: Data projected onto 2 dimensions')
plt.xlabel('LD1')
plt.ylabel('LD2')
plt.show()

#3D
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(data_3d_lda[:, 0], data_3d_lda[:, 1], data_3d_lda[:, 2], c=int_labels, cmap='rainbow', alpha=0.5)
ax.set_title('LDA: Data projected onto 3 dimensions')
ax.set_xlabel('LD1')
ax.set_ylabel('LD2')
ax.set_zlabel('LD3')
plt.show()

在LDA中，SVD用于分解矩阵以提取出最重要的方向。特别是在LDA求解过程中，我们计算了矩阵
𝐴，这个矩阵是类内散度矩阵的逆乘上类间散度矩阵。SVD的目的是分解这个矩阵，提取其特征值和特征向量，从而找到数据在不同类别之间分离最好的方向。
LDA的核心思想是找到一组线性投影，将数据从高维空间投影到低维空间，同时保证在低维空间中，不同类别的样本能够很好地分离开。具体步骤如下：
1.类内散度矩阵 SW：计算每一类样本的类内差异，衡量同一类别样本之间的散布情况。
2.类间散度矩阵 SB：计算不同类别之间的差异，衡量类别之间的中心点的差异。
3.优化目标：LDA的目标是最大化类间散度和类内散度的比率，投影后类间差异越大，类内差异越小，说明投影效果越好。
4.求解矩阵 ：这个矩阵的特征向量表示最佳的投影方向。
5.奇异值分解：通过SVD分解来选择前几个最大的特征值方向，即这些方向是最能区分类别的。

在代码中只要左奇异向量的前n个向量，就是我们要的。然后和PCA一样，可以用knn来进行计算。

4.GMM

GMM是无监督学习算法，它通过聚类来发现数据中的潜在结构，不需要事先知道标签。KNN等算法是有监督学习算法，它们通过已知的标签进行分类。因此，GMM不会直接利用标签，而是通过数据的特征来生成分类。

import numpy as np
from scipy.stats import multivariate_normal
class GMM:
    def __init__(self, n_components=3, tol=1e-4, max_iter=100):
        self.n_components = n_components
        self.tol = tol
        self.max_iter = max_iter

    def fit(self, X):
        n_samples, n_features = X.shape
        self.weights_ = np.full(self.n_components, 1 / self.n_components)
        self.means_ = X[np.random.choice(n_samples, self.n_components, replace=False)]
        self.covariances_ = [np.cov(X.T) for _ in range(self.n_components)]
        
        log_likelihood = 0
        self.converged_ = False
        self.log_likelihoods_ = []
        
        for _ in range(self.max_iter):
            responsibilities = self._e_step(X)
            self._m_step(X, responsibilities)
            
            new_log_likelihood = np.sum(np.log(np.dot(responsibilities, self.weights_)))
            self.log_likelihoods_.append(new_log_likelihood)
            
            if abs(new_log_likelihood - log_likelihood) <= self.tol:
                self.converged_ = True
                break
                
            log_likelihood = new_log_likelihood
            
    def _e_step(self, X):
        likelihood = np.zeros((X.shape[0], self.n_components))
        for i in range(self.n_components):
            likelihood[:, i] = self.weights_[i] * multivariate_normal.pdf(X, self.means_[i], self.covariances_[i])
        responsibilities = likelihood / likelihood.sum(axis=1, keepdims=True)
        return responsibilities
    
    def _m_step(self, X, responsibilities):
        n_samples = X.shape[0]
        for i in range(self.n_components):
            weight = responsibilities[:, i].sum()
            mean = np.dot(responsibilities[:, i], X) / weight
            covariance = (np.dot((responsibilities[:, i] * (X - mean).T), (X - mean)) / weight) + self.tol * np.identity(X.shape[1])
            
            self.weights_[i] = weight / n_samples
            self.means_[i] = mean
            self.covariances_[i] = covariance
            
    def predict_proba(self, X):
        likelihood = np.zeros((X.shape[0], self.n_components))
        for i in range(self.n_components):
            likelihood[:, i] = self.weights_[i] * multivariate_normal.pdf(X, self.means_[i], self.covariances_[i])
        return likelihood / likelihood.sum(axis=1, keepdims=True)

目标是期望最大化，E是期望，M是最大化，通过EM算法在最大迭代次数或达到收敛条件（对数似然变化小于 tol）之间进行迭代：
E步: 计算每个数据点属于每个高斯分布的责任度（即属于某个簇的概率）。
M步: 根据责任度重新估计每个簇的参数（权重、均值和协方差）。
具体如下：
1.初始定义，定义聚类类别数，收敛域值和最大化次数
2.E步中，计算责任度，即每个数据点属于每个簇的概率，用似然值除以所有簇的似然值之和，看看谁最大
3.在 M 步中，算法根据责任度更新模型的参数，具体分为：
权重: 每个簇的权重 𝜋𝑘由该簇的责任度之和决定：
均值: 每个簇的均值 𝜇𝑘是该簇责任度加权后的均值
协方差: 每个簇的协方差矩阵也根据责任度进行更新
这样在最后的时候就可以收敛了，就不会再大的变化了。

问题：
优势：
柔性较强：与K-Means不同，GMM允许每个簇具有不同的形状（由协方差矩阵决定），因此能更好地处理复杂的数据分布。
软分类：GMM为每个数据点分配概率，而不是硬性分配到某个簇，适合处理一些数据边界不清晰的情况。
无监督学习：GMM可以在没有标签的情况下发现数据的潜在结构。
劣势：
对初始参数敏感：GMM的结果可能依赖于初始参数的选择，特别是初始均值的选择。
计算复杂度较高：每次迭代中都需要计算高斯分布的概率，尤其是在处理高维数据时，计算开销较大。
需要假设数据来自高斯分布：如果数据分布与高斯分布假设差距较大，GMM的效果可能不理想。

5.SVM支持向量机

VM 的核心思想是寻找一个能够最大化分类间隔（Margin）的超平面。分类间隔是指超平面到最近的训练样本（支持向量）的距离，SVM 尽可能选择让这个间隔最大的超平面来划分不同类别的数据点。

1.线性可分的情况：在二维空间里，超平面就是一条直线。在高维空间，超平面是一个 d−1 维的平面，用来将不同类别的数据分开。
支持向量：指那些位于分类间隔边界上的样本点，这些点对超平面的最终位置有决定性影响。
SVM 希望找到如下形式的超平面：
𝑤𝑇𝑥+b=0 其中，𝑤是权重向量，b 是偏置，𝑥是输入向量。目标是让不同类别的样本点尽可能远离这个超平面。

2.软间隔与正则化参数 𝐶在实际问题中，数据可能并不是线性可分的，这时就需要允许一定的误分类。这就引入了软间隔和正则化参数
𝐶。
软间隔：允许部分数据点处于超平面的错误一侧，即允许一定量的误分类。C 参数控制误分类的容忍度。具体来说：
大 C：对误分类的容忍度低，会导致模型更加严格地拟合训练数据，可能导致过拟合。
小 C：对误分类的容忍度高，允许更多的误分类，以换取更平滑的分类边界，可能导致欠拟合。

from sklearn.svm import SVC

def train_and_evaluate_svm(X_train, X_test, y_train, y_test, C):
    svm = SVC(C=C, kernel='linear', random_state=42)
    svm.fit(X_train, y_train)
    y_pred = svm.predict(X_test)
    accuracy = accuracy_score(y_test, y_pred)
    return accuracy

C_values = [0.01, 0.1, 1]
X_train_raw, X_test_raw, y_train_raw, y_test_raw = split_data(data, label_counts)

for C in C_values:
    accuracy_raw = train_and_evaluate_svm(X_train_raw, X_test_raw, y_train_raw, y_test_raw, C)
    accuracy_80d = train_and_evaluate_svm(X_train_80d, X_test_80d, y_train_80d, y_test_80d, C)
    accuracy_200d = train_and_evaluate_svm(X_train_200d, X_test_200d, y_train_200d, y_test_200d, C)
    
    print(f"Classification Accuracy with raw vector and C={C}: {accuracy_raw}")
    print(f"Classification Accuracy with 80 dimensions and C={C}: {accuracy_80d}")
    print(f"Classification Accuracy with 200 dimensions and C={C}: {accuracy_200d}")