1. 从公式到代码：为什么你需要亲手推导算法

很多朋友刚学机器学习的时候，总想着赶紧调包跑通代码，看到模型跑出结果就觉得大功告成了。我以前也这样，直到有一次面试，面试官让我在白板上手推逻辑回归的梯度下降公式，我当场卡壳。那次经历让我明白，只会调sklearn的fit和predict，就像只会开车却不懂发动机原理，一旦车子抛锚，你就束手无策了。

理解算法的数学推导，绝不是为了炫耀数学能力。它的实际好处太多了。第一，能帮你真正理解参数的意义。比如线性回归里的正则化项，你看公式就知道L1正则（Lasso）会让一些系数直接变成0，从而实现特征选择；而L2正则（Ridge）只是把系数压缩得更小。光看代码，你可能只知道加了个penalty参数，但为什么这么选，心里是没底的。

第二，能让你高效地调试模型。模型效果不好，是过拟合还是欠拟合？损失函数不收敛，是学习率太大还是数据没归一化？如果你清楚算法每一步在数学上是怎么更新的，你就能像老中医一样，一眼看出问题所在，而不是盲目地乱调参数。第三，能让你有能力改造和创新。很多顶会的论文，其实就是对经典算法损失函数或优化过程的一个小改进。你不懂底层数学，连人家的论文都看不懂，更别提自己创新了。

所以，这篇文章我会带你用Python，把18个核心机器学习算法的数学公式“翻译”成代码。我们不只满足于调用sklearn，更要搞清楚sklearn在背后到底为我们做了什么。我会分享很多我踩过的坑，比如矩阵维度对不齐、梯度爆炸、公式推导符号搞反了等等。相信我，跟着走完这一趟，你对机器学习的理解会深刻得多。

2. 线性模型家族：从回归到正则化

线性模型是机器学习的基石，看似简单，却蕴含着最重要的思想。我们从最基础的开始，一点点增加复杂度。

2.1 一元与多元线性回归：最小二乘法的本质

一元线性回归的公式y = ax + b，初中生都学过。但它的核心思想——最小二乘法，才是关键。我们的目标是找到一条直线，让所有数据点到这条直线的垂直距离的平方和最小。这个距离就是误差，平方是为了消除正负号影响，并且对大的误差给予更大的惩罚。

用数学表达，损失函数 J(a, b) 就是： J(a, b) = Σ(y_i - (a*x_i + b))² 我们要找到使 J(a, b) 最小的 a 和 b。怎么找？在大学微积分里，我们对 a 和 b 分别求偏导数，并令其等于0，解方程组就能得到解析解（也叫闭式解）。这个解是： a = cov(x, y) / var(x) b = mean(y) - a * mean(x) 其中 cov 是协方差，var 是方差。这个公式很美，因为它直接用数据的统计量表达了结果。

但现实中，我们几乎不会手动去算协方差和方差来求参数。我们会用梯度下降，这是一种更通用、能解决没有解析解问题的方法。它的思想是：随机初始化 a 和 b，然后计算损失函数对它们的梯度（即偏导数），沿着梯度下降最快的方向，用小步长（学习率）一点点更新参数，直到收敛。

import numpy as np

def linear_regression_gd(X, y, lr=0.01, epochs=1000):
    """
    使用梯度下降实现一元线性回归
    X: 特征， shape (n_samples,)
    y: 标签， shape (n_samples,)
    """
    # 初始化参数
    a, b = 0.0, 0.0
    n = len(X)
    
    for epoch in range(epochs):
        # 预测值
        y_pred = a * X + b
        # 计算误差
        error = y_pred - y
        # 计算梯度 (对a和b的偏导)
        grad_a = (2/n) * np.dot(X, error)
        grad_b = (2/n) * np.sum(error)
        # 更新参数
        a -= lr * grad_a
        b -= lr * grad_b
        
        # 每100轮打印一次损失
        if epoch % 100 == 0:
            loss = np.mean(error ** 2)
            print(f"Epoch {epoch}, loss: {loss:.4f}, a: {a:.4f}, b: {b:.4f}")
    return a, b

# 生成示例数据
np.random.seed(42)
X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X + np.random.randn(100, 1)

# 调用我们的梯度下降函数
a, b = linear_regression_gd(X.flatten(), y.flatten(), lr=0.1, epochs=1000)
print(f"\n最终参数: a(斜率)={a:.4f}, b(截距)={b:.4f}")

多元线性回归无非是把 a 变成了一个向量 w（权重），b 变成了标量 b（偏置）。公式变成 y = w1*x1 + w2*x2 + ... + wn*xn + b。此时，用矩阵表示会非常简洁：y = Xw + b。其中 X 是数据矩阵，每一行是一个样本，每一列是一个特征。梯度下降的更新公式也变成了矩阵运算，效率高得多。这里的关键是理解向量化编程，用numpy的矩阵运算代替for循环，速度能快几十倍。

2.2 岭回归与Lasso：正则化如何对抗过拟合

当你用线性回归拟合一个复杂数据集时，可能会发现模型在训练集上表现完美，但在测试集上一塌糊涂。这就是过拟合——模型把数据中的噪声也当规律学进去了。怎么办？正则化来了。

正则化的核心思想是在损失函数中加入一个对模型复杂度的惩罚项。模型越复杂（通常表现为权重 w 的绝对值很大），惩罚就越大。这样，模型就会在“拟合数据”和“保持简单”之间找一个平衡。

岭回归用的是L2正则化，惩罚项是权重的平方和 λ * Σ(w_i²)。它的损失函数是： J(w) = Σ(y_i - ŷ_i)² + λ * Σ(w_i²) 这里的 λ 是超参数，控制惩罚力度。λ 越大，模型越倾向于让权重趋近于0（但不会等于0），模型就越简单。从几何上看，L2正则化相当于给权重向量 w 加了一个圆形约束。

Lasso回归用的是L1正则化，惩罚项是权重的绝对值之和 λ * Σ|w_i|。它的损失函数是： J(w) = Σ(y_i - ŷ_i)² + λ * Σ|w_i| L1正则化更厉害，它不仅能防止过拟合，还能做特征选择。因为它的约束区域是菱形，最优解更容易落在菱形的角上，这使得一部分权重精确地为0。相当于模型自动帮你筛选掉了不重要的特征。

from sklearn.linear_model import Ridge, Lasso
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline

# 假设我们有一个高维数据集，特征很多，样本相对少，容易过拟合
# 使用Pipeline方便地组合标准化和回归模型
ridge_model = make_pipeline(StandardScaler(), Ridge(alpha=1.0)) # alpha就是λ
lasso_model = make_pipeline(StandardScaler(), Lasso(alpha=0.1))

# 拟合模型
# ridge_model.fit(X_train, y_train)
# lasso_model.fit(X_train, y_train)

# 查看Lasso的系数，会发现很多是0
# print(lasso_model.named_steps['lasso'].coef_)

注意：使用L1或L2正则化前，一定要对特征进行标准化。因为惩罚项是对所有权重施加同样的惩罚力度。如果特征A的范围是0-1，特征B的范围是0-1000，那么特征B的权重自然会很小，惩罚项对它的影响就微乎其微，这就不公平了。标准化能让所有特征处于同一量纲。

2.3 逻辑回归：分类世界的“回归”算法

别被名字骗了，逻辑回归是个分类算法，而且是二分类的扛把子。那为什么叫“回归”呢？因为它用线性回归的公式 z = w*x + b 作为输入，然后套上一个Sigmoid函数，把这个线性输出 z 映射到 (0, 1) 之间，解释为概率。

Sigmoid函数长这样：σ(z) = 1 / (1 + e^{-z}) 它的图像是个漂亮的S形曲线。当 z 很大时，输出接近1；当 z 很小时，输出接近0；当 z=0 时，输出正好是0.5。

所以逻辑回归的完整公式是：P(y=1|x) = σ(w*x + b) = 1 / (1 + e^{-(w*x + b)}) 我们预测时，如果 P > 0.5，就认为是正类（1），否则是负类（0）。

训练逻辑回归不能用最小二乘法了，因为标签是0/1。这里我们使用交叉熵损失函数，它衡量的是预测概率分布与真实标签分布的差异。对于单个样本，损失是： Loss = -[y * log(p) + (1-y) * log(1-p)] 当真实标签 y=1 时，损失是 -log(p)，预测概率 p 越接近1，损失越小；当 y=0 时，损失是 -log(1-p)，预测概率 p 越接近0，损失越小。这个函数是凸函数，可以用梯度下降优化。

def sigmoid(z):
    """Sigmoid激活函数"""
    return 1 / (1 + np.exp(-z))

def logistic_regression_gd(X, y, lr=0.01, epochs=30000):
    """
    手动实现逻辑回归梯度下降
    X: 特征矩阵，shape (n_samples, n_features)
    y: 标签，shape (n_samples,)，取值为0或1
    """
    # 初始化参数
    n_samples, n_features = X.shape
    w = np.zeros(n_features)
    b = 0.0
    
    # 梯度下降循环
    for epoch in range(epochs):
        # 线性部分
        z = np.dot(X, w) + b
        # 通过sigmoid得到概率
        p = sigmoid(z)
        
        # 计算梯度
        # 这是交叉熵损失对z的导数，推导后形式非常简单：p - y
        dz = p - y
        dw = (1/n_samples) * np.dot(X.T, dz)
        db = (1/n_samples) * np.sum(dz)
        
        # 更新参数
        w -= lr * dw
        b -= lr * db
        
        if epoch % 5000 == 0:
            # 计算损失
            loss = -np.mean(y * np.log(p + 1e-15) + (1-y) * np.log(1-p + 1e-15))
            # 计算准确率
            y_pred = (p > 0.5).astype(int)
            acc = np.mean(y_pred == y)
            print(f"Epoch {epoch}, Loss: {loss:.4f}, Accuracy: {acc:.4f}")
    return w, b

# 生成简单的二分类数据
from sklearn.datasets import make_classification
X, y = make_classification(n_samples=100, n_features=2, n_informative=2, n_redundant=0, random_state=42)
# 训练我们自己的逻辑回归
w, b = logistic_regression_gd(X, y, lr=0.1, epochs=30000)

3. 树模型与集成学习：从简单到强大的进化

如果说线性模型是“全局”思维，试图用一条直线或平面去拟合所有数据，那么树模型就是“分而治之”的局部思维，它把特征空间切分成一个个矩形区域，在每个区域里做简单的预测（比如平均值）。

3.1 决策树：如何做出一个又一个选择

决策树的核心就两个问题：1. 选哪个特征来分裂？ 2. 在哪里分裂（分裂阈值是多少）？ 答案是要找到一个分裂方式，能让分裂后的子集“纯度”更高。衡量纯度的指标主要有三个：

• 基尼不纯度（Gini Impurity）：随机从数据集中抽两个样本，它们类别不同的概率。概率越低，纯度越高。CART树常用这个。
• 信息增益（Information Gain）：用分裂前的信息熵减去分裂后的加权平均信息熵。信息熵表示混乱程度，减少得越多越好。ID3和C4.5算法用这个。
• 方差减少（Variance Reduction）：用于回归树，目标是让分裂后子集的标签方差最小。

以分类树常用的基尼指数为例，对于一个数据集 D，其基尼指数计算为： Gini(D) = 1 - Σ (p_i)²，其中 p_i 是第 i 类样本在 D 中的比例。 假设我们根据特征 A 的某个阈值 t 将 D 分裂成 D_left 和 D_right，那么分裂后的加权基尼指数为： Gini_split = (|D_left|/|D|) * Gini(D_left) + (|D_right|/|D|) * Gini(D_right) 我们遍历所有特征和所有可能的分裂点，选择使 Gini_split 最小的那个特征和阈值。这个过程递归进行，直到满足停止条件（如树达到最大深度、节点样本数太少等）。

from sklearn.tree import DecisionTreeClassifier, plot_tree
import matplotlib.pyplot as plt

# 使用sklearn的决策树，我们可以直观看到其分裂过程
dt_clf = DecisionTreeClassifier(max_depth=3, criterion='gini', random_state=42)
# 假设用鸢尾花数据集
from sklearn.datasets import load_iris
iris = load_iris()
X_iris, y_iris = iris.data[:, :2], iris.target # 只取前两个特征方便可视化

dt_clf.fit(X_iris, y_iris)

# 可视化决策树
plt.figure(figsize=(12,8))
plot_tree(dt_clf, feature_names=iris.feature_names[:2], class_names=iris.target_names, filled=True)
plt.show()

# 解读：图中每个节点都显示了使用的特征、分裂阈值、基尼指数、样本数和类别分布。
# 例如，第一个节点用“花瓣长度”是否小于等于2.45来分裂。这个分裂点就是算法遍历所有可能值后，选出的能使基尼指数下降最多的点。

决策树很容易过拟合，因为它会一直分裂到每个叶子节点都“纯”为止。所以我们需要剪枝，或者使用它的集成版本。

3.2 随机森林：群体的智慧

“三个臭皮匠，顶个诸葛亮。”随机森林就是这个思想的完美体现。它构建很多棵决策树（比如500棵），然后让它们投票做决定。关键有两点：Bagging 和 随机特征子集。

Bagging（Bootstrap Aggregating）：每棵树训练时，不是用全部数据，而是从原始数据集中有放回地随机抽取一个子集（Bootstrap样本）。这样每棵树看到的数据都略有不同，减少了模型方差。

随机特征子集：每棵树在分裂节点时，不是从所有特征里找最佳分裂点，而是先随机选取一个特征子集（比如总特征数的平方根），然后只在这个子集里找。这进一步增加了树之间的差异性。

最终预测时，分类问题就投票，回归问题就取平均。因为每棵树可能都会犯不同的错误，但大家一平均，错误就抵消了，模型的泛化能力大大增强。随机森林几乎是我做项目的首选基线模型，因为它不容易过拟合，对参数不敏感，还自带特征重要性评估。

from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split

# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X_iris, y_iris, test_size=0.3, random_state=42)

rf_clf = RandomForestClassifier(n_estimators=100, # 树的数量
                                 max_depth=5,      # 每棵树的最大深度，控制单棵树复杂度
                                 max_features='sqrt', # 分裂时考虑的特征数，sqrt是常用选择
                                 random_state=42)
rf_clf.fit(X_train, y_train)
print(f"随机森林在测试集上的准确率: {rf_clf.score(X_test, y_test):.4f}")

# 查看特征重要性
importances = rf_clf.feature_importances_
for name, importance in zip(iris.feature_names[:2], importances):
    print(f"{name}: {importance:.4f}")
# 特征重要性是基于“该特征在森林中所有树上，带来的不纯度减少的总和”来计算的，非常直观有用。

3.3 XGBoost与梯度提升：迭代的艺术

如果说随机森林是“并行”集成（树之间独立训练），那么以XGBoost为代表的梯度提升树（Gradient Boosting Trees） 就是“串行”集成。它的思想更巧妙：每一棵树都在学习前一棵树留下的残差。

想象一下，我们先用一棵很浅的树（比如深度为3）去拟合数据，预测结果肯定不准，会有残差（真实值减去预测值）。然后，我们训练第二棵树，但它的目标不是原始标签 y，而是第一棵树预测的残差。第二棵树拟合了部分残差后，我们把它的预测结果加到第一棵树上，得到一个新的、更准的模型。这个过程不断重复，每一棵新树都致力于纠正前面所有树组合起来的错误。

XGBoost是梯度提升的一个高效、工程化实现。它的核心在于其目标函数的设计。目标函数不仅包含衡量预测好坏的损失函数（如均方误差、交叉熵），还包含正则化项来控制模型的复杂度（叶子节点权重和树结构）。在每一轮迭代中，XGBoost会使用目标函数的二阶泰勒展开（用到了一阶梯度 g 和二阶海森矩阵 h）来近似损失，从而更聪明地确定如何生长一棵树，以及分裂后叶子节点的最优权重是多少。

# 使用XGBoost需要先安装：pip install xgboost
import xgboost as xgb
from sklearn.metrics import accuracy_score

# 将数据转换为XGBoost专用的DMatrix格式，效率更高
dtrain = xgb.DMatrix(X_train, label=y_train)
dtest = xgb.DMatrix(X_test, label=y_test)

# 设置参数
params = {
    'objective': 'multi:softprob',  # 多分类任务
    'num_class': 3,
    'max_depth': 3,                 # 每棵树的最大深度
    'eta': 0.1,                     # 学习率，相当于每棵树的贡献权重
    'subsample': 0.8,               # 每棵树使用的样本比例，类似Bagging
    'colsample_bytree': 0.8,        # 每棵树使用的特征比例
    'seed': 42
}
num_rounds = 50  # 迭代轮数，即树的数量

# 训练模型
bst = xgb.train(params, dtrain, num_rounds)

# 预测
y_pred_prob = bst.predict(dtest)
y_pred = np.argmax(y_pred_prob, axis=1)  # 取概率最大的类别
print(f"XGBoost测试集准确率: {accuracy_score(y_test, y_pred):.4f}")

# XGBoost同样可以输出特征重要性
importance_dict = bst.get_score(importance_type='weight') # 'weight'表示特征被用作分裂点的总次数
print("特征重要性(按分裂次数):", importance_dict)

提示：XGBoost中的 eta 参数（学习率）非常重要。它控制了每棵树对最终模型的贡献程度。较小的 eta（如0.01-0.1）意味着我们需要更多的树（num_rounds）来达到好的效果，但模型会更稳健，不容易过拟合。这被称为Shrinkage（收缩）技术，是提升方法稳定性的关键。

4. 无监督学习与经典算法：发现数据的内在结构

无监督学习没有标签，它的任务是探索数据本身的结构，比如把相似的数据聚在一起（聚类），或者把高维数据压缩到低维以便可视化（降维）。

4.1 K-Means聚类：迭代寻找中心点

K-Means的想法非常直观：我想把数据分成K个簇，那我就先随机找K个点作为“中心点”（质心），然后进行以下两步的迭代：

1. 分配：计算每个数据点到所有质心的距离，把它分配给距离最近的那个质心所在的簇。
2. 更新：重新计算每个簇的质心，方法就是取这个簇里所有数据点的平均值。

重复这两步，直到质心的位置不再变化（或者变化很小）。它的目标是最小化所有数据点到其所属簇质心的距离平方和，这个指标叫簇内平方和（Within-Cluster Sum of Squares, WCSS）。

K-Means最大的挑战是K值的选择和初始质心的敏感性。K值选不对，结果可能没意义。常用的方法是“肘部法则”：画出不同K值对应的WCSS，WCSS会随着K增大而减小，但减小的幅度会突然变缓，那个拐点就像手肘，对应的K值可能是一个好选择。

from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs
import matplotlib.pyplot as plt

# 生成模拟的聚类数据
X_blobs, y_true = make_blobs(n_samples=300, centers=4, cluster_std=0.6, random_state=42)

# 应用K-Means
kmeans = KMeans(n_clusters=4, init='k-means++', n_init=10, random_state=42) # init='k-means++'是智能初始化，比纯随机好
y_kmeans = kmeans.fit_predict(X_blobs)

# 可视化结果
plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.scatter(X_blobs[:, 0], X_blobs[:, 1], c=y_true, cmap='viridis', s=50, alpha=0.7)
plt.title("真实簇分布")
plt.subplot(1, 2, 2)
plt.scatter(X_blobs[:, 0], X_blobs[:, 1], c=y_kmeans, cmap='viridis', s=50, alpha=0.7)
plt.scatter(kmeans.cluster_centers_[:, 0], kmeans.cluster_centers_[:, 1], 
            c='red', marker='X', s=200, label='质心')
plt.title("K-Means聚类结果")
plt.legend()
plt.show()

# 肘部法则选择K
wcss = []
for k in range(1, 11):
    kmeans_test = KMeans(n_clusters=k, init='k-means++', n_init=10, random_state=42)
    kmeans_test.fit(X_blobs)
    wcss.append(kmeans_test.inertia_)  # inertia_属性就是WCSS

plt.figure(figsize=(8,5))
plt.plot(range(1, 11), wcss, marker='o')
plt.xlabel('簇的数量 K')
plt.ylabel('簇内平方和 (WCSS)')
plt.title('肘部法则')
plt.grid(True)
plt.show()
# 从图中可以看到，当K=4时，WCSS下降的斜率发生明显变化，这就是“肘部”。

4.2 DBSCAN：基于密度的聚类高手

K-Means假设簇是球形的，并且大小差不多。但现实中的数据簇可能是任意形状的，比如一个圆环套着另一个圆环。这时就需要DBSCAN。

DBSCAN的核心概念就三个：核心点、边界点和噪声点。

• 核心点：在它的邻域半径 eps 内，至少有 min_samples 个其他点。
• 边界点：在某个核心点的邻域内，但自身不是核心点。
• 噪声点：既不是核心点，也不是边界点。

算法过程就是：随机找一个未访问的核心点，把它和所有由它密度可达的点（包括其他核心点和边界点）划为一个簇。然后找下一个未访问的核心点，重复这个过程。所有没被划入任何簇的点就是噪声。

DBSCAN的强大之处在于它能发现任意形状的簇，并且能识别噪声。但它对参数 eps 和 min_samples 比较敏感，且不适用于密度差异很大的数据集。

from sklearn.cluster import DBSCAN
from sklearn.datasets import make_moons

# 生成半月形数据，K-Means对此无能为力
X_moons, y_moons = make_moons(n_samples=300, noise=0.05, random_state=42)

# 尝试用K-Means，效果会很差
kmeans_moons = KMeans(n_clusters=2, random_state=42)
y_kmeans_moons = kmeans_moons.fit_predict(X_moons)

# 使用DBSCAN
dbscan = DBSCAN(eps=0.2, min_samples=5)
y_dbscan = dbscan.fit_predict(X_moons)

# 可视化对比
plt.figure(figsize=(12, 4))
plt.subplot(1, 3, 1)
plt.scatter(X_moons[:, 0], X_moons[:, 1], c=y_moons, cmap='viridis', s=50)
plt.title("真实分布（半月形）")
plt.subplot(1, 3, 2)
plt.scatter(X_moons[:, 0], X_moons[:, 1], c=y_kmeans_moons, cmap='viridis', s=50)
plt.title("K-Means聚类结果")
plt.subplot(1, 3, 3)
plt.scatter(X_moons[:, 0], X_moons[:, 1], c=y_dbscan, cmap='viridis', s=50)
# DBSCAN中标签为-1的是噪声点，我们可以用不同标记标出
noise_mask = (y_dbscan == -1)
plt.scatter(X_moons[noise_mask, 0], X_moons[noise_mask, 1], c='black', marker='x', s=100, label='噪声')
plt.title("DBSCAN聚类结果")
plt.legend()
plt.show()
# 可以看到DBSCAN完美地分开了两个半月形，而K-Means则按空间位置切了一刀。

4.3 PCA主成分分析：化繁为简的魔法

我们经常遇到成百上千个特征的数据集，不仅计算慢，而且可能存在多重共线性。PCA就是一种降维技术，它通过线性变换，把原始高维数据投影到新的低维坐标系中，并且要求新坐标系的轴（主成分）能最大程度地保留原始数据的方差（即信息）。

数学上，PCA的步骤是：

1. 中心化：每个特征减去其均值，使数据以原点为中心。
2. 计算协方差矩阵：协方差矩阵反映了特征之间的线性关系。
3. 特征值分解：对协方差矩阵进行特征值分解，得到特征值和特征向量。特征值的大小代表了对应特征向量方向上方差的大小。
4. 选择主成分：将特征值从大到小排序，选择前k个最大的特征值对应的特征向量，组成投影矩阵 W。
5. 降维：将原始数据 X 乘以投影矩阵 W，得到降维后的数据 X_pca = X * W。

第一个主成分是方差最大的方向，第二个主成分是与第一个正交的、方差次大的方向，以此类推。PCA之后的数据，各维度之间是线性无关的（协方差为0）。

from sklearn.decomposition import PCA
from sklearn.datasets import load_iris

iris = load_iris()
X = iris.data
y = iris.target

# 应用PCA，降到2维以便可视化
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X)  # 这一步同时完成了拟合（计算投影矩阵）和转换（降维）

print(f"原始数据形状: {X.shape}")
print(f"降维后数据形状: {X_pca.shape}")
print(f"各主成分解释的方差比例: {pca.explained_variance_ratio_}")
print(f"前两个主成分累计解释方差: {pca.explained_variance_ratio_.sum():.4f}")

# 可视化降维后的数据
plt.figure(figsize=(8,6))
scatter = plt.scatter(X_pca[:, 0], X_pca[:, 1], c=y, cmap='viridis', s=70, alpha=0.8)
plt.xlabel('第一主成分 (PC1)')
plt.ylabel('第二主成分 (PC2)')
plt.title('鸢尾花数据集PCA降维可视化')
plt.colorbar(scatter, label='鸢尾花类别')
plt.grid(True)
plt.show()

# 查看主成分（特征向量）与原始特征的关系
print("\n主成分（特征向量）:")
for i, component in enumerate(pca.components_):
    print(f"PC{i+1}: {component}")
# 例如，PC1 = 0.36*花萼长 + -0.08*花萼宽 + 0.86*花瓣长 + 0.36*花瓣宽
# 这说明第一主成分主要由花瓣长度主导。

理解PCA的数学推导，关键在于理解特征值分解和方差最大化的优化目标。它本质上是在寻找一个正交基，使得数据在这个新基下的表示，其协方差矩阵是对角矩阵（各维度无关），并且对角线上的值（方差）尽可能大。