1. 项目概述：从理论到代码，打通预测模型的任督二脉

在数学建模竞赛和实际数据分析项目中，预测模型是当之无愧的“顶流”。无论是预测销量、股价，还是分析用户行为、评估风险，一个得心应手的预测模型库就是你的“瑞士军刀”。然而，很多同学和初入行的朋友常常面临一个尴尬的断层：理论课上听得头头是道，公式推导也能看懂，但一到自己动手，面对空白的代码编辑器就两眼一抹黑，不知道如何把那些优美的数学公式变成一行行能跑出结果的代码。更常见的是，虽然能调用 sklearn 的几行接口跑出结果，但对模型内部的参数调整、结果评估的深层逻辑、乃至不同模型间的核心差异，依然是一知半解。

这篇内容，就是为你弥合这道鸿沟而准备的。我们不空谈理论，也不只做简单的代码搬运。我将结合自己多年打比赛和做项目的经验，带你深入线性回归、决策树、随机森林这几个最经典、最常用的预测模型内部。目标是让你不仅“知其然”（知道代码怎么写），更能“知其所以然”（理解为什么这么写，以及背后的数学和统计逻辑）。我们会从最基础的模型假设和原理讲起，然后立刻用 Python 代码实现其核心计算过程，最后再与成熟的机器学习库（如 scikit-learn ）进行对比和衔接。你会发现，自己动手实现一遍，对于理解模型、调试参数、甚至创新改进，有着不可替代的价值。无论你是正在备战数学建模竞赛的学生，还是希望夯实基础的算法工程师，这篇文章都能提供一条从理论直通实战的清晰路径。

2. 基石模型：线性回归的原理拆解与手动实现

线性回归堪称预测世界的“Hello World”。它形式简单，却蕴含着最小二乘估计、最大似然估计等重要的统计思想。很多人觉得它太简单而忽视其深度，但恰恰是吃透线性回归，才能为理解更复杂的模型打下坚实基础。

2.1 模型核心：最小二乘的几何与代数视角

线性回归试图用一个线性方程来拟合数据： y = w1*x1 + w2*x2 + ... + wn*xn + b 。其中， y 是因变量（我们要预测的）， x1, x2,..., xn 是自变量（特征）， w1, w2,..., wn 是权重（系数）， b 是截距。我们的目标是找到一组 w 和 b ，使得预测值 y_hat 与真实值 y 之间的差距最小。

这个“差距最小”在数学上通常定义为 残差平方和（RSS）最小 ，即最小二乘法。从代数角度看，我们是在求解一个优化问题。从几何角度看，我们是在高维空间中，寻找一个超平面，使得所有数据点到这个超平面的垂直距离（残差）的平方和最小。这个解有一个漂亮的闭式解（解析解），可以通过矩阵运算直接求得： W = (X^T * X)^(-1) * X^T * y 。这里 X 是增加了全为1的列（对应截距b）的特征矩阵， W 是包含所有权重和截距的向量。

这个公式很美，但它隐藏了两个关键前提：一是 X^T * X 矩阵必须是可逆的（满秩），这意味着特征之间不能存在完全的线性相关（即多重共线性）；二是它假设误差项服从均值为0、方差恒定的正态分布，且相互独立。在实际中，这些假设经常被违反，这就引出了岭回归（Ridge）、Lasso回归等正则化变体。

2.2 从零实现：梯度下降与矩阵求导

虽然闭式解很优雅，但在特征维度很高（ n 很大）时，计算矩阵的逆 (X^T * X)^(-1) 会非常耗时，甚至可能因为矩阵病态而无法计算。因此，工业界和实际编程中更常用的是 梯度下降法 ，这是一种迭代逼近的数值优化方法。

梯度下降的核心思想很直观：想象你站在一座山上，要最快下到山谷（找到损失函数的最小值）。你每走一步，都沿着当前所在位置最陡峭的下山方向（负梯度方向）前进。步长就是学习率。对于线性回归，损失函数是均方误差（MSE）： L = (1/m) * Σ(y_i - y_hat_i)^2 。对其求关于权重 W 的偏导数，得到梯度： ∇L = -(2/m) * X^T * (y - X*W) 。

下面，我们抛开 sklearn ，用 NumPy 从头实现一个基于批量梯度下降的线性回归：

import numpy as np

