数学建模模型与方法——Python模板

一、数据处理方法

1. 数据标准化

（1）Z - score 标准化

数学公式：

\[X_{\text{norm}} = \frac{X - \mu}{\sigma} \]

其中，\(\mu\) 是均值，\(\sigma\) 是标准差

解释：

将数据转换为均值为 0、标准差为 1 的分布。每个特征的值减去该特征的均值，再除以标准差

（2）最小 - 最大归一化（Min-Max Scaling）

数学公式：

\[X_{\text{norm}} = \frac{X - \min(X)}{\max(X) - \min(X)} \]

结果范围通常为 [0, 1]

解释：

将数据缩放到固定区间（如 [0,1]）。每个特征的值减去该特征的最小值，再除以最大值与最小值的差

（3）小数定标标准化（Decimal Scaling）

数学公式：

\[X_{\text{norm}} = \frac{X}{10^k} \]

其中，\(k\) 满足 \(\max(|X|) < 10^k\) 的最小整数（即 \(k = \lceil \log_{10}(\max(|X|)) \rceil\) ）

解释：

通过移动小数点位置（除以 10 的幂）将数据缩放到 [-1,1] 之间。缩放因子取决于数据的最大绝对值

Code（Python）：

def data_normalization(X, method='z-score'):
    """
    数据标准化
    X: 原始数据矩阵
    method: 标准化方法('z-score','minmax','decimal')
    返回: 标准化后数据
    """
    if method == 'z-score':
        # Z-score标准化
        mean = np.mean(X, axis=0)
        std = np.std(X, axis=0)
        return (X - mean) / (std + 1e-8)
    
    elif method == 'minmax':
        # 最小-最大归一化
        min_vals = np.min(X, axis=0)
        max_vals = np.max(X, axis=0)
        return (X - min_vals) / (max_vals - min_vals + 1e-8)
    
    elif method == 'decimal':
        # 小数定标标准化
        max_abs = np.max(np.abs(X), axis=0)
        k = np.ceil(np.log10(max_abs + 1))
        return X / (10**k)
    
    else:
        raise ValueError("Unsupported normalization method")

2. 缺失值处理

（1）均值填充（Mean Imputation）

数学公式：

对于包含缺失值的特征 \(j\)，用该特征所有非缺失值的均值 \(\bar{x}_j\) 填充：

\(x_{i,j} = \bar{x}_j \quad \text{if } x_{i,j} \text{ is missing}\)

其中，\(\bar{x}_j = \frac{1}{n} \sum_{i=1}^{n} x_{i,j}\)（\(n\) 为非缺失值数量）

解释：

计算每个特征的均值，并用该均值替换该特征中的所有缺失值。

（2）中位数填充（Median Imputation）

数学公式：

对于包含缺失值的特征 \(j\)，用该特征所有非缺失值的中位数 \(\text{median}(x_j)\) 填充：

\(x_{i,j} = \text{median}(x_j) \quad \text{if } x_{i,j} \text{ is missing}\)

解释：

计算每个特征的中位数，并用该中位数替换该特征中的所有缺失值

（3）K 近邻填充（KNN Imputation）

数学公式：

对于缺失值 \(x_{i,j}\)，找到与其最相似的 \(k\) 个样本（基于其他特征的距离），用这 \(k\) 个样本在特征 \(j\) 上的取值的均值 / 中位数填充：

\[x_{i,j} = \frac{1}{k} \sum_{n \in \mathcal{N}_k(i)} x_{n,j} \]

其中，\({N}_k(i)\) 表示与样本 \(i\) 最近的 \(k\) 个样本的索引集合

解释：

基于其他特征的相似度（如欧氏距离）找到最近的 k 个样本，用它们的值来估计缺失值

（4）删除缺失值（Listwise Deletion）

数学公式：

直接移除包含任何缺失值的样本： \(X_{\text{clean}} = \{ x_i \mid x_i \text{ has no missing values} \}\)

解释：

删除所有包含缺失值的行，仅保留完整样本

Code（Python）：

def handle_missing_data(X, method='mean'):
    """
    缺失值处理
    X: 原始数据矩阵(含NaN)
    method: 处理方法('mean','median','knn','drop')
    返回: 处理后的数据
    """
    if method == 'mean':
        # 均值填充
        col_means = np.nanmean(X, axis=0)
        nan_indices = np.where(np.isnan(X))
        X_filled = X.copy()
        for i, j in zip(*nan_indices):
            X_filled[i, j] = col_means[j]
        return X_filled
    
    elif method == 'median':
        # 中位数填充
        col_medians = np.nanmedian(X, axis=0)
        nan_indices = np.where(np.isnan(X))
        X_filled = X.copy()
        for i, j in zip(*nan_indices):
            X_filled[i, j] = col_medians[j]
        return X_filled
    
    elif method == 'knn':
        # K近邻填充
        from sklearn.impute import KNNImputer
        imputer = KNNImputer(n_neighbors=5)
        return imputer.fit_transform(X)
    
    elif method == 'drop':
        # 删除含缺失值的行
        return X[~np.isnan(X).any(axis=1)]
    
    else:
        raise ValueError("Unsupported missing data method")

3. 异常值处理

（1）Z-score 法（基于正态分布）

数学公式：

Z-score（标准分数）计算： \(Z = \frac{X - \mu}{\sigma}\) 其中，\(\mu\) 是均值，\(\sigma\) 是标准差。

异常值判定条件： \(|Z| > \text{threshold}\)（默认阈值为 3）

解释：

该方法基于数据的分布特性，将每个数据点转换为标准分数（Z-score），通过设定阈值识别远离均值的异常值，并提供删除或替换两种处理方式。

Code（Python）：

import numpy as np
import pandas as pd

def z_score_outlier_handle(data, threshold=3, method='drop', replace_value=None):
    """
    基于Z-score的异常值处理
    :param data: 输入数据（numpy数组或DataFrame，支持多特征）
    :param threshold: Z分数阈值（默认3）
    :param method: 处理方式：'drop'删除异常值，'replace'替换异常值
    :param replace_value: 替换值（默认用该特征的中位数）
    :return: 处理后的数据（与输入格式一致）
    """
    # 统一转换为numpy数组处理
    if isinstance(data, pd.DataFrame):
        data_np = data.values
        is_df = True
        cols = data.columns
    else:
        data_np = data.copy()
        is_df = False
    
    # 计算每个特征的Z分数：Z = (x - μ) / σ
    mean = np.mean(data_np, axis=0)
    std = np.std(data_np, axis=0)
    std[std == 0] = 1e-8  # 避免除以0
    z_scores = (data_np - mean) / std
    
    # 标记异常值（任何特征的Z分数超出阈值）
    outliers = np.any(np.abs(z_scores) > threshold, axis=1)
    
    # 处理异常值
    if method == 'drop':
        processed = data_np[~outliers]
    elif method == 'replace':
        processed = data_np.copy()
        for col in range(processed.shape[1]):
            # 替换当前特征的异常值
            col_outliers = np.abs(z_scores[:, col]) > threshold
            replace_val = np.median(processed[:, col]) if replace_value is None else replace_value
            processed[col_outliers, col] = replace_val
    else:
        raise ValueError("method必须为'drop'或'replace'")
    
    # 转回DataFrame（如果输入是DataFrame）
    return pd.DataFrame(processed, columns=cols) if is_df else processed

（2） IQR 法（四分位距法）

数学公式：

1. 四分位数计算：
   • Q1：第 25 百分位数（下四分位数）
   • Q3：第 75 百分位数（上四分位数）
   • IQR（四分位距）：IQR = Q3 - Q1
2. 异常值判定边界：
   • 下界：lower = Q1 - 1.5 × IQR
   • 上界：upper = Q3 + 1.5 × IQR
   • 异常值条件：x <lower 或 x> upper

解释：

该方法基于数据的四分位数范围识别异常值。IQR 表示数据中间 50% 的范围，通过设定上下界（通常为 IQR 的 1.5 倍），将远离中间区域的数据点视为异常。相比 Z-score 方法，IQR 对数据分布没有严格要求，更适合处理非正态分布数据。

Code（Python）：

import numpy as np
import pandas as pd

def iqr_outlier_handle(data, method='drop', replace_value=None):
    """
    基于四分位距（IQR）的异常值处理
    :param data: 输入数据（numpy数组或DataFrame）
    :param method: 处理方式：'drop'删除，'replace'替换
    :param replace_value: 替换值（默认用中位数）
    :return: 处理后的数据
    """
    if isinstance(data, pd.DataFrame):
        data_np = data.values
        is_df = True
        cols = data.columns
    else:
        data_np = data.copy()
        is_df = False
    
    # 计算四分位数：Q1（25%）、Q3（75%）、IQR = Q3 - Q1
    q1 = np.percentile(data_np, 25, axis=0)
    q3 = np.percentile(data_np, 75, axis=0)
    iqr = q3 - q1
    
    # 异常值边界：[Q1 - 1.5*IQR, Q3 + 1.5*IQR]（1.5为经验系数）
    lower = q1 - 1.5 * iqr
    upper = q3 + 1.5 * iqr
    
    # 标记异常值（超出边界）
    outliers = np.any((data_np < lower) | (data_np > upper), axis=1)
    
    # 处理异常值
    if method == 'drop':
        processed = data_np[~outliers]
    elif method == 'replace':
        processed = data_np.copy()
        for col in range(processed.shape[1]):
            col_outliers = (processed[:, col] < lower[col]) | (processed[:, col] > upper[col])
            replace_val = np.median(processed[:, col]) if replace_value is None else replace_value
            processed[col_outliers, col] = replace_val
    else:
        raise ValueError("method必须为'drop'或'replace'")
    
    return pd.DataFrame(processed, columns=cols) if is_df else processed

（3）DBSCAN（基于密度的聚类异常检测）

方法核心功能

该函数通过 DBSCAN（密度聚类算法）识别多变量数据中的异常值，并提供 “删除” 或 “替换” 两种处理方式。其核心逻辑是：利用 DBSCAN 对数据进行聚类，将聚类结果中标签为-1的点判定为异常值（非核心点且不被任何核心点 “密度可达”），再通过删除异常值或用中位数（或自定义值）替换异常值，实现数据清洗。

数学基础（DBSCAN 核心概念与公式）

DBSCAN 的异常值判定基于 “密度” 定义，核心概念包括：

1. ε- 邻域（ε-neighborhood） 对于数据点x，其 ε- 邻域是所有与x的距离≤ε 的点的集合，公式为： \(N_\varepsilon(x) = \{ y \mid \text{distance}(x, y) \leq \varepsilon \}\) （注：默认使用欧氏距离，即 \(\text{distance}(x,y) = \sqrt{\sum_{i=1}^d (x_i - y_i)^2}\)，其中d为特征维度）
2. 核心点（Core Point） 若一个点x的 ε- 邻域内至少包含min_samples个点（含自身），则x为核心点： \(|N_\varepsilon(x)| \geq \text{min\_samples}\)
3. 密度可达（Density-Reachable） 若存在核心点序列 \(x_1, x_2, ..., x_k\)，满足 \(x_1=x\)， \(x_k=y\)，且对任意 \(i\) ，\(x_{i+1} \in N_\varepsilon(x_i)\)，则称y从x密度可达。
4. 异常值（Outlier） 既不是核心点，也不被任何核心点 “密度可达” 的点，即为异常值。在 DBSCAN 中，这类点会被标记为-1。

Code（Python）：

import numpy as np
import pandas as pd
from sklearn.cluster import DBSCAN

def dbscan_outlier_handle(data, eps=0.5, min_samples=5, method='drop', replace_value=None):
    """
    基于DBSCAN的异常值处理（适用于多变量）
    :param data: 输入数据（numpy数组或DataFrame）
    :param eps: DBSCAN的邻域半径
    :param min_samples: 形成密集区域的最小样本数
    :param method: 处理方式：'drop'删除，'replace'替换
    :return: 处理后的数据
    """
    if isinstance(data, pd.DataFrame):
        data_np = data.values
        is_df = True
        cols = data.columns
    else:
        data_np = data.copy()
        is_df = False
    
    # DBSCAN聚类：标签为-1的是异常值
    dbscan = DBSCAN(eps=eps, min_samples=min_samples)
    labels = dbscan.fit_predict(data_np)
    outliers = labels == -1  # 标记异常值
    
    # 处理异常值
    if method == 'drop':
        processed = data_np[~outliers]
    elif method == 'replace':
        processed = data_np.copy()
        for col in range(processed.shape[1]):
            replace_val = np.median(processed[:, col]) if replace_value is None else replace_value
            processed[outliers, col] = replace_val
    else:
        raise ValueError("method必须为'drop'或'replace'")
    
    return pd.DataFrame(processed, columns=cols) if is_df else processed

（4）孤立森林（Isolation Forest）

核心原理：

孤立森林是一种基于树结构的无监督异常检测算法，其核心思想是：异常值与正常样本相比，具有 “更容易被孤立” 的特性。

• 正常样本往往在数据空间中形成密集的聚类，需要更多次划分才能被分离；
• 异常样本由于偏离整体分布，通常可以通过更少的随机划分被孤立出来。

算法通过构建多棵随机决策树（“孤立树”），计算每个样本在这些树中的 “平均路径长度”（即从根节点到叶子节点所需的划分次数），路径越短的样本越可能是异常值。

数学公式：

孤立森林的异常判断依赖于异常分数（Anomaly Score） 的计算，其核心公式如下：

a. 路径长度（Path Length）

对于样本 \(x\)，在一棵孤立树中，其路径长度 \(h(x)\) 定义为：从树的根节点到 \(x\) 所在叶子节点经过的边数（划分次数）。

• 若节点为叶子节点（无法再划分），则 \(h(x) = 0\)；
• 否则，\(h(x) = 1 + h(\text{child})\)，其中 \(\text{child}\) 是 \(x\) 进入的子节点（左 / 右子树）。

b. 平均路径长度的期望值（Normalization）

为了消除树结构差异的影响，需对路径长度进行标准化。对于包含 \(n\) 个样本的数据集，单棵树中平均路径长度的期望值 \(c(n)\) 定义为：

\(c(n) = \begin{cases} 0 & \text{if } n \leq 1, \\ 2(\ln(n-1) + \gamma) - \frac{2(n-1)}{n} & \text{if } n \geq 2, \end{cases}\)

其中 \(\gamma \approx 0.5772\) 是欧拉常数（Euler-Mascheroni constant），用于修正离散样本的偏差。

c. 异常分数（Anomaly Score）

对于样本 \(x\)，在由 t 棵孤立树组成的森林中，其异常分数 \(s(x, n)\) 定义为：\(s(x, n) = 2^{-\frac{\mathbb{E}[h(x)]}{c(n)}}\) 其中 \(\mathbb{E}[h(x)]\) 是 \(x\) 在所有树中的平均路径长度。

• 若 \(s(x) \approx 1\)：样本 \(x\) 的平均路径长度远小于正常样本，判定为异常值；
• 若 \(s(x) \approx 0.5\)：样本 \(x\) 的路径长度接近正常样本的平均水平，判定为正常样本；
• 若 \(s(x) \approx 0\)：样本 \(x\) 路径长度较长，属于密集分布的正常样本。

Code（Python）：

import numpy as np
import pandas as pd
from sklearn.ensemble import IsolationForest

def isolation_forest_outlier_handle(data, n_estimators=100, contamination='auto', method='drop', replace_value=None):
    """
    基于孤立森林的异常值处理
    :param data: 输入数据（numpy数组或DataFrame）
    :param n_estimators: 树的数量
    :param contamination: 异常值比例（'auto'或0~0.5的浮点数）
    :param method: 处理方式：'drop'删除，'replace'替换
    :return: 处理后的数据
    """
    if isinstance(data, pd.DataFrame):
        data_np = data.values
        is_df = True
        cols = data.columns
    else:
        data_np = data.copy()
        is_df = False
    
    # 孤立森林：标签为-1的是异常值
    iso_forest = IsolationForest(n_estimators=n_estimators, contamination=contamination, random_state=42)
    labels = iso_forest.fit_predict(data_np)
    outliers = labels == -1
    
    # 处理异常值
    if method == 'drop':
        processed = data_np[~outliers]
    elif method == 'replace':
        processed = data_np.copy()
        for col in range(processed.shape[1]):
            replace_val = np.median(processed[:, col]) if replace_value is None else replace_value
            processed[outliers, col] = replace_val
    else:
        raise ValueError("method必须为'drop'或'replace'")
    
    return pd.DataFrame(processed, columns=cols) if is_df else processed

二、优化算法

1. 线性规划

数学模型

线性规划（Linear Programming, LP）是一种求解线性目标函数在线性约束条件下最优解的数学方法。其标准形式（以最小化目标为例）如下：

目标函数：

\(\min_{\mathbf{x}} \quad \mathbf{c}^T \mathbf{x}\)

约束条件：

\[\begin{cases} \mathbf{A}_{ub} \mathbf{x} \leq \mathbf{b}_{ub} \quad \text{（不等式约束）} \\ \mathbf{A}_{eq} \mathbf{x} = \mathbf{b}_{eq} \quad \text{（等式约束）} \\ l_i \leq x_i \leq u_i \quad \text{（变量边界约束，由bounds定义）} \end{cases} \]

其中，\(\mathbf{x} = (x_1, x_2, ..., x_n)^T\) 为决策变量向量；\(\mathbf{c} = (c_1, c_2, ..., c_n)^T\) 为目标函数系数向量；\(\mathbf{A}_{ub}\) 为不等式约束系数矩阵，\(\mathbf{b}_{ub}\) 为不等式约束右端项向量；\(\mathbf{A}_{eq}\) 为等式约束系数矩阵，\(\mathbf{b}_{eq}\) 为等式约束右端项向量；\(l_i, u_i\) 为变量 \(x_i\) 的下界和上界。

Code（Python）：

from scipy.optimize import linprog

def linear_programming(c, A_ub, b_ub, A_eq, b_eq, bounds):
    """
    线性规划求解
    c: 目标函数系数
    A_ub, b_ub: 不等式约束
    A_eq, b_eq: 等式约束
    bounds: 变量边界
    返回: 最优解、最优值
    """
    res = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=bounds)
    return res.x, res.fun

2. 遗传算法

方法概述：

遗传算法（Genetic Algorithm, GA），通过模拟生物进化过程（选择、交叉、变异），在连续搜索空间中寻找目标函数的最优解（最小化）。核心逻辑是：初始化随机种群，通过迭代保留优质个体（精英保留）、选择父代（锦标赛选择）、交叉产生子代（SBX 交叉）、引入变异（多项式变异），逐步逼近最优解。

核心步骤及数学公式

1. 种群初始化

• 目标：生成初始解的随机集合，确保解空间覆盖的均匀性。
• 数学公式： 对于每个变量 \(x_i\)（维度 \(i=1,2,\dots,n\)），从均匀分布中采样： \(x_i \sim \text{Uniform}(l_i, u_i)\) 其中 \(l_i\) 和 \(u_i\) 是变量 \(x_i\) 的下界和上界（由 bounds 参数定义）。

2. 适应度评估

• 目标：量化每个解的优劣程度（此处为最小化目标函数）。
• 数学公式： 适应度函数直接使用目标函数值： \(\text{Fitness}(x) = f(x)\) 其中 \(f(x)\) 是待优化的函数（如最小化 \(x^2\)）。

3. 选择操作（锦标赛选择）

• 目标：基于适应度概率性地选择父代，保留优质解。
• 数学公式： 从种群中随机选择 \(k\) 个个体（代码中 \(k=3\)），选择适应度最优的个体： \(\text{Parent} = \arg\min_{x \in \text{Candidates}} \text{Fitness}(x)\) 其中 \(\text{Candidates} = \{x_1, x_2, x_3\}\) 是随机抽取的候选集。

4. 交叉操作（模拟二进制交叉，SBX）

• 目标：通过父代基因交换生成新解，模拟生物遗传。
• 数学公式： 对于父代 \(x_1\) 和 \(x_2\)，生成子代 \(c_1\) 和 \(c_2\) 的步骤如下：
  1. 计算参数 \(\beta\)（控制交叉强度）： \(\beta = 1 + \frac{2 \cdot \min(x_1 - l_i, u_i - x_2)}{|x_1 - x_2|}\) 其中 \(l_i\) 和 \(u_i\) 是变量边界，确保子代不超出边界。
  2. **计算参数 \(\alpha\)： \(\alpha = 2 - \beta^{-(\eta_c + 1)}\) 其中 \(\eta_c\) 是分布指数，值越大表示子代越接近父代。
  3. 生成随机因子 \(\beta_q\)： \(\beta_q = \begin{cases} (u \cdot \alpha)^{1/(\eta_c + 1)} & \text{若 } u \leq 1/\alpha \\ \left(\frac{1}{2 - u \cdot \alpha}\right)^{1/(\eta_c + 1)} & \text{否则} \end{cases}\) 其中 \(u \sim \text{Uniform}(0,1)\)。
  4. 计算子代： \(c_1 = 0.5 \left[(x_1 + x_2) - \beta_q \cdot |x_1 - x_2|\right],c_2 = 0.5 \left[(x_1 + x_2) + \beta_q \cdot |x_1 - x_2|\right]\)
  5. 边界修正： \(c_1 = \max(l_i, \min(c_1, u_i)) ,c_2 = \max(l_i, \min(c_2, u_i))\)

5. 变异操作（多项式变异）

• 目标：引入随机扰动，避免算法陷入局部最优。
• 数学公式： 对于变量 \(x\)，变异后的值 \(x'\) 的计算步骤：
  1. 计算相对位置： \(\delta_1 = \frac{x - l_i}{u_i - l_i}, \quad \delta_2 = \frac{u_i - x}{u_i - l_i}\)
  2. 生成变异因子 \(\delta_q\)： \(\delta_q = \begin{cases} \left(2u + (1-2u)(1-\delta_1)^{\eta_m + 1}\right)^{1/(\eta_m + 1)} - 1 & \text{若 } u \leq 0.5 \\ 1 - \left(2(1-u) + 2(u-0.5)(1-\delta_2)^{\eta_m + 1}\right)^{1/(\eta_m + 1)} & \text{否则} \end{cases}\) 其中 \(\eta_m\) 是分布指数，控制变异幅度。
  3. 更新变量： \(x' = x + \delta_q \cdot (u_i - l_i)\)
  4. 边界修正： \(x' = \max(l_i, \min(x', u_i))\)

6. 精英保留策略

• 目标：确保每代中最优解不被破坏。
• 数学公式： 直接保留适应度前 \(k\) 个个体： \(\text{Elite} = \left\{x \mid \text{Fitness}(x) \leq \text{Fitness}(x') \text{ for all } x' \in \text{Population}\right\}\)

Code（Python）：

import numpy as np
import matplotlib.pyplot as plt

def genetic_algorithm(objective_func, bounds, pop_size=50, elite_size=2, 
                      mutation_rate=0.1, crossover_rate=0.9, generations=100,
                      verbose=True, plot_progress=False):
    """
    遗传算法优化函数
    
    参数:
    objective_func -- 目标函数（需最小化的函数）
    bounds -- 变量边界列表，例如[(min1, max1), (min2, max2), ...]
    pop_size -- 种群大小 (默认50)
    elite_size -- 直接保留的精英个体数量 (默认2)
    mutation_rate -- 变异概率 (默认0.1)
    crossover_rate -- 交叉概率 (默认0.9)
    generations -- 迭代次数 (默认100)
    verbose -- 是否打印进度 (默认True)
    plot_progress -- 是否绘制收敛曲线 (默认False)
    
    返回:
    best_solution -- 找到的最优解
    best_fitness -- 最优解对应的适应度值
    """
    
    # 1. 初始化参数
    n_var = len(bounds)  # 变量维度
    lower_bounds = np.array([b[0] for b in bounds])
    upper_bounds = np.array([b[1] for b in bounds])
    
    # 2. 初始化种群
    population = np.random.uniform(low=lower_bounds, 
                                  high=upper_bounds, 
                                  size=(pop_size, n_var))
    
    # 3. 存储最优解历史
    best_fitness_history = []
    
    # 4. 主循环
    for gen in range(generations):
        # 计算适应度
        fitness = np.array([objective_func(ind) for ind in population])
        
        # 记录最佳解
        best_idx = np.argmin(fitness)
        best_fitness = fitness[best_idx]
        best_solution = population[best_idx].copy()
        best_fitness_history.append(best_fitness)
        
        if verbose and gen % 10 == 0:
            print(f"Generation {gen}: Best Fitness = {best_fitness:.6f}")
        
        # 选择精英
        elite_indices = np.argsort(fitness)[:elite_size]
        elite_population = population[elite_indices].copy()
        
        # 初始化下一代种群
        new_population = np.empty((pop_size, n_var))
        new_population[:elite_size] = elite_population
        
        # 锦标赛选择
        for i in range(elite_size, pop_size):
            # 随机选择3个候选个体
            candidates = np.random.choice(pop_size, size=3, replace=False)
            winner_idx = candidates[np.argmin(fitness[candidates])]
            new_population[i] = population[winner_idx]
        
        # 交叉操作 (SBX)
        for i in range(elite_size, pop_size, 2):
            if np.random.rand() < crossover_rate and i+1 < pop_size:
                parent1, parent2 = new_population[i], new_population[i+1]
                child1, child2 = sbx_crossover(parent1, parent2, lower_bounds, upper_bounds)
                new_population[i] = child1
                new_population[i+1] = child2
        
        # 变异操作 (多项式变异)
        for i in range(elite_size, pop_size):
            if np.random.rand() < mutation_rate:
                new_population[i] = polynomial_mutation(new_population[i], 
                                                       lower_bounds, 
                                                       upper_bounds)
        
        # 更新种群
        population = new_population
    
    # 最终评估
    final_fitness = np.array([objective_func(ind) for ind in population])
    best_idx = np.argmin(final_fitness)
    best_solution = population[best_idx]
    best_fitness = final_fitness[best_idx]
    
    if verbose:
        print(f"\nFinal Generation: Best Fitness = {best_fitness:.6f}")
        print(f"Optimal Solution: {best_solution}")
    
    if plot_progress:
        plt.plot(best_fitness_history)
        plt.title('Convergence Curve')
        plt.xlabel('Generation')
        plt.ylabel('Best Fitness')
        plt.grid(True)
        plt.show()
    
    return best_solution, best_fitness

def sbx_crossover(parent1, parent2, lower_bounds, upper_bounds, eta_c=20):
    """模拟二进制交叉 (SBX)"""
    child1, child2 = parent1.copy(), parent2.copy()
    for i in range(len(parent1)):
        if np.random.rand() <= 0.5:
            continue
            
        x1, x2 = parent1[i], parent2[i]
        xl, xu = lower_bounds[i], upper_bounds[i]
        
        if abs(x1 - x2) > 1e-10:
            if x2 > x1:
                x1, x2 = x2, x1
                
            beta = 1.0 + (2.0 * (x1 - xu) / (x2 - x1))
            alpha = 2.0 - beta**(-(eta_c + 1.0))
            u = np.random.rand()
            
            if u <= 1.0/alpha:
                beta_q = (u * alpha)**(1.0/(eta_c + 1.0))
            else:
                beta_q = (1.0/(2.0 - u*alpha))**(1.0/(eta_c + 1.0))
                
            c1 = 0.5 * (x1 + x2 - beta_q * (x1 - x2))
            
            beta = 1.0 + (2.0 * (xl - x2) / (x2 - x1))
            alpha = 2.0 - beta**(-(eta_c + 1.0))
            
            if u <= 1.0/alpha:
                beta_q = (u * alpha)**(1.0/(eta_c + 1.0))
            else:
                beta_q = (1.0/(2.0 - u*alpha))**(1.0/(eta_c + 1.0))
                
            c2 = 0.5 * (x1 + x2 + beta_q * (x1 - x2))
            
            c1 = min(max(c1, xl), xu)
            c2 = min(max(c2, xl), xu)
            
            if np.random.rand() <= 0.5:
                child1[i] = c2
                child2[i] = c1
            else:
                child1[i] = c1
                child2[i] = c2
                
    return child1, child2

def polynomial_mutation(individual, lower_bounds, upper_bounds, eta_m=20):
    """多项式变异"""
    mutated = individual.copy()
    for i in range(len(individual)):
        if np.random.rand() < 1.0/len(individual):  # 每个基因变异概率
            x = individual[i]
            xl, xu = lower_bounds[i], upper_bounds[i]
            
            delta1 = (x - xl) / (xu - xl)
            delta2 = (xu - x) / (xu - xl)
            u = np.random.rand()
            
            if u <= 0.5:
                delta_q = (2*u + (1-2*u)*(1-delta1)**(eta_m+1))**(1/(eta_m+1)) - 1
            else:
                delta_q = 1 - (2*(1-u) + 2*(u-0.5)*(1-delta2)**(eta_m+1))**(1/(eta_m+1))
                
            x = x + delta_q * (xu - xl)
            mutated[i] = min(max(x, xl), xu)
            
    return mutated

3. 模拟退火算法

1. 核心原理

模拟退火算法（Simulated Annealing, SA）灵感来源于物理中金属退火的过程：高温时金属原子可自由运动（允许较大范围探索），随温度降低，原子逐渐稳定在低能量状态（收敛到最优解）。算法通过概率性接受 “较差解” 避免陷入局部最优，最终在低温时收敛到全局最优。

2. 数学公式与核心步骤

代码实现的核心步骤对应以下数学逻辑：

（1）新解生成

每次迭代通过加入与当前温度相关的随机扰动生成新解： \(\text{new\_sol} = \text{current\_sol} + \mathcal{N}(0, T^2) \tag{1}\) 其中，\(\mathcal{N}(0, T^2)\) 是均值为 \(0\)、方差为 \(T^2\) 的正态分布噪声（温度 \(T\) 越高，扰动越大，探索范围越广）。

（2）接受准则（Metropolis 准则）

• 若新解代价更低（\(\text{new\_cost} < \text{current\_cost}\)），直接接受： \(\text{current\_sol} \leftarrow \text{new\_sol} \quad (\text{if } \Delta \text{cost} < 0)\)
• 若新解代价更高（\(\text{new\_cost} \geq \text{current\_cost}\)，以概率接受： \(P(\text{accept}) = \exp\left(-\frac{\Delta \text{cost}}{T}\right) \tag{2}\) 其中 \(\Delta \text{cost} = \text{new\_cost} - \text{current\_cost}\)（代价差），\(T\) 为当前温度。

（3）降温策略

温度随迭代逐步降低，控制探索向利用过渡： \(T_{t+1} = \alpha \cdot T_t \tag{3}\) 其中 \(\alpha\) 为冷却系数（\(0 < \alpha < 1\)，代码中 \(\alpha=0.9\)），温度越低，接受较差解的概率越小。

Code（Python）：

import numpy as np

def simulated_annealing(cost_func, initial_solution, T=1.0, T_min=1e-5, alpha=0.9, max_iter=1000):
    """
    模拟退火优化算法
    cost_func: 代价函数
    initial_solution: 初始解
    T: 初始温度
    T_min: 最小温度
    alpha: 冷却系数
    max_iter: 最大迭代次数
    返回: 最优解、最优值
    """
    current_sol = initial_solution
    best_sol = current_sol
    current_cost = cost_func(current_sol)
    best_cost = current_cost
    
    while T > T_min and max_iter > 0:
        # 生成新解
        new_sol = current_sol + np.random.randn(*current_sol.shape) * T
        new_cost = cost_func(new_sol)
        
        # 接受更优解或以概率接受劣解
        delta_cost = new_cost - current_cost
        if delta_cost < 0 or np.random.rand() < np.exp(-delta_cost / T):
            current_sol = new_sol
            current_cost = new_cost
            
            # 更新最优解
            if current_cost < best_cost:
                best_sol = current_sol
                best_cost = current_cost
        
        # 降温
        T *= alpha
        max_iter -= 1
    
    return best_sol, best_cost

三、预测模型

1. 线性回归模型

数学原理

线性回归是最简单的监督学习模型，假设目标变量与特征之间存在线性关系，通过最小化预测值与真实值的平方误差拟合模型。

模型表达式：

对于特征矩阵 \(X \in {R}^{m \times n}\)（\(m\) 为样本数，\(n\) 为特征数）和目标向量 \(y \in {R}^m\)，线性回归模型为： \(\hat{y} = Xw + b \tag{1}\) 其中，\(w \in \mathbb{R}^n\) 是特征系数（权重），\(b \in {R}\)) 是截距，\(\hat{y}\) 是预测值。

目标函数（损失函数）： 通过最小化残差平方和（RSS）求解 \(w\) 和 \(b\)： \(\min_{w, b} \quad \text{RSS} = \sum_{i=1}^m (\hat{y}_i - y_i)^2 = \|y - Xw - b\|_2^2 \tag{2}\)

求解结果： 系数 \(w\) 和截距 \(b\) 可通过最小二乘法解析求解（要求 \(X^TX\) 可逆）

Code（Python）：

from sklearn.linear_model import LinearRegression

def linear_regression(X, y):
    """
    线性回归模型
    X: 特征矩阵(m×n)
    y: 目标向量(m)
    返回: 模型对象、系数、截距
    """
    model = LinearRegression()
    model.fit(X, y)
    coef = model.coef_
    intercept = model.intercept_
    return model, coef, intercept

2. 灰色预测GM(1,1)

核心原理

灰色预测是针对 “少数据、信息不完全” 系统的预测方法，GM (1,1)（一阶单变量灰色模型）是最常用的类型。其核心思想是通过对原始数据进行 “累加生成”（AGO）弱化随机性，转化为近似指数增长的有规律序列，再通过微分方程拟合规律，最后还原得到预测结果。

数学公式

1. 原始序列与累加生成（AGO） 设原始非负序列为 \(x^{(0)} = [x^{(0)}(1), x^{(0)}(2), ..., x^{(0)}(n)]\)，累加生成序列 \(x^{(1)}\) 定义为：\(x^{(1)}(k) = \sum_{i=1}^k x^{(0)}(i) \quad (k=1,2,...,n)\) 累加后，\(x^{(1)}\) 通常呈现近似指数增长趋势，随机性降低。
2. 微分方程模型 对 \(x^{(1)}\) 构建一阶线性微分方程（白化方程）：\(\frac{dx^{(1)}}{dt} + a x^{(1)} = b\) 其中 \(a\) 为发展系数（反映增长趋势），\(b\) 为灰作用量（反映外部驱动）。
3. 参数求解 通过最小二乘法估计参数 \(\hat{a}, \hat{b}\)。构造矩阵 \(B\) 和向量 \(Y\)：\(B = \begin{bmatrix} -0.5(x^{(1)}(1)+x^{(1)}(2)) & 1 \\ -0.5(x^{(1)}(2)+x^{(1)}(3)) & 1 \\ ... & ... \\ -0.5(x^{(1)}(n-1)+x^{(1)}(n)) & 1 \end{bmatrix}, \quad Y = \begin{bmatrix} x^{(0)}(2) \\ x^{(0)}(3) \\ ... \\ x^{(0)}(n) \end{bmatrix}\) 参数估计为：\(\begin{bmatrix} \hat{a} \\ \hat{b} \end{bmatrix} = (B^T B)^{-1} B^T Y\)
4. 预测与还原 微分方程的解（拟合 / 预测公式）为：\(\hat{x}^{(1)}(k+1) = \left( x^{(0)}(1) - \frac{\hat{b}}{\hat{a}} \right) e^{-\hat{a}k} + \frac{\hat{b}}{\hat{a}}\) 通过累减还原为原始序列预测值：\(\hat{x}^{(0)}(k+1) = \hat{x}^{(1)}(k+1) - \hat{x}^{(1)}(k)\)

Code（Python）：

def grey_model(data, predict_num=1):
    """
    灰色预测GM(1,1)模型
    data: 原始序列(1维数组)
    predict_num: 预测步数
    返回: 预测序列、模型参数[a,b]
    """
    n = len(data)
    # 累加生成
    AGO = np.cumsum(data)
    
    # 构造矩阵B,Y
    B = np.zeros((n-1, 2))
    Y = data[1:].reshape(-1, 1)
    for i in range(n-1):
        B[i, 0] = -0.5*(AGO[i] + AGO[i+1])
        B[i, 1] = 1
    
    # 求解参数
    a, b = np.linalg.inv(B.T @ B) @ B.T @ Y
    a, b = float(a), float(b)
    
    # 预测
    predict = []
    for k in range(n + predict_num):
        if k < n:  # 拟合值
            predict.append((data[0] - b/a) * np.exp(-a*k) + b/a)
        else:  # 预测值
            predict.append((data[0] - b/a) * np.exp(-a*k) + b/a)
    
    # 还原结果
    predict_restore = [predict[0]]
    for i in range(1, len(predict)):
        predict_restore.append(predict[i] - predict[i-1])
    
    return predict_restore, (a, b)

3.ARIMA 时间序列预测

核心原理

ARIMA（自回归积分移动平均模型）是处理非平稳时间序列的经典方法，通过 “差分”（消除趋势使序列平稳）、“自回归”（利用历史值预测）、“移动平均”（消除随机波动）结合，捕捉时间序列的动态规律。其参数表示为 \(\text{ARIMA}(p, d, q)\)，其中：

• \(p\)：自回归项数（使用过去 \(p\) 个值预测当前值）；
• \(d\)：差分次数（使序列平稳的差分阶数）；
• \(q\)：移动平均项数（使用过去 \(q\) 个误差项修正预测）。

数学公式

1. 差分操作：对非平稳序列 \(x_t\) 进行 \(d\) 次差分，得到平稳序列 \(y_t = \nabla^d x_t\)（\(\nabla\) 为差分算子，\(\nabla x_t = x_t - x_{t-1}\)）。
2. ARIMA 模型表达式： 对平稳序列 \(y_t\)，模型为：\(y_t = c + \sum_{i=1}^p \phi_i y_{t-i} + \sum_{j=1}^q \theta_j \varepsilon_{t-j} + \varepsilon_t\) 其中：
   • \(c\) 为常数项；
   • \(\phi_i\) 为自回归系数（反映历史值对当前值的影响）；
   • \(\theta_j\) 为移动平均系数（反映历史误差对当前值的修正）；
   • \(\varepsilon_t\) 为白噪声（均值为 \(0\)、方差恒定的随机误差）。

Code（Python）：

from statsmodels.tsa.arima.model import ARIMA

def arima_model(data, order=(1,0,1), predict_num=5):
    """
    ARIMA时间序列预测
    data: 时间序列数据(1维数组)
    order: (p,d,q)参数
    predict_num: 预测点数
    返回: 预测结果、模型对象z
    """
    model = ARIMA(data, order=order)
    model_fit = model.fit()
    forecast = model_fit.forecast(steps=predict_num)
    return forecast, model_fit

四、分类与聚类

1. 逻辑回归（Logistic Regression）分类

核心原理

逻辑回归是一种用于分类任务的线性模型，尽管名字含 “回归”，但其本质是通过 Sigmoid 函数将线性回归的输出映射到 [0,1] 区间，作为样本属于某一类别的概率，最终通过阈值判断类别（二分类或多分类）。核心是学习特征与类别概率之间的线性关系。

数学公式

1. 线性得分计算： 对输入特征 \(X\)（样本数\(m \times\)特征数\(n\)），先计算线性得分： \(z = w \cdot X + b \tag{1}\) 其中，\(w \in {R}^n\) 是特征系数，\(b \in {R}\) 是偏置项。

2. Sigmoid 函数（概率映射）： 通过 Sigmoid 函数将线性得分 \(z\) 映射为类别概率（二分类场景，输出为样本属于 “正类” 的概率）： \(\hat{y} = \sigma(z) = \frac{1}{1 + e^{-z}} \tag{2}\)

   Sigmoid 函数的特性：当 \(z \to +\infty\) 时，\(\hat{y} \to 1\)；当 \(z \to -\infty\) 时，\(\hat{y} \to 0\)；当 \(z=0\) 时，\(\hat{y}=0.5\)。

3. 决策边界： 设定阈值（通常为 0.5），当 \(\hat{y} \geq 0.5\) 时，预测为正类；否则为负类。对应的线性决策边界为： \(w \cdot X + b = 0 \tag{3}\)

4. 损失函数（交叉熵损失）： 逻辑回归通过最小化交叉熵损失学习参数 \(w\) 和 \(b\)，公式为： \(\mathcal{L}(w, b) = -\frac{1}{m} \sum_{i=1}^m \left[ y_i \log(\hat{y}_i) + (1 - y_i) \log(1 - \hat{y}_i) \right] \tag{4}\) 其中，\(y_i \in \{0,1\}\) 是真实标签，\(\hat{y}_i\) 是预测概率。

Code（Python）：

from sklearn.linear_model import LogisticRegression

def logistic_classification(X, y, C=1.0):
    """
    逻辑回归分类
    C: 正则化强度的倒数（越小正则化越强）
    返回: 模型对象、预测概率
    """
    model = LogisticRegression(C=C, max_iter=1000)
    model.fit(X, y)
    proba = model.predict_proba(X)
    return model, proba

2. K-means 聚类

基本原理

K-means 是一种无监督学习的聚类算法，核心目标是将数据集自动划分为k个互不重叠的簇（cluster），使得簇内样本的相似度尽可能高（簇内距离小），簇间样本的相似度尽可能低（簇间距离大）。其核心流程为：

1. 随机初始化k个簇中心（初始质心）；
2. 计算每个样本到k个中心的距离，将样本分配到最近的簇；
3. 重新计算每个簇的质心（簇内所有样本的均值）；
4. 重复步骤 2-3，直到簇中心不再显著变化（收敛）或达到最大迭代次数。

数学公式

K-means 的目标函数是最小化簇内平方误差和（SSE），即最小化所有样本到其所属簇中心的平方距离之和：

\(J = \sum_{j=1}^{k} \sum_{x_i \in C_j} \| x_i - \mu_j \|^2\)

其中：

• \(k\) 为预设的簇数量；
• \(C_j\) 为第 \(j\) 个簇；
• \(x_i\) 为簇 \(C_j\) 中的样本；
• \(\mu_j\) 为簇 \(C_j\) 的质心（\(\mu_j = \frac{1}{|C_j|} \sum_{x_i \in C_j} x_i\)）；
• \(\| \cdot \|^2\) 为欧氏距离的平方。

Code（Python）：

from sklearn.cluster import KMeans

def kmeans_clustering(X, k=3):
    """
    K-means聚类
    X: 特征矩阵(m×n)
    k: 聚类数量
    返回: 聚类标签、聚类中心
    """
    kmeans = KMeans(n_clusters=k, random_state=0)
    labels = kmeans.fit_predict(X)
    centers = kmeans.cluster_centers_
    return labels, centers

3. SVM 分类（支持向量机）

基本原理

SVM（Support Vector Machine）是一种监督学习的分类算法，核心目标是在特征空间中找到一个最优超平面，使得两类样本的 “间隔”（超平面到两侧最近样本的距离）最大。该超平面由 “支持向量”（距离超平面最近的样本）唯一决定。

• 对于线性可分数据：直接寻找最大化间隔的超平面；
• 对于线性不可分数据：通过核函数（如 RBF 核、多项式核）将低维数据映射到高维空间，使其线性可分，再在高维空间中寻找最优超平面。

数学公式

1. 线性 SVM： 设超平面为\(w \cdot x + b = 0\)（\(w\) 为法向量，\(b\) 为偏置），对样本\((x_i, y_i)\)（\(y_i \in \{+1, -1\}\)为标签），需满足：

   • 正样本：\(y_i=1\)时，\(w \cdot x_i + b \geq 1\)；
   • 负样本：\(y_i=-1\)时，\(w \cdot x_i + b \leq -1\)。

   目标是最小化超平面的 “复杂度”（即\(\frac{1}{2}\|w\|^2\)），同时满足上述约束，即：$ \min_{w,b} \frac{1}{2}|w|^2 \quad \text{s.t.} \quad y_i(w \cdot x_i + b) \geq 1 \quad (\forall i)$

2. 非线性 SVM（核函数）： 引入核函数\(K(x_i, x_j) = \phi(x_i) \cdot \phi(x_j)\)（\(\phi\)为高维映射函数），将低维内积替换为核函数，目标函数变为： \(\min_{\alpha} \frac{1}{2}\sum_{i,j} \alpha_i \alpha_j y_i y_j K(x_i, x_j) - \sum_i \alpha_i \quad \text{s.t.} \quad \sum_i \alpha_i y_i = 0, \alpha_i \geq 0\) （\(\alpha_i\)为拉格朗日乘子，非零\(\alpha_i\)对应的样本即为支持向量）。

Code（Python）：

from sklearn.svm import SVC
from sklearn.model_selection import train_test_split

def svm_classification(X, y, test_size=0.2):
    """
    SVM分类模型
    X: 特征矩阵
    y: 标签向量
    test_size: 测试集比例
    返回: 准确率、模型对象
    """
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=test_size)
    model = SVC(kernel='rbf', C=1.0)
    model.fit(X_train, y_train)
    accuracy = model.score(X_test, y_test)
    return accuracy, model

五、降维方法

1. 主成分分析法（PCA）

定义

主成分分析（Principal Component Analysis, PCA）是一种线性降维方法，通过将高维数据投影到一组正交的 “主成分” 上，在保留数据中最关键信息（方差最大）的同时降低维度。核心思想是：找到数据中方差最大的方向（主成分），用这些方向构成的低维空间近似表示原始高维数据。

数学核心

1. 数据预处理：对原始数据矩阵\(X \in {R}^{m \times n}\)（\(m\)为样本数，\(n\)为特征数）进行中心化（减去均值），得到\(\tilde{X}\)。
2. 协方差矩阵：计算中心化数据的协方差矩阵，描述特征间的相关性：\(C = \frac{1}{m-1} \tilde{X}^T \tilde{X}\)
3. 特征值分解：对协方差矩阵 \(C\) 进行特征值分解，得到特征值 \(\lambda_1 \geq \lambda_2 \geq \dots \geq \lambda_n\)和对应的特征向量\(v_1, v_2, \dots, v_n\)。
   • 特征值 \(\lambda_i\) 表示第 \(i\) 个主成分的方差（信息含量）；
   • 特征向量 \(v_i\) 是第 \(i\) 个主成分的方向（正交向量）。
4. 降维投影：选择前 \(k\) 个最大特征值对应的特征向量 \(v_1, \dots, v_k\) ，构成投影矩阵 \(V_k \in {R}^{n \times k}\)，将高维数据投影到低维空间：\(Y = \tilde{X} V_k\) 其中 \(Y \in {R}^{m \times k}\)为降维后的数据。

Code（Python）：

from sklearn.decomposition import PCA

def pca_reduction(X, n_components=2):
    """
    主成分分析降维
    X: 原始数据矩阵(m×n)
    n_components: 保留主成分数
    返回: 降维后数据、解释方差比
    """
    pca = PCA(n_components=n_components)
    X_reduced = pca.fit_transform(X)
    explained_variance = pca.explained_variance_ratio_
    return X_reduced, explained_variance

2. t-SNE 非线性降维

定义

t 分布邻域嵌入（t-Distributed Stochastic Neighbor Embedding, t-SNE）是一种非线性降维方法，通过构建高维与低维空间中数据点的 “相似性分布”，并最小化两者的差异，重点保留数据的局部邻域结构（适合可视化高维数据的聚类模式）。

数学核心

1. 高维相似性：对高维空间中的点 \(x_i\) 和 \(x_j\)，用高斯分布定义条件概率 \(p_{j|i}\)，表示 \(x_j\) 是 \(x_i\) 邻域的概率：\(p_{j|i} = \frac{\exp\left(-\frac{\|x_i - x_j\|^2}{2\sigma_i^2}\right)}{\sum_{k \neq i} \exp\left(-\frac{\|x_i - x_k\|^2}{2\sigma_i^2}\right)}\) 其中\(\sigma_i\)通过 “困惑度”（perplexity）控制（困惑度≈邻域内的有效样本数，通常取 5~50）。
2. 低维相似性：对低维空间中的对应点 \(y_i\) 和 \(y_j\)，用 \(t\) 分布（自由度为 1）定义对称概率 \(q_{ij}\)，避免高维中 “拥挤问题”：\(q_{ij} = \frac{\left(1 + \|y_i - y_j\|^2\right)^{-1}}{\sum_{k < l} \left(1 + \|y_k - y_l\|^2\right)^{-1}}\)。
3. 优化目标：通过梯度下降最小化高维分布P（\(p_{ij} = \frac{p_{j|i} + p_{i|j}}{2m}\)）与低维分布Q（\(q_{ij}\)）的 KL 散度：\(\text{KL}(P \| Q) = \sum_{i < j} p_{ij} \log\left(\frac{p_{ij}}{q_{ij}}\right)\)。

Code（Python）：

from sklearn.manifold import TSNE

def tsne_visualization(X, n_components=2, perplexity=30):
    """
    t-SNE非线性降维
    X: 原始数据矩阵
    n_components: 目标维度
    perplexity: 困惑度参数
    返回: 降维后数据
    """
    tsne = TSNE(n_components=n_components, perplexity=perplexity)
    X_tsne = tsne.fit_transform(X)
    return X_tsne

六、综合评价方法

1. 熵权法

方法概述 核心思想：通过指标变异程度（信息熵）客观确定权重

特点：数据驱动、避免主观偏差 适用场景：多指标综合评价/决策（如方案优选、绩效评估等）

输入数据要求

1. 原始数据矩阵\({X}_{m×n}\)： - \(m\)个评价对象（行） - \(n\)个评价指标（列）
2. 指标类型标注： - ✅ 效益型（越大越好） - ❌ 成本型（越小越好）

计算步骤

Step 1：数据标准化

\[p_{ij}= \begin{cases} \dfrac{x_{ij} - \min(x_j)}{\max(x_j) - \min(x_j)} & \text{效益型} \\ \dfrac{\max(x_j) - x_{ij}}{\max(x_j) - \min(x_j)} & \text{成本型} \end{cases} \]

Step 2：计算指标比重

\[r_{ij} = \frac{p_{ij}}{\sum_{i=1}^{m} p_{ij}} \]

Step 3：计算信息熵

\[e_j = -\frac{1}{\ln m} \sum_{i=1}^{m} r_{ij} \ln r_{ij} \]

**注：当 \(r_{ij}=0\) ，定义 $ r_{ij} \ln r_{ij} = 0 $ **

Step 4：计算差异系数

\[g_j = 1 - e_j \]

Step 5：计算权重

\[w_j = \frac{g_j}{\sum_{j=1}^{n} g_j} \]

Step 6：综合得分

\[S_i = \sum_{j=1}^{n} w_j \cdot p_{ij} \]

使用注意事项

1. 数据要求：

   样本量 \(m ≥ 5\)

   指标区分度：\(g_j \geq 0.1\) （若差异度 \(g_j=1-e_j<0.1\) 需要评估是否保留）

2. 权重修正：当结果与常识不符时，建议组合AHP（层次分析法）等主观赋权法

3. 结果验证：敏感性分析（\(±10\%\)数据扰动检验排名稳定性）

模板使用说明：

1. **填充原始数据矩阵 \(X\) 和指标类型 \(types\) **
2. 运行标准化→熵值→权重→得分的计算流程
3. 输出权重表和得分排名
4. 结合问题背景解释结果
5. 若指标数量过多可在Step2前加入主成分分析（PCA）预处理模块

Code（Python） :

def entropy_weight(X, types):
    """
    熵权法计算函数 - 模板
    X: 原始数据矩阵(m×n)
    types: 指标类型列表[1,0,...] (1=效益型, 0=成本型)
    返回: 权重、综合得分、信息熵、差异系数
    """
    m, n = X.shape
    P = np.zeros((m, n))

    # 数据标准化
    for j in range(n):
        val = X[:, j]
        mx = np.max(val)
        mn = np.min(val)
        if types[j] == 1:  # 效益型
            if mx == mn:
                P[:, j] = 0.0
            else:
                P[:, j] = (val - mn) / (mx - mn)
        else:  # 成本型
            if mx == mn:
                P[:, j] = 0.0
            else:
                P[:, j] = (mx - val) / (mx - mn)

    # 计算指标比重
    sum = np.sum(P, axis=0)
    sum[sum == 0] = 1e-6
    R = P / sum
    R[R == 0] = 1e-6 

    # 计算信息熵
    k = 1 / np.log(m) if m > 1 else 0
    e = -k * np.sum(R * np.log(R), axis=0)

    # 计算差异系数与权重
    g = 1 - e
    sum_g = np.sum(g)
    W = g / sum_g if sum_g != 0 else np.ones(n) / n

    # 计算综合得分
    res = P @ W

    return W, res, e, g

2. TOPSIS 法

原理与数学公式

TOPSIS（Technique for Order Preference by Similarity to Ideal Solution，逼近理想解排序法）是一种多属性决策方法，通过计算评价对象与 “正理想解”（最优方案）和 “负理想解”（最劣方案）的距离，确定其综合优劣程度。距离越接近正理想解、越远离负理想解的对象，综合评价越好。

步骤与公式（结合提供的代码逻辑）：

1. 数据标准化：

   对原始数据矩阵 \(X = (x_{ij})_{m \times n}\) 进行向量归一化（代码中采用的方法）：\(X_{\text{norm}}[i,j] = \frac{x_{ij}}{\sqrt{\sum_{i=1}^m x_{ij}^2}}\) （注：标准化方法不唯一，也可采用 min-max 归一化或 z-score 标准化，需根据场景选择）。

2. 确定正负理想解：

   • 正理想解 \(Z_{\text{pos}}\)：各指标的 “最优值”（效益型指标取最大值，成本型指标取最小值）：\(Z_{\text{pos}}[j] = \begin{cases} \max(X_{\text{norm}}[:,j]) & \text{若指标为效益型（types[j]=1）} \\ \min(X_{\text{norm}}[:,j]) & \text{若指标为成本型（types[j]=0）} \end{cases}\)
   • 负理想解 \(Z_{\text{neg}}\)：各指标的 “最劣值”（与正理想解相反）：\(Z_{\text{neg}}[j] = \begin{cases} \min(X_{\text{norm}}[:,j]) & \text{若指标为效益型} \\ \max(X_{\text{norm}}[:,j]) & \text{若指标为成本型} \end{cases}\)

3. 计算距离：

   用欧氏距离计算第 \(i\) 个对象与正、负理想解的距离：\(D_{\text{pos}}[i] = \sqrt{\sum_{j=1}^n (X_{\text{norm}}[i,j] - Z_{\text{pos}}[j])^2} \quad \text{（与正理想解的距离）}\) \(D_{\text{neg}}[i] = \sqrt{\sum_{j=1}^n (X_{\text{norm}}[i,j] - Z_{\text{neg}}[j])^2} \quad \text{（与负理想解的距离）}\)

4. 计算综合得分（贴近度）：

   第 \(i\) 个对象的综合得分 \(C[i]\) 为：\(C[i] = \frac{D_{\text{neg}}[i]}{D_{\text{pos}}[i] + D_{\text{neg}}[i]}\) 其中 \(C[i] \in [0,1]\)，越接近 \(1\) 说明评价对象越优。

Code（Python）：

def topsis(X, types):
    """
    TOPSIS综合评价法
    X: 原始数据矩阵(m×n)
    types: 指标类型列表[1,0,...] (1=效益型, 0=成本型)
    返回: 综合得分、正理想解、负理想解
    """
    m, n = X.shape
    # 标准化
    X_norm = X / np.sqrt(np.sum(X**2, axis=0))
    
    # 确定正负理想解
    Z_pos = np.zeros(n)
    Z_neg = np.zeros(n)
    for j in range(n):
        if types[j] == 1:  # 效益型
            Z_pos[j] = np.max(X_norm[:, j])
            Z_neg[j] = np.min(X_norm[:, j])
        else:  # 成本型
            Z_pos[j] = np.min(X_norm[:, j])
            Z_neg[j] = np.max(X_norm[:, j])
    
    # 计算距离
    D_pos = np.sqrt(np.sum((X_norm - Z_pos)**2, axis=1))
    D_neg = np.sqrt(np.sum((X_norm - Z_neg)**2, axis=1))
    
    # 计算贴近度
    C = D_neg / (D_pos + D_neg)
    return C, Z_pos, Z_neg

七、图论与网络分析

1. 最短路径算法

（1）Dijkstra 最短路径算法

原理

单源最短路径算法，用于寻找从一个起始节点到其他所有节点的最短路径，要求图中所有边的权重非负。 核心思想是：通过贪心策略，维护一个 “已确定最短路径的节点集”，每次从 “未确定节点” 中选择距离起始节点最近的节点，并用该节点更新其邻居的距离，直到所有节点都被处理。

数学公式

• 设起始节点为 \(s\)，节点 \(u\) 到 \(s\) 的当前最短距离为 \(dist[u]\)，边 \((u, v)\) 的权重为 \(w(u, v)\)。
• 松弛操作（更新距离）：当通过节点 \(u\) 到达 \(v\) 的路径更短时，更新 \(dist[v]:dist[v] = \min(dist[v],\ dist[u] + w(u, v))\)

（2）Bellman-Ford 最短路径算法

原理

单源最短路径算法，用于寻找从一个起始节点到其他所有节点的最短路径，支持负权重边，且能检测图中是否存在 “负权重环”（即总权重为负的环，此时最短路径无意义，因为可以无限绕环减小距离）。 核心思想是：通过对所有边进行 \(V-1\) 次 “松弛操作”（V 为节点数），逐步更新最短距离；最后再额外检查一次边，若仍能松弛则说明存在负权重环。

数学公式

• 设起始节点为 \(s\)，节点 \(u\) 到 \(s\) 的距离为 \(dist[u]\)，边 \((u, v)\) 的权重为 \(w(u, v)\)。
• 松弛操作（重复 \(V-1\) 次）：\(\forall (u, v) \in E,\ dist[v] = \min(dist[v],\ dist[u] + w(u, v))\)
• 负环检测（第 \(V\) 次检查）：若存在边 \((u, v)\) 满足 \(dist[v] > dist[u] + w(u, v)\)，则存在负权重环。

（3）Floyd-Warshall 算法

原理

全源最短路径算法，用于计算图中所有节点对之间的最短路径（即任意两个节点 \((i, j)\) 之间的最短距离），支持负权重边，但不支持含负权重的环（会导致路径长度无限小）。 核心思想是：基于动态规划，通过 “中间节点” 逐步优化路径。设 \(dist[k][i][j]\) 为 “允许使用前 \(k\) 个节点作为中间节点时，\(i\) 到 \(j\) 的最短距离”，则通过递推关系更新距离。

数学公式

• 初始状态：\(dist[0][i][j] = w(i, j)\)（若 \(i\) 和 \(j\) 直接相连），否则为 \(\infty\)，且 \(dist[0][i][i] = 0\)。
• 动态规划递推（\(k\) 为中间节点）：\(dist[k][i][j] = \min(dist[k-1][i][j],\ dist[k-1][i][k] + dist[k-1][k][j])\)
• 最终结果：\(dist[V][i][j]\) 即为 \(i\) 到 \(j\) 的最短距离（\(V\) 为节点总数）。

Code（Python）：

def dijkstra_shortest_path(graph, start, end, weight='weight'):
    """
    Dijkstra最短路径算法
    graph: 图对象 (networkx.Graph 或 networkx.DiGraph)
    start: 起始节点
    end: 目标节点
    weight: 边权重属性
    返回: 最短路径列表, 路径长度
    """
    try:
        path = nx.dijkstra_path(graph, start, end, weight=weight)
        length = nx.dijkstra_path_length(graph, start, end, weight=weight)
        return path, length
    except nx.NetworkXNoPath:
        return [], float('inf')
    except Exception as e:
        print(f"Dijkstra算法错误: {str(e)}")
        return [], float('inf')

def bellman_ford_shortest_path(graph, start, end, weight='weight'):
    """
    Bellman-Ford最短路径算法 (支持负权重)
    graph: 图对象
    start: 起始节点
    end: 目标节点
    weight: 边权重属性
    返回: 最短路径列表, 路径长度
    """
    try:
        path = nx.bellman_ford_path(graph, start, end, weight=weight)
        length = nx.bellman_ford_path_length(graph, start, end, weight=weight)
        return path, length
    except nx.NetworkXNoPath:
        return [], float('inf')
    except nx.NetworkXUnbounded:
        print("存在负权重环")
        return [], float('-inf')
    except Exception as e:
        print(f"Bellman-Ford算法错误: {str(e)}")
        return [], float('inf')

def floyd_warshall_all_pairs(graph, weight='weight'):
    """
    Floyd-Warshall所有节点对最短路径
    graph: 图对象
    weight: 边权重属性
    返回: 距离矩阵, 前驱矩阵
    """
    try:
        dist, path = nx.floyd_warshall_predecessor_and_distance(graph, weight=weight)
        return dist, path
    except Exception as e:
        print(f"Floyd-Warshall算法错误: {str(e)}")
        return {}, {}

2. 最小生成树

（1）Kruskal 算法（最小生成树）

核心原理

Kruskal 算法是一种基于「边排序」和「避环」策略的最小生成树（Minimum Spanning Tree, MST）算法。其核心步骤为：

1. 排序所有边：将图中所有边按权重从小到大排序；
2. 贪心选边：依次选取权重最小的边，若该边连接的两个顶点不在同一个连通分量中（即不会形成环），则将其加入 MST；
3. 终止条件：当加入的边数为「顶点数 - 1」时，生成树构建完成（包含所有顶点且无环）。

实现中需用并查集（Union-Find） 数据结构高效判断两顶点是否连通（避免环）。

数学表达

设图 \(G = (V, E)\)，其中 \(V\) 为顶点集，\(E\) 为边集，边 \(e \in E\) 的权重为 \(w(e)\)。 MST 的目标是找到子集 \(T \subseteq E\)，满足：

• \(T\) 包含所有顶点（即 \(T\) 连接 \(V\) 中所有顶点）；
• \(T\) 无环（树结构）；
• 总权重最小：\(\min \sum_{e \in T} w(e)\)。

（2）Prim 算法（最小生成树）

核心原理

Prim 算法是一种基于「顶点扩展」的 MST 算法，从一个起始顶点开始，逐步「生长」生成树：

1. 初始化：任选一个顶点作为起点，标记为「已加入 MST」，初始 MST 仅包含该顶点；
2. 贪心选边：在所有连接「MST 内顶点」和「MST 外顶点」的边中，选择权重最小的边，将对应的外顶点加入 MST；
3. 迭代扩展：重复步骤 2，直到所有顶点都被加入 MST（加入的边数为「顶点数 - 1」）。

数学表达

与 Kruskal 算法一致，目标是最小化生成树的总权重：\(\min \sum_{e \in T} w(e)\)，其中 \(T\) 为 MST 的边集，满足包含所有顶点且无环。

Code（Python）：

def kruskal_mst(graph, weight='weight'):
    """
    Kruskal最小生成树算法
    graph: 图对象
    weight: 边权重属性
    返回: 最小生成树图
    """
    try:
        mst = nx.minimum_spanning_tree(graph, weight=weight, algorithm='kruskal')
        total_weight = sum(data[weight] for u, v, data in mst.edges(data=True))
        print(f"最小生成树总权重: {total_weight:.2f}")
        return mst
    except Exception as e:
        print(f"Kruskal算法错误: {str(e)}")
        return None

def prim_mst(graph, weight='weight'):
    """
    Prim最小生成树算法
    graph: 图对象
    weight: 边权重属性
    返回: 最小生成树图
    """
    try:
        mst = nx.minimum_spanning_tree(graph, weight=weight, algorithm='prim')
        total_weight = sum(data[weight] for u, v, data in mst.edges(data=True))
        print(f"最小生成树总权重: {total_weight:.2f}")
        return mst
    except Exception as e:
        print(f"Prim算法错误: {str(e)}")
        return None

八、统计分析方法

1. 主成分回归 (PCR)

核心原理

主成分回归（Principal Component Regression, PCR）是一种将主成分分析（PCA）与线性回归相结合的降维回归方法。其核心思想是：

1. 数据降维：通过 PCA 将原始高维特征矩阵 X 转换为低维主成分 Z，主成分是原始特征的线性组合，且彼此正交（不相关）；
2. 回归建模：在降维后的主成分 Z 上构建线性回归模型，预测目标变量 y。

PCR 通过保留方差最大的主成分，在减少特征维度的同时保留数据的主要信息，从而缓解多重共线性问题，提高模型泛化能力。

数学公式

1. PCA 降维： 设原始特征矩阵 \(X \in {R}^{n \times p}\)（\(n\) 为样本数，\(p\) 为特征数），通过 PCA 转换为：\(Z = XW\) 其中 \(W \in \mathbb{R}^{p \times k}\) 是主成分载荷矩阵（\(k\) 为保留的主成分数，\(k \leq p\)），\(Z \in {R}^{n \times k}\) 是降维后的主成分矩阵。
2. 线性回归： 在主成分 Z 上建立线性模型：\(y = Z\beta + \epsilon = (XW)\beta + \epsilon\) 其中 \(\beta \in \mathbb{R}^{k}\) 是回归系数，\(\epsilon\) 是误差项。 通过最小二乘法估计 \(\beta：\hat{\beta} = (Z^TZ)^{-1}Z^Ty\) 最终原特征的回归系数为：\(\hat{\alpha} = W\hat{\beta}\)

Code（Python）：

from sklearn.decomposition import PCA
from sklearn.linear_model import LinearRegression

def pcr_regression(X, y, n_components=5):
    """
    主成分回归
    X: 自变量矩阵
    y: 因变量
    n_components: 保留主成分数
    返回: 回归系数、预测值
    """
    pca = PCA(n_components=n_components)
    X_pca = pca.fit_transform(X)
    model = LinearRegression()
    model.fit(X_pca, y)
    coef = model.coef_
    predictions = model.predict(X_pca)
    return coef, predictions

2. 因子分析 (FA)

核心原理

因子分析（Factor Analysis, FA）是一种数据降维与潜在变量提取方法，其核心假设是：观测变量 \(X\) 的变异可由少数不可观测的潜在因子 \(F\) 和特殊因子 \(\epsilon\) 解释。FA 通过寻找潜在因子，揭示数据的内在结构和变量间的关联性。

与 PCA 不同，FA 是一种生成模型，强调数据的生成机制（而非仅数据变换）。

数学公式

1. 因子模型： 对观测变量 \(X \in {R}^{p}\)，假设存在 \(k\) 个潜在因子 \(F \in {R}^{k}\)（\(k < p\)），满足：\(X = \Lambda F + \epsilon\) 其中：
   • \(\Lambda \in {R}^{p \times k}\) 是因子载荷矩阵，表示观测变量与潜在因子的关联强度；
   • \(\epsilon \in {R}^{p}\) 是特殊因子（每个变量特有，互不相关），通常假设 \(\epsilon \sim N(0, \Psi)\)，$\Psi $ 是对角协方差矩阵。
2. 参数估计： 通过最大化似然函数或最小化重构误差估计 \(\Lambda\) 和 \(\Psi\)，常用方法包括：
   • 最大似然估计（MLE）；
   • 主成分法（通过 PCA 载荷矩阵初始化）；
   • 迭代最小二乘法（如 MINRES）。
3. 因子得分： 估计潜在因子 \(F\) 的值（即样本在各因子上的得分），常用方法为：\(\hat{F} = (X - \hat{\mu}) \hat{\Lambda}(\hat{\Lambda}^T\hat{\Lambda})^{-1}\) 其中 \(\hat{\mu}\) 是观测变量的均值。

Code（Python）：

from sklearn.decomposition import FactorAnalysis

def factor_analysis(X, n_factors=3):
    """
    因子分析
    X: 原始数据矩阵
    n_factors: 提取因子数
    返回: 因子载荷矩阵、因子得分
    """
    fa = FactorAnalysis(n_components=n_factors)
    fa.fit(X)
    loadings = fa.components_.T
    scores = fa.transform(X)
    return loadings, scores

九、深度学习模型

1. LSTM 时间序列预测

核心原理

长短期记忆网络（Long Short-Term Memory, LSTM）是一种特殊的循环神经网络（RNN），专门设计用于处理序列数据中的长期依赖关系。传统 RNN 在处理长序列时易出现梯度消失或爆炸问题，而 LSTM 通过引入门控机制解决了这一问题。

LSTM 的核心组件是记忆单元和三个门控结构：

1. 输入门（Input Gate）：控制新输入信息的多少进入记忆单元；
2. 遗忘门（Forget Gate）：控制上一时刻的记忆单元状态保留多少；
3. 输出门（Output Gate）：控制当前记忆单元状态输出多少到隐藏状态。

数学公式

对于时刻 \(t\) 的输入 \(x_t\) 和上一时刻的隐藏状态 \(h_{t-1}\)，LSTM 的计算流程为：

1. 门控计算：\(\begin{aligned} i_t &= \sigma(W_i x_t + U_i h_{t-1} + b_i) \\ f_t &= \sigma(W_f x_t + U_f h_{t-1} + b_f) \\ o_t &= \sigma(W_o x_t + U_o h_{t-1} + b_o) \\ \tilde{C}_t &= \tanh(W_C x_t + U_C h_{t-1} + b_C) \end{aligned}\) 其中 \(\sigma\) 是 sigmoid 函数，将值压缩到 [0,1] 区间；\(\tanh\) 将值压缩到 [-1,1] 区间。
2. 记忆单元更新：\(C_t = f_t \odot C_{t-1} + i_t \odot \tilde{C}_t\) \(\odot\) 表示逐元素乘法，遗忘门 \(f_t\) 决定保留多少历史状态 \(C_{t-1}\)，输入门 \(i_t\) 决定添加多少新候选状态 \(\tilde{C}_t\)。
3. 隐藏状态输出：\(h_t = o_t \odot \tanh(C_t)\)

Code（Python）：

from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import LSTM, Dense

def lstm_forecast(data, look_back=10, epochs=50):
    """
    LSTM时间序列预测
    data: 时间序列数据
    look_back: 回溯步数
    epochs: 训练轮数
    返回: 预测模型、预测结果
    """
    # 数据预处理
    X, y = [], []
    for i in range(len(data)-look_back):
        X.append(data[i:i+look_back])
        y.append(data[i+look_back])
    X = np.array(X)
    y = np.array(y)
    
    # 构建LSTM模型
    model = Sequential()
    model.add(LSTM(50, activation='relu', input_shape=(look_back, 1)))
    model.add(Dense(1))
    model.compile(optimizer='adam', loss='mse')
    
    # 训练模型
    model.fit(X, y, epochs=epochs, verbose=0)
    
    # 预测
    predictions = model.predict(X)
    return model, predictions.flatten()

2. 卷积神经网络 (CNN)

核心原理

卷积神经网络（Convolutional Neural Network, CNN）是一种专门为处理网格结构数据（如图像、音频）设计的深度学习模型。其核心思想是通过卷积层自动提取数据的局部特征，减少参数数量并提高模型效率。

CNN 的主要组件包括：

1. 卷积层（Convolutional Layer）：使用多个卷积核（滤波器）滑动扫描输入，提取局部特征；
2. 池化层（Pooling Layer）：降维并保留主要特征（如最大池化、平均池化）；
3. 全连接层（Fully Connected Layer）：将提取的特征映射到分类或回归输出。

数学公式

1. 卷积操作： 对输入 \(X\) 和卷积核 \(W\)，输出特征图 \(Y\) 的计算为：\(Y[i,j] = \sum_{m,n} X[i+m,j+n] \cdot W[m,n] + b\) 其中 \(b\) 是偏置，卷积核 \(W\) 在输入上滑动并逐元素相乘求和。
2. 激活函数： 通常在卷积后应用非线性激活函数（如 ReLU）：\(Y = \max(0, X \ast W + b)\)
3. 池化操作（以最大池化为例）： 将输入划分为不重叠的区域，取每个区域的最大值作为输出：\(Y[i,j] = \max_{m,n \in \text{window}} X[i+m,j+n]\)

Code（Python）：

from tensorflow.keras.layers import Conv2D, MaxPooling2D, Flatten

def cnn_classifier(input_shape, num_classes):
    """
    CNN图像分类模板
    input_shape: 输入图像尺寸
    num_classes: 分类数量
    返回: CNN模型
    """
    model = Sequential([
        Conv2D(32, (3,3), activation='relu', input_shape=input_shape),
        MaxPooling2D((2,2)),
        Conv2D(64, (3,3), activation='relu'),
        MaxPooling2D((2,2)),
        Flatten(),
        Dense(128, activation='relu'),
        Dense(num_classes, activation='softmax')
    ])
    model.compile(optimizer='adam',
                 loss='sparse_categorical_crossentropy',
                 metrics=['accuracy'])
    return model

十、创新新模型应用——多模型结合

1. 神经网络特征工程 + 线性回归

核心原理

这个方法结合了神经网络的特征提取能力和线性回归的可解释性，通过以下步骤实现：

1. 数据标准化：使用StandardScaler将输入特征缩放为均值 0、方差 1，加速神经网络训练；
2. 神经网络特征提取：利用多层感知机（MLP）自动学习输入特征的非线性变换，生成高级抽象特征；
3. 线性回归建模：将神经网络的输出作为新特征，输入到线性回归模型中进行最终预测。

这种组合方式既利用了神经网络的强大表达能力，又保留了线性模型的可解释性和计算效率。

数学公式

1. 神经网络特征提取： 设输入特征为 \(X \in \mathbb{R}^{n \times p}\)（\(n\) 为样本数，\(p\) 为特征数），MLP 通过 L 层非线性变换生成新特征 Z：\(\begin{aligned} Z_1 &= \sigma(W_1 X + b_1) \\ Z_2 &= \sigma(W_2 Z_1 + b_2) \\ &\vdots \\ Z_L &= W_L Z_{L-1} + b_L \end{aligned}\) 其中 \(W_i\) 和 \(b_i\) 是第 \(i\) 层的权重和偏置，\(\sigma\) 是 ReLU 激活函数（\(\sigma(x) = \max(0, x)\)）。
2. 线性回归预测： 将神经网络输出 \(Z_L\) 作为新特征，通过线性回归预测目标变量 \(y\)：\(\hat{y} = Z_L \beta + \epsilon\) 其中 \(\beta\) 是回归系数，\(\epsilon\) 是误差项。

Code（Python）：

from sklearn.neural_network import MLPRegressor
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

def nn_feature_regression(X, y, hidden_layers=(50,)):
    """
    神经网络特征工程+线性回归
    1. 使用神经网络提取高级特征
    2. 在线性回归中使用这些特征
    返回: 组合模型、预测结果
    """
    # 创建组合模型
    model = make_pipeline(
        StandardScaler(),
        MLPRegressor(hidden_layer_sizes=hidden_layers, activation='relu', max_iter=1000),
        LinearRegression()
    )
    model.fit(X, y)
    predictions = model.predict(X)
    return model, predictions

2. 梯度提升回归树（GBRT）

核心原理

梯度提升回归树（Gradient Boosting Regression Tree, GBRT）是一种集成学习方法，通过迭代训练多个弱学习器（决策树）并组合它们的预测结果，形成一个强大的回归模型。其核心思想是迭代地拟合前一轮预测的残差，逐步减小误差。

GBRT 属于Boosting框架，与随机森林（Bagging）的并行训练不同，GBRT 采用串行训练，每棵树专注于修正前一棵树的错误。

数学公式

给定训练数据 \(\{(x_i, y_i)\}_{i=1}^n\)，GBRT 的目标是构建一个加法模型：\(F(x) = \sum_{m=1}^M \gamma_m h_m(x)\) 其中：

• \(h_m(x)\) 是第 \(m\) 棵决策树（弱学习器）；
• \(\gamma_m\) 是对应的权重（学习率）；
• \(M\) 是树的总数（迭代次数）。

迭代过程：

1. 初始化模型：\(F_0(x) = \text{argmin}_c \sum_{i=1}^n L(y_i, c)\) （通常取目标值的均值，\(L\) 为损失函数，如均方误差 \(MSE\)）。
2. 对 \(m = 1, 2, \dots, M\)：
   • 计算负梯度（即残差）：\(r_{im} = -\left[ \frac{\partial L(y_i, F(x_i))}{\partial F(x_i)} \right]_{F=F_{m-1}}\)
   • 拟合一棵决策树 \(h_m(x)\) 到残差 \(r_{im}\)；
   • 通过线搜索确定最优步长 \(\gamma_m\)：\(\gamma_m = \text{argmin}_\gamma \sum_{i=1}^n L(y_i, F_{m-1}(x_i) + \gamma h_m(x_i))\)
   • 更新模型：\(F_m(x) = F_{m-1}(x) + \gamma_m h_m(x)\)
3. 最终预测：\(\hat{y} = F_M(x)\)

Code（Python）：

from sklearn.ensemble import GradientBoostingRegressor

def gradient_boosting_regression(X, y, n_estimators=100):
    """
    梯度提升回归树（GBRT）
    结合多个决策树形成强模型
    返回: 模型对象、特征重要性
    """
    model = GradientBoostingRegressor(n_estimators=n_estimators)
    model.fit(X, y)
    feature_importance = model.feature_importances_
    return model, feature_importance

3. KMeans 聚类 + 线性回归

核心原理

这个方法将聚类分析与回归建模相结合，通过以下步骤实现：

1. 数据分群：使用 K-means 算法将样本分为多个簇（clusters），使得簇内样本相似度高，簇间差异大；
2. 分群建模：对每个簇单独训练一个线性回归模型，捕捉簇内数据的局部模式；
3. 预测整合：根据样本所属簇，使用对应的回归模型进行预测。

这种方法假设数据中存在异质性（即不同子群体的回归关系不同），通过分群建模可提高模型的拟合精度和解释性。

数学公式

1. K-means 聚类： 将 \(n\) 个样本 \(X = \{x_1, x_2, \dots, x_n\}\) 划分为 \(k\) 个簇 \(C = \{C_1, C_2, \dots, C_k\}\)，使得簇内平方和最小：\(\min_{C} \sum_{i=1}^k \sum_{x \in C_i} \| x - \mu_i \|^2\) 其中 \(\mu_i\) 是簇 \(C_i\) 的质心（均值向量）。
2. 分群回归： 对每个簇 \(C_i\)，训练线性回归模型：\(y = X_i \beta_i + \epsilon_i \quad (x \in C_i)\) 其中 \(\beta_i\) 是第 \(i\) 个簇的回归系数，通过最小化残差平方和估计：\(\hat{\beta}_i = (X_i^T X_i)^{-1} X_i^T y_i\)
3. 预测： 对新样本 \(x\)，先确定其所属簇 \(C_j\)，再用对应模型预测：\(\hat{y} = x^T \hat{\beta}_j\)

Code（Python）:

from sklearn.cluster import KMeans
from sklearn.linear_model import LinearRegression

def cluster_based_regression(X, y, n_clusters=3):
    """
    聚类+回归分群建模
    1. 对样本聚类
    2. 每个聚类单独建立回归模型
    返回: 聚类模型字典、预测结果
    """
    kmeans = KMeans(n_clusters=n_clusters)
    clusters = kmeans.fit_predict(X)
    
    models = {}
    predictions = np.zeros_like(y)
    
    for i in range(n_clusters):
        mask = (clusters == i)
        if np.sum(mask) > 1:  # 确保有足够样本
            cluster_model = LinearRegression()
            cluster_model.fit(X[mask], y[mask])
            models[i] = cluster_model
            predictions[mask] = cluster_model.predict(X[mask])
    
    return models, predictions

4. CNN-SVM 组合模型

核心原理

这个方法将卷积神经网络（CNN）的特征提取能力与支持向量机（SVM）的分类能力相结合，通过以下步骤实现：

1. CNN 特征提取：使用预训练或自定义的 CNN 架构（如 LeNet）从图像中提取高级特征；
2. 特征向量化：通过 Flatten 层将 CNN 输出的特征图转换为一维向量；
3. SVM 分类：将提取的特征向量输入 SVM，训练分类器。

这种组合利用了 CNN 在图像特征提取上的优势（自动学习层次化特征）和 SVM 在小样本分类上的优势（寻找最大间隔超平面）。

数学公式

1. CNN 特征提取： 输入图像 \(X \in {R}^{h \times w \times c}\)，通过卷积、池化和激活函数生成特征向量 z：\(z = \text{Flatten}(\text{MaxPool}(\text{ReLU}(W_2 \ast \text{MaxPool}(\text{ReLU}(W_1 \ast X + b_1)) + b_2)))\) 其中 \(W_1, W_2\) 是卷积核权重，\(b_1, b_2\) 是偏置，\(\ast\) 表示卷积操作。
2. SVM 分类： 在特征空间中寻找最优超平面 \((w, b)\)，使得不同类别的样本间隔最大：\(\min_{w, b} \frac{1}{2} \|w\|^2 + C \sum_{i=1}^n \xi_i\) 约束条件：\(y_i(w^T \phi(z_i) + b) \geq 1 - \xi_i, \quad \xi_i \geq 0\) 其中 \(\phi(z_i)\) 是特征映射（在 RBF 核中隐式定义），\(\xi_i\) 是松弛变量，\(C\) 是正则化参数。
3. 核函数（RBF）：\(K(z_i, z_j) = \exp\left(-\gamma \|z_i - z_j\|^2\right)\) 其中 \(\gamma\) 是核系数。

Code（Python）：

from tensorflow.keras.models import Model
from tensorflow.keras.layers import Input, Conv2D, MaxPooling2D, Flatten
from sklearn.svm import SVC

def cnn_svm_classifier(X, y, input_shape):
    """
    CNN特征提取 + SVM分类
    1. 使用CNN提取图像特征
    2. 将特征输入SVM进行分类
    返回: 组合模型、准确率
    """
    # 构建CNN特征提取器
    inputs = Input(shape=input_shape)
    x = Conv2D(32, (3,3), activation='relu')(inputs)
    x = MaxPooling2D((2,2))(x)
    x = Conv2D(64, (3,3), activation='relu')(x)
    x = MaxPooling2D((2,2))(x)
    x = Flatten()(x)
    feature_extractor = Model(inputs, x)
    
    # 提取特征
    features = feature_extractor.predict(X)
    
    # SVM分类
    svm = SVC(kernel='rbf', C=1.0)
    svm.fit(features, y)
    accuracy = svm.score(features, y)
    
    return (feature_extractor, svm), accuracy