class LinearRegressionGD:
    def __init__(self, learning_rate=0.01, n_iters=1000):
        """
        初始化线性回归模型。
        :param learning_rate: 学习率，控制每一步更新的幅度。太大可能震荡，太小收敛慢。
        :param n_iters: 梯度下降迭代次数。
        """
        self.lr = learning_rate
        self.n_iters = n_iters
        self.weights = None # 权重向量（包含截距）
        self.loss_history = [] # 记录每次迭代的损失，用于可视化检查收敛

    def _add_intercept(self, X):
        """ 在特征矩阵X前添加一列1，用于计算截距项b。"""
        intercept = np.ones((X.shape[0], 1))
        return np.hstack((intercept, X))

    def fit(self, X, y):
        """
        训练模型，使用梯度下降法拟合数据。
        :param X: 训练特征，形状 (m_samples, n_features)
        :param y: 训练标签，形状 (m_samples,)
        """
        # 添加截距项
        X_b = self._add_intercept(X)
        m, n = X_b.shape

        # 初始化权重，通常用小的随机数或零初始化
        self.weights = np.random.randn(n) * 0.01

        for i in range(self.n_iters):
            # 计算当前预测值
            y_pred = X_b.dot(self.weights)
            # 计算误差
            error = y_pred - y
            # 计算梯度：X_b.T.dot(error) / m
            gradients = (2 / m) * X_b.T.dot(error)
            # 更新权重：向负梯度方向移动
            self.weights -= self.lr * gradients

            # 计算并记录当前损失（MSE）
            loss = np.mean(error ** 2)
            self.loss_history.append(loss)

            # 可选：每100次迭代打印一次损失，监控训练过程
            if i % 100 == 0:
                print(f"Iteration {i}: Loss {loss:.4f}")

    def predict(self, X):
        """ 使用训练好的权重进行预测。"""
        X_b = self._add_intercept(X)
        return X_b.dot(self.weights)

# 使用示例
if __name__ == "__main__":
    # 生成一些简单的线性数据，加一点噪声
    np.random.seed(42)
    m = 100
    X = 2 * np.random.rand(m, 1) # 一个特征
    y = 4 + 3 * X + np.random.randn(m, 1) # 真实关系: y = 4 + 3x + noise

    # 实例化并训练模型
    lr_model = LinearRegressionGD(learning_rate=0.1, n_iters=1000)
    lr_model.fit(X, y.flatten()) # y需要展平

    # 查看学到的参数 (截距和权重)
    print(f"Learned intercept (b): {lr_model.weights[0]:.4f}")
    print(f"Learned coefficient (w): {lr_model.weights[1]:.4f}")
    # 预测新数据
    X_new = np.array([[1.5]])
    print(f"Prediction for X=1.5: {lr_model.predict(X_new)[0]:.4f}")

注意 ：手动实现梯度下降时，学习率 learning_rate 的选择至关重要。上述代码用了0.1，对于这个简单数据集可行。但对于不同尺度的数据，可能需要先进行特征标准化（如 StandardScaler ），否则梯度下降可能难以收敛。这是手动实现教会我们的第一课： 数据预处理与模型训练密不可分 。

2.3 与Sklearn对比：理解封装背后的逻辑

运行上面的代码，你会得到接近真实值（4和3）的参数。现在，我们用 sklearn 做同样的事：

from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error

sk_model = LinearRegression()
sk_model.fit(X, y)
print(f"Sklearn intercept: {sk_model.intercept_[0]:.4f}")
print(f"Sklearn coefficient: {sk_model.coef_[0][0]:.4f}")

y_pred_sk = sk_model.predict(X)
y_pred_manual = lr_model.predict(X)

print(f"Sklearn MSE: {mean_squared_error(y, y_pred_sk):.6f}")
print(f"Manual GD MSE: {mean_squared_error(y, y_pred_manual):.6f}")

你会发现， sklearn 的结果通常更精确、更稳定，因为它内部可能使用更高效的数值算法（如基于奇异值分解SVD的解析解计算），并且做了大量的数值稳定性处理。我们自己实现的梯度下降版本，MSE可能会略高一点，这是迭代逼近的特性。

这里的核心收获是 ： sklearn 的 .fit() 和 .predict() 只是一个高度封装的黑箱。通过手动实现，你明白了这个黑箱里在发生什么：它可能在计算梯度、迭代更新权重、并监控损失。当你在使用 sklearn 的 LinearRegression 遇到奇异矩阵报错时，你就能立刻想到可能是多重共线性问题，从而考虑使用 Ridge 回归（它在损失函数中加入了权重的L2范数惩罚项，使得 X^T * X + αI 总是可逆）。这种从内到外的理解，是调参和模型选择的底气所在。

3. 树模型初探：决策树的核心分裂逻辑与实现

如果说线性回归是试图用一条直线（或平面）去拟合世界，那么决策树则是用一系列“如果...那么...”的规则来划分世界。它更直观，更容易解释，也能处理非线性关系。决策树的核心在于“分裂”：如何选择最佳的特征和分割点，将数据分成更“纯”的子集。

3.1 分裂准则：信息增益、基尼系数与方差减少

决策树的学习过程是一个递归的“选择特征-分割数据”的过程。关键问题是：我们依据什么标准来选择“最好”的分割？对于 分类树 ，常用指标有：

1. 信息增益（Information Gain） ：基于信息熵。熵表示集合的混乱程度。信息增益是父节点的熵减去子节点熵的加权和。增益越大，说明用该特征分割后，子集的“纯度”提升越多。ID3算法使用此准则。
2. 基尼不纯度（Gini Impurity） ：衡量从数据集中随机抽取两个样本，其类别标签不一致的概率。基尼系数越小，纯度越高。CART分类树使用此准则。

对于 回归树 （预测连续值），常用指标是： 3. 方差减少（Variance Reduction） ：选择那个能使分割后子节点目标值方差加权和减少最多的特征和分割点。直观理解就是让每个子组内的数据尽可能相似。

以基尼系数为例，其计算公式为： Gini(p) = 1 - Σ (p_i)^2 ，其中 p_i 是第 i 类在节点中出现的概率。假设一个节点有10个样本，6个属于A类，4个属于B类。则该节点的基尼系数为 1 - ((6/10)^2 + (4/10)^2) = 1 - (0.36 + 0.16) = 0.48 。如果我们找到一个分割点，将数据分成两个子节点：左节点（8个样本，6个A，2个B）和右节点（2个样本，0个A，2个B）。我们可以计算分割后的加权基尼系数，并与父节点的0.48比较，得出基尼系数的减少量，即该分割的“增益”。

3.2 手动构建一棵简单的CART分类树

实现一棵完整的、支持剪枝的决策树比较复杂，但我们可以实现其核心：根据基尼系数寻找最佳分割。这能让你彻底理解 .fit() 方法在做什么。

import numpy as np
from collections import Counter

class SimpleDecisionTree:
    """ 一个极简的CART分类树实现，仅用于演示核心分裂逻辑。"""
    def __init__(self, max_depth=5, min_samples_split=2):
        self.max_depth = max_depth
        self.min_samples_split = min_samples_split
        self.tree = None

    def _gini(self, y):
        """ 计算一个节点中样本集合的基尼不纯度。"""
        counter = Counter(y)
        m = len(y)
        if m == 0:
            return 0
        gini = 1.0
        for count in counter.values():
            p = count / m
            gini -= p ** 2
        return gini

    def _split(self, X_column, split_val):
        """ 根据特征列和分割值，返回左右子集的索引。"""
        left_indices = np.where(X_column <= split_val)[0]
        right_indices = np.where(X_column > split_val)[0]
        return left_indices, right_indices

    def _best_split(self, X, y):
        """ 遍历所有特征和所有可能的分割点，找到最佳分割。"""
        best_gini_gain = -1
        best_feature_idx = None
        best_split_val = None
        best_left_idx, best_right_idx = None, None

        m, n = X.shape
        if m <= 1:
            return None, None, None, None # 无法继续分割

        parent_gini = self._gini(y)

        for feature_idx in range(n):
            # 获取该特征列的所有唯一值，并排序，取相邻值的中点作为候选分割点
            feature_vals = np.unique(X[:, feature_idx])
            split_points = (feature_vals[:-1] + feature_vals[1:]) / 2

            for split_val in split_points:
                left_idx, right_idx = self._split(X[:, feature_idx], split_val)

                if len(left_idx) == 0 or len(right_idx) == 0:
                    continue # 分割无效，跳过

                # 计算加权子节点基尼系数
                left_gini = self._gini(y[left_idx])
                right_gini = self._gini(y[right_idx])
                weighted_gini = (len(left_idx)/m)*left_gini + (len(right_idx)/m)*right_gini

                # 计算基尼增益
                gini_gain = parent_gini - weighted_gini

                # 记录最佳增益
                if gini_gain > best_gini_gain:
                    best_gini_gain = gini_gain
                    best_feature_idx = feature_idx
                    best_split_val = split_val
                    best_left_idx, best_right_idx = left_idx, right_idx

        # 如果增益太小或为负，则认为没有有效分割
        if best_gini_gain < 1e-7:
            return None, None, None, None

        return best_feature_idx, best_split_val, best_left_idx, best_right_idx

    def _build_tree(self, X, y, depth=0):
        """ 递归构建决策树。"""
        # 终止条件：达到最大深度、样本数过少、或节点纯度已很高（基尼系数为0）
        if (depth >= self.max_depth or
            len(y) < self.min_samples_split or
            self._gini(y) < 1e-7):
            # 返回叶子节点，其值为该节点中最常见的类别
            most_common = Counter(y).most_common(1)[0][0]
            return {'type': 'leaf', 'value': most_common}

        # 寻找最佳分割
        feat_idx, split_val, left_idx, right_idx = self._best_split(X, y)

        # 如果找不到有效分割，也作为叶子节点
        if feat_idx is None:
            most_common = Counter(y).most_common(1)[0][0]
            return {'type': 'leaf', 'value': most_common}

        # 递归构建左右子树
        left_tree = self._build_tree(X[left_idx], y[left_idx], depth+1)
        right_tree = self._build_tree(X[right_idx], y[right_idx], depth+1)

        # 返回一个内部节点
        return {
            'type': 'internal',
            'feature_idx': feat_idx,
            'split_val': split_val,
            'left': left_tree,
            'right': right_tree
        }

    def fit(self, X, y):
        """ 训练决策树。"""
        self.tree = self._build_tree(np.array(X), np.array(y))

    def _predict_one(self, x, node):
        """ 对单个样本进行预测。"""
        if node['type'] == 'leaf':
            return node['value']
        # 根据节点的分裂规则，决定走左子树还是右子树
        if x[node['feature_idx']] <= node['split_val']:
            return self._predict_one(x, node['left'])
        else:
            return self._predict_one(x, node['right'])

    def predict(self, X):
        """ 对数据集进行预测。"""
        return np.array([self._predict_one(x, self.tree) for x in np.array(X)])

# 使用示例：鸢尾花数据集
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

iris = load_iris()
X, y = iris.data, iris.target
# 为了简化，我们只取前两个特征和两个类别
X = X[y != 2, :2]
y = y[y != 2]

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

simple_tree = SimpleDecisionTree(max_depth=3)
simple_tree.fit(X_train, y_train)
y_pred_simple = simple_tree.predict(X_test)
print(f"Simple Tree Accuracy: {accuracy_score(y_test, y_pred_simple):.4f}")

# 与sklearn的决策树对比
from sklearn.tree import DecisionTreeClassifier
sk_tree = DecisionTreeClassifier(max_depth=3, criterion='gini', random_state=42)
sk_tree.fit(X_train, y_train)
y_pred_sk = sk_tree.predict(X_test)
print(f"Sklearn Tree Accuracy: {accuracy_score(y_test, y_pred_sk):.4f}")

这段代码实现了一个非常基础的决策树。它清晰地展示了递归分裂的过程：从根节点开始，计算所有可能分割的基尼增益，选择增益最大的那个特征和分割点，将数据分成两份，然后对左右子集递归地重复这个过程，直到满足停止条件（深度、最小样本数、纯度）。

实操心得 ：自己实现一遍后，你会深刻理解决策树的两个关键超参数 max_depth 和 min_samples_split 的意义。 max_depth 控制树的复杂度，防止过拟合； min_samples_split 规定了一个节点必须有多少样本才允许继续分裂，避免产生只包含极少样本的无效节点。在数学建模中，理解这些参数比盲目调参更重要。

3.3 决策树的优势、劣势与可视化

决策树最大的优点是 可解释性极强 。你可以直接把训练好的树画出来，看到从根到叶子的每一条判断路径，这非常符合人类的决策逻辑。在需要向非技术人员解释模型原因的场合（如金融风控、医疗诊断辅助），这是巨大的优势。

但其劣势也很明显： 非常容易过拟合 。一棵不加限制的树会一直分裂，直到每个叶子节点都只包含一个样本（纯度100%），这相当于把训练数据完全背了下来，对噪声极度敏感，在新数据上表现会很差。这就是为什么我们需要剪枝（Pruning），或者使用它的集成版本——随机森林。

你可以使用 sklearn.tree.plot_tree 或更强大的 graphviz 库来可视化你的树，直观地检查模型的学习逻辑。

4. 集成力量：随机森林的“群众智慧”与并行实现

“三个臭皮匠，顶个诸葛亮。”随机森林就是这句话在机器学习中的完美体现。它通过构建多棵决策树，并将它们的预测结果进行综合（分类用投票，回归用平均），来获得比单棵决策树更稳定、更准确的模型。其核心思想是 Bootstrap Aggregating（Bagging） 和 随机特征子空间 。

4.1 Bagging与随机性：为何有效？

Bagging的核心步骤是：

1. Bootstrap抽样 ：从原始训练集中有放回地随机抽取 m 个样本（通常 m 等于原始集大小），形成一个自助采样集。由于是有放回，一些样本可能被抽到多次，一些样本可能一次都没被抽到。未被抽到的样本称为“袋外样本”（Out-Of-Bag, OOB），可用于模型验证。
2. 并行训练 ：用这个采样集独立训练一个基学习器（比如一棵决策树）。
3. 重复与聚合 ：重复以上过程 n_estimators 次，得到 n 个基学习器。对于预测，分类问题采用投票法，回归问题采用平均法。

随机森林在Bagging的基础上，增加了一层随机性：在每棵决策树进行节点分裂时，不是从所有 d 个特征中选择最优特征，而是先随机选取一个特征子集（通常大小为 sqrt(d) 或 log2(d) ），然后只在这个子集中寻找最优分裂特征。

这两重随机性带来了巨大好处 ：

• 降低方差 ：通过平均多棵树的预测，平滑了单棵树因数据扰动而产生的方差，有效抑制了过拟合。
• 提升泛化能力 ：特征随机性使得每棵树关注数据的不同侧面，树之间的差异性（多样性）增大，集成的效果更好。
• 提供OOB估计 ：无需额外的验证集，即可用OOB样本评估模型性能，这是一个非常方便的内置交叉验证。

4.2 手动实现一个基础随机森林（回归）

下面我们实现一个回归版本的随机森林，重点展示Bagging和特征随机性的过程。为了简化，我们使用 sklearn 的决策树作为基学习器，但融入我们自己的随机采样逻辑。

import numpy as np
from sklearn.tree import DecisionTreeRegressor
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split
from sklearn.datasets import make_regression

class SimpleRandomForestRegressor:
    def __init__(self, n_estimators=100, max_depth=None, max_features='sqrt', random_state=None):
        """
        初始化随机森林回归器。
        :param n_estimators: 森林中树的数量。
        :param max_depth: 每棵树的最大深度。
        :param max_features: 每棵树分裂时考虑的最大特征数。可以是'int', 'float', 'sqrt', 'log2'。
        :param random_state: 随机种子，用于复现结果。
        """
        self.n_estimators = n_estimators
        self.max_depth = max_depth
        self.max_features = max_features
        self.random_state = random_state
        self.trees = []
        self.feature_indices_for_trees = [] # 记录每棵树使用的特征索引，用于理解模型
        self.oob_predictions = None # 存储OOB预测，用于评估

    def _get_max_features(self, n_features):
        """ 根据参数确定每棵树实际使用的特征数量。"""
        if isinstance(self.max_features, int):
            return min(self.max_features, n_features)
        elif isinstance(self.max_features, float):
            return max(1, int(self.max_features * n_features))
        elif self.max_features == 'sqrt':
            return int(np.sqrt(n_features))
        elif self.max_features == 'log2':
            return int(np.log2(n_features))
        else:
            return n_features # 默认使用所有特征（这时就退化为Bagging了）

    def fit(self, X, y):
        """
        训练随机森林。
        核心：为每棵树进行Bootstrap采样，并记录其使用的特征子集。
        """
        X = np.array(X)
        y = np.array(y).flatten()
        n_samples, n_features = X.shape
        rng = np.random.RandomState(self.random_state)

        # 确定每棵树使用的特征数
        n_features_per_tree = self._get_max_features(n_features)

        # 初始化OOB预测容器：一个列表的列表，每个内列表存储一个样本的所有OOB预测
        oob_pred_list = [[] for _ in range(n_samples)]

        for i in range(self.n_estimators):
            # 1. Bootstrap采样
            # 有放回地随机抽取样本索引
            sample_indices = rng.choice(n_samples, size=n_samples, replace=True)
            # 获取未被抽中的样本索引（OOB样本）
            oob_indices = np.setdiff1d(np.arange(n_samples), np.unique(sample_indices))

            # 2. 随机选择特征子集
            feature_indices = rng.choice(n_features, size=n_features_per_tree, replace=False)
            feature_indices.sort() # 排序便于后续索引
            self.feature_indices_for_trees.append(feature_indices)

            # 3. 准备该树的训练数据（仅使用选中的特征）
            X_bootstrap = X[sample_indices][:, feature_indices]
            y_bootstrap = y[sample_indices]

            # 4. 训练一棵决策树（回归树）
            tree = DecisionTreeRegressor(max_depth=self.max_depth, random_state=rng)
            tree.fit(X_bootstrap, y_bootstrap)
            self.trees.append((tree, feature_indices))

            # 5. 用这棵树对它的OOB样本进行预测，并存储
            if len(oob_indices) > 0:
                X_oob = X[oob_indices][:, feature_indices]
                y_oob_pred = tree.predict(X_oob)
                for idx, pred in zip(oob_indices, y_oob_pred):
                    oob_pred_list[idx].append(pred)

        # 计算每个样本的OOB预测（平均所有包含该样本的树的预测）
        self.oob_predictions = np.zeros(n_samples)
        self.oob_sample_count = np.zeros(n_samples)
        for i in range(n_samples):
            if oob_pred_list[i]: # 如果该样本至少在一棵树的OOB中
                self.oob_predictions[i] = np.mean(oob_pred_list[i])
                self.oob_sample_count[i] = len(oob_pred_list[i])

    def predict(self, X):
        """ 对输入X进行预测，取所有树预测结果的平均值。"""
        X = np.array(X)
        n_samples = X.shape[0]
        all_predictions = np.zeros((self.n_estimators, n_samples))

        for i, (tree, feat_idx) in enumerate(self.trees):
            # 每棵树只使用它训练时对应的特征子集
            X_subset = X[:, feat_idx]
            all_predictions[i, :] = tree.predict(X_subset)

        # 对所有树的预测结果按行取平均
        return np.mean(all_predictions, axis=0)

    def score_oob(self, y_true):
        """ 使用OOB预测计算R^2分数（仅针对有OOB预测的样本）。"""
        mask = self.oob_sample_count > 0
        if np.sum(mask) == 0:
            print("Warning: No OOB samples found. Consider increasing n_estimators.")
            return None
        y_true_masked = y_true[mask]
        y_pred_masked = self.oob_predictions[mask]
        # 计算R^2
        ss_res = np.sum((y_true_masked - y_pred_masked) ** 2)
        ss_tot = np.sum((y_true_masked - np.mean(y_true_masked)) ** 2)
        r2 = 1 - (ss_res / ss_tot)
        return r2

# 生成一个回归数据集进行测试
X, y = make_regression(n_samples=500, n_features=10, noise=0.1, random_state=42)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# 使用我们的简单随机森林
rf_manual = SimpleRandomForestRegressor(n_estimators=50, max_depth=5, max_features='sqrt', random_state=42)
rf_manual.fit(X_train, y_train)

# OOB评分
oob_r2 = rf_manual.score_oob(y_train)
print(f"Manual Random Forest OOB R^2 Score: {oob_r2:.4f}")

# 测试集预测
y_pred_manual = rf_manual.predict(X_test)
test_mse_manual = mean_squared_error(y_test, y_pred_manual)
print(f"Manual Random Forest Test MSE: {test_mse_manual:.6f}")

# 与sklearn的随机森林对比
from sklearn.ensemble import RandomForestRegressor
rf_sk = RandomForestRegressor(n_estimators=50, max_depth=5, max_features='sqrt',
                               oob_score=True, random_state=42)
rf_sk.fit(X_train, y_train)
print(f"Sklearn Random Forest OOB R^2 Score: {rf_sk.oob_score_:.4f}")
y_pred_sk = rf_sk.predict(X_test)
test_mse_sk = mean_squared_error(y_test, y_pred_sk)
print(f"Sklearn Random Forest Test MSE: {test_mse_sk:.6f}")

这个手动实现版本虽然功能上远不如 sklearn 的 RandomForestRegressor 完整和优化（例如缺少对缺失值的处理、并行计算等），但它清晰地揭示了随机森林的工作原理：循环创建多棵树，每棵树基于不同的数据子集和特征子集进行训练，最终预测时集体投票或取平均。

4.3 随机森林的调参要点与特征重要性

理解了原理，调参就不再是玄学。随机森林的关键参数包括：

• n_estimators ：树的数量。越多越好，但计算成本也越高。通常在一定数量后（如100-200），性能提升会趋于平缓。
• max_depth ：单棵树的最大深度。限制深度可以防止过拟合，但可能欠拟合。通常通过交叉验证选择。
• max_features ：每棵树分裂时考虑的最大特征数。这是控制树之间相关性的关键参数。较小的值（如 sqrt 或 log2 ）能增加多样性，可能提升效果，但太小也可能导致信息不足。
• min_samples_split / min_samples_leaf ：分裂节点所需的最小样本数/叶节点所需的最小样本数。增大这些值可以正则化模型，防止过拟合。

随机森林还有一个极其有用的副产品： 特征重要性 。其计算方式通常基于基尼重要性或平均不纯度减少。简单说，一个特征在所有树的所有分裂节点中，被用于分裂时所带来的不纯度减少量的总和越大，它的重要性就越高。 sklearn 通过 model.feature_importances_ 属性直接提供。这在数学建模中用于特征筛选和理解问题驱动因素时，价值连城。

# 接上例，使用sklearn模型获取特征重要性
importances = rf_sk.feature_importances_
indices = np.argsort(importances)[::-1] # 按重要性降序排列

print("Feature ranking:")
for i, idx in enumerate(indices[:5]): # 打印前5个重要特征
    print(f"{i+1}. Feature {idx} ({importances[idx]:.4f})")

5. 模型评估与选择：超越准确率的思考

实现模型只是第一步，如何评估和选择模型才是建模工作的核心。在数学建模竞赛和实际项目中，切忌只看一个指标（如准确率、R平方）。

5.1 回归任务评估矩阵

对于回归问题，常用的评估指标有：

• 均方误差（MSE） ： (1/n) * Σ(y_i - y_hat_i)^2 。平方项放大了大误差的惩罚，是最常用的损失函数。
• 均方根误差（RMSE） ： sqrt(MSE) 。其量纲与原始数据一致，更易于解释。
• 平均绝对误差（MAE） ： (1/n) * Σ|y_i - y_hat_i| 。对异常值不如MSE敏感。
• 决定系数（R²） ： 1 - (SS_res / SS_tot) 。表示模型解释的数据方差比例，越接近1越好。但要注意，在特征很多时，R²会天然偏高。

在 sklearn 中，可以方便地调用这些指标：

from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score

mse = mean_squared_error(y_true, y_pred)
rmse = np.sqrt(mse)
mae = mean_absolute_error(y_true, y_pred)
r2 = r2_score(y_true, y_pred)

5.2 分类任务评估矩阵

对于分类问题，准确率（Accuracy）常常具有误导性，特别是在类别不平衡的数据集上（比如99%的样本是A类，1%是B类，一个把所有样本都预测为A类的傻瓜模型也有99%的准确率）。

• 混淆矩阵（Confusion Matrix） ：这是所有评估的基础。它展示了真正例（TP）、假正例（FP）、真反例（TN）、假反例（FN）的数量。
• 精确率（Precision） ： TP / (TP + FP) 。在所有预测为正的样本中，真正为正的比例。关注预测的“准不准”。
• 召回率（Recall） ： TP / (TP + FN) 。在所有真实为正的样本中，被正确预测出来的比例。关注“找得全不全”。
• F1分数（F1-Score） ：精确率和召回率的调和平均数， 2 * (Precision * Recall) / (Precision + Recall) 。是综合考量。
• ROC曲线与AUC ：通过变化分类阈值，绘制真正例率（TPR） vs. 假正例率（FPR）的曲线，其下面积（AUC）用于评估模型整体排序能力，对类别不平衡不敏感。

from sklearn.metrics import confusion_matrix, classification_report, roc_auc_score

# 对于二分类
print(confusion_matrix(y_true, y_pred))
print(classification_report(y_true, y_pred))
# 注意：roc_auc_score通常需要预测的概率值，而非类别标签
# auc = roc_auc_score(y_true, y_pred_proba)

5.3 交叉验证：稳健评估的黄金准则

永远不要只在一个固定的训练-测试集上评估模型，这有很大的偶然性。 K折交叉验证（K-Fold Cross Validation） 是标准做法。它将数据分成K份，轮流将其中一份作为验证集，其余作为训练集，重复K次，最后取K次评估结果的平均值。

from sklearn.model_selection import cross_val_score
from sklearn.ensemble import RandomForestClassifier

model = RandomForestClassifier(n_estimators=100, random_state=42)
# 使用准确率作为评估指标，进行5折交叉验证
cv_scores = cross_val_score(model, X, y, cv=5, scoring='accuracy')
print(f"Cross-validation scores: {cv_scores}")
print(f"Mean CV accuracy: {cv_scores.mean():.4f} (+/- {cv_scores.std()*2:.4f})")

交叉验证的结果更能反映模型的泛化能力。在数学建模论文中，汇报交叉验证的结果比单一划分的结果更有说服力。

5.4 模型选择实战：以波士顿房价数据集为例

让我们用一个完整的流程，对比线性回归、决策树和随机森林在一个经典数据集上的表现。

import numpy as np
import pandas as pd
from sklearn.datasets import fetch_openml
from sklearn.model_selection import train_test_split, cross_val_score, GridSearchCV
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error, r2_score
import warnings
warnings.filterwarnings('ignore')

# 加载波士顿房价数据集（注意：由于伦理问题，sklearn已移除，这里从openml获取）
boston = fetch_openml(name='boston', version=1, as_frame=True, parser='pandas')
X = boston.data
y = boston.target

# 数据划分
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# 特征标准化（对线性模型很重要，对树模型无所谓但做了也无害）
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# 初始化模型
models = {
    'Linear Regression': LinearRegression(),
    'Ridge Regression (alpha=1.0)': Ridge(alpha=1.0),
    'Decision Tree (max_depth=5)': DecisionTreeRegressor(max_depth=5, random_state=42),
    'Random Forest (n_estimators=100)': RandomForestRegressor(n_estimators=100, random_state=42)
}

# 训练、预测、评估
results = {}
for name, model in models.items():
    model.fit(X_train_scaled, y_train)
    y_pred = model.predict(X_test_scaled)
    mse = mean_squared_error(y_test, y_pred)
    r2 = r2_score(y_test, y_pred)
    results[name] = {'MSE': mse, 'R2': r2}
    print(f"{name:30} - Test MSE: {mse:.4f}, Test R2: {r2:.4f}")

# 使用交叉验证比较（以Random Forest为例）
rf = RandomForestRegressor(random_state=42)
cv_scores_mse = -cross_val_score(rf, X, y, cv=5, scoring='neg_mean_squared_error')
cv_scores_r2 = cross_val_score(rf, X, y, cv=5, scoring='r2')
print(f"\nRandom Forest 5-Fold CV:")
print(f"  MSE: {cv_scores_mse.mean():.4f} (+/- {cv_scores_mse.std()*2:.4f})")
print(f"  R2:  {cv_scores_r2.mean():.4f} (+/- {cv_scores_r2.std()*2:.4f})")

# 简单的网格搜索调参（以决策树为例）
param_grid = {
    'max_depth': [3, 5, 10, None],
    'min_samples_split': [2, 5, 10],
    'min_samples_leaf': [1, 2, 4]
}
tree = DecisionTreeRegressor(random_state=42)
grid_search = GridSearchCV(tree, param_grid, cv=5, scoring='neg_mean_squared_error', n_jobs=-1)
grid_search.fit(X_train_scaled, y_train)
print(f"\nBest Decision Tree Parameters: {grid_search.best_params_}")
print(f"Best CV MSE: {-grid_search.best_score_:.4f}")

运行这段代码，你会直观地看到不同模型的表现差异。通常，随机森林会表现最好，因为它通过集成降低了方差。线性模型可能因为数据中的非线性关系而表现稍差，但它的可解释性最强。决策树单独使用容易过拟合，表现不稳定。

关键经验 ：没有“最好”的模型，只有“最合适”的模型。选择模型时，必须权衡预测精度、计算效率、模型复杂度和可解释性。在数学建模中，如果问题要求对预测结果给出解释（比如“哪些因素对房价影响最大”），那么线性回归或带特征重要性的随机森林可能比一个黑箱的深度神经网络更合适。