【D题解题思路】2025深圳杯&东三省D题解题思路+可运行代码参考（无偿分享）

最新推荐文章于 2025-05-01 16:43:55 发布

睿森竞赛

最新推荐文章于 2025-05-01 16:43:55 发布

阅读量637

点赞数 14

分类专栏：深圳杯数学建模东三省数学建模文章标签：数学建模深圳杯数学建模东三省数学建模

本文链接：https://blog.csdn.net/RS_handy/article/details/147644196

版权

深圳杯数学建模同时被 2 个专栏收录

3 篇文章

订阅专栏

东三省数学建模

3 篇文章

订阅专栏

注：该内容由“数模加油站”原创，无偿分享，可以领取参考但不要利用该内容倒卖，谢谢！

D 题法医物证多人身份鉴定问题

问题1 贡献者人数识别建模

问题 1 分析

在多人混合的STR图谱分析中，判断贡献者人数是最基础也是最关键的步骤之一。由于每位个体在每个位点上最多拥有两个等位基因，因此当STR图谱中某个位点所观测到的等位基因峰数超过两个时，便可能暗示该样本为多人混合体。STR图谱的每个主峰由其片段长度（size）和峰高（height）组成，其中峰数的多少与混合人数之间存在潜在的正相关关系。因此，本问题可转化为对多个基因座上峰值数量的统计建模问题。具体而言，我们可通过构造统计特征（如每个位点的等位基因数目、所有位点的最大/平均峰数等）来进行人数估计。此外，为了克服数据中的波动性和测序误差，可以引入概率模型或基于聚类的方法，如使用高斯混合模型（GMM）在height维度上分析不同峰值类别的聚集情况，进一步通过信息准则（如AIC或BIC）选择最优的人数估计。在此基础上，可结合部分已知标签数据构建监督学习模型进行验证与校准，最终实现对未知混合图谱的贡献者人数的准确识别。

解题思路：

1.1 问题背景与建模目标

在法医物证DNA分析中，STR（Short Tandem Repeat）图谱的解读是关键环节。当存在多人混合样本时，首要任务即是判断该混合样本中究竟由几位个体贡献，这一判断直接关系到后续的比例估计与基因型溯源准确性。给定一个混合STR图谱样本，其每个位点上会出现多个等位基因峰，每个峰由其片段长度（size）和峰高（height）共同定义。由于每个正常个体在一个基因座上最多只表现两个等位基因，因此若观察到某个位点的峰数超过两个，通常表明该样本为混合样本。基于此，我们可以统计多个基因座上的峰数特征，构建数学模型以估计样本中的贡献者人数。

1.2 数据结构与预处理

附件1提供了多个混合STR图谱样本文件，其编码中已嵌入真实贡献者人数，例如 "40_41-1;4" 表示该样本由编号为40和41的两人混合构成。在数据读取后，我们需对每一个样本，提取其在各个marker（基因座）下的主峰数（即等位基因数）与峰高（height）信息。设定样本 S 包含 M 个marker，记第 i 个marker上的峰数为 $k_i$ ，其峰高向量为 $\mathbf{h}i = [h{i1}, h_{i2}, \dots, h_{ik_i}]$ 。

1.3 特征构建与统计建模

我们以每个样本 S 的marker峰数分布为基础，构建特征向量 $\mathbf{x} \in \mathbb{R}^d$ 表征该样本的混合特征。常见的构造方式如下：

平均峰数： $\mu_k = \frac{1}{M} \sum_{i=1}^M k_i$

最大峰数： $k_{\text{max}} = \max { k_1, k_2, \dots, k_M }$

峰数方差： $\sigma_k^2 = \frac{1}{M} \sum_{i=1}^M (k_i - \mu_k)^2$

在观察中，随着人数 n 的增加，每个位点的峰数趋于增加。因此我们初步假设混合人数 n 与平均峰数 $\mu_k$ 呈正相关关系，可建立如下线性估计模型：

$\hat{n} = a \cdot \mu_k + b$

其中 a 与 b 为通过训练样本拟合得到的参数。为提升鲁棒性，我们也可引入非线性建模方式，如多项式回归或决策树模型。

1.4 峰高分布的辅助建模

除了峰数信息外，峰高的分布亦可提供有效判别信息。在多人混合中，因个体贡献不等，部分峰高显著不同，故可通过峰高间方差判断是否存在多源干扰。定义每个位点峰高的变异系数（Coefficient of Variation）为：

$CV_i = \frac{\text{std}(\mathbf{h}_i)}{\text{mean}(\mathbf{h}_i)}$

然后求出整个样本的平均变异系数：

$\overline{CV} = \frac{1}{M} \sum_{i=1}^{M} CV_i$

我们可以构建一个二元特征向量 $\mathbf{x} = [\mu_k, \overline{CV}]$ 输入分类模型，如支持向量机（SVM）或随机森林（RF），训练以预测贡献者人数 n。目标函数形式为：

$\hat{n} = f(\mu_k, \overline{CV}; \theta)$

其中 \theta 为模型参数，f 为监督学习模型结构。

1.5 无监督学习方法与信息准则选择

考虑到部分样本缺乏真实标签或存在不确定性，我们还可引入无监督聚类建模。对每个marker的峰高向量 $\mathbf{h}_i$ ，使用高斯混合模型（GMM）进行聚类，估计峰群中隐含的个体成分数。设GMM的混合分布为：

$p(h) = \sum_{j=1}^{K} \pi_j \cdot \mathcal{N}(h|\mu_j, \sigma_j^2)$

其中 K 即为拟合人数的估计值， $\pi_j$ 为权重参数。最终选择最优 K 可基于BIC准则进行：

$\text{BIC}(K) = -2 \cdot \log \mathcal{L}_K + K \cdot \log(N)$

其中 $\mathcal{L}_K$ 是在 K 个分布下的最大似然值，N 为观测数据数。

1.6 评估与精度指标

为了验证模型的准确性，我们将在已知人数的样本集上评估预测准确率、精度、召回率和F1-score等标准分类指标。若实际人数为 n，预测结果为 \hat{n}，则定义准确率为：

$\text{Accuracy} = \frac{1}{T} \sum_{i=1}^{T} \mathbb{I}(\hat{n}_i = n_i)$

其中 $\mathbb{I}(\cdot)$ 为指示函数，T 为测试样本数。

Python代码：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import re
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import classification_report, confusion_matrix, accuracy_score

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# Excel 路径：请根据本地路径自行修改
file_path = '附件1：不同人数的STR图谱数据.xlsx'
xls = pd.ExcelFile(file_path)
sheet_names = xls.sheet_names

# 定义等位基因和峰高字段
allele_cols = [f'Allele {i}' for i in range(1, 7)]
height_cols = [f'Height {i}' for i in range(1, 7)]

# 存储所有样本特征
data_records = []

for name in sheet_names:
    try:
        df = xls.parse(name)
        if not set(allele_cols).issubset(df.columns):
            continue  # 忽略字段不完整的sheet

        # 每行统计该marker上非空等位基因数作为峰数
        df['PeakCount'] = df[allele_cols].notna().sum(axis=1)

        # 样本统计特征
        avg_peak = df['PeakCount'].mean()
        max_peak = df['PeakCount'].max()
        std_peak = df['PeakCount'].std()

        # 从sheet名中提取贡献者人数
        match = re.search(r'-([\d_]+)', name)
        if match:
            ids = match.group(1).replace('_', ';').split(';')
            num_contributors = len([i for i in ids if i.isdigit()])
        else:
            num_contributors = 1

        data_records.append({
            'Sample': name,
            'AvgPeak': avg_peak,
            'MaxPeak': max_peak,
            'StdPeak': std_peak,
            'NumContributors': num_contributors
        })

    except Exception as e:
        print(f"跳过 {name}，错误：{e}")

# 构建DataFrame
df_features = pd.DataFrame(data_records)
print(df_features.head())

# 特征 + 标签划分
X = df_features[['AvgPeak', 'MaxPeak', 'StdPeak']]
y = df_features['NumContributors']

# 若样本数量允许，划分训练集和测试集
if len(df_features) > 5:
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

    # 建模训练
    clf = RandomForestClassifier(n_estimators=100, random_state=42)
    clf.fit(X_train, y_train)
    y_pred = clf.predict(X_test)

    # 输出结果
    print(f"模型准确率：{accuracy_score(y_test, y_pred):.4f}")
    print("分类报告：\n", classification_report(y_test, y_pred))

    # 混淆矩阵可视化
    plt.figure(figsize=(8, 6))
    sns.heatmap(confusion_matrix(y_test, y_pred), annot=True, fmt='d', cmap='Blues',
                xticklabels=sorted(y.unique()), yticklabels=sorted(y.unique()))
    plt.title('混淆矩阵：预测STR混合人数')
    plt.xlabel('预测值')
    plt.ylabel('真实值')
    plt.tight_layout()
    plt.savefig('混淆矩阵-预测STR混合人数.png')

# 可视化特征分布图
plt.figure(figsize=(10, 6))
sns.scatterplot(data=df_features, x='AvgPeak', y='StdPeak', hue='NumContributors', palette='viridis', s=100)
plt.title('STR样本峰数统计特征分布')
plt.xlabel('平均峰数')
plt.ylabel('峰数标准差')
plt.tight_layout()
plt.savefig('STR样本峰数统计特征分布.png')

问题2贡献者比例识别建模

问题 2 分析

在识别出混合样本中参与者人数后，进一步精确判断各个贡献者的混合比例对于后续基因型的还原和个体溯源具有重要意义。STR图谱中每个峰的height反映了该等位基因对应DNA片段的相对丰度，理论上，其值应与该等位基因所属个体的贡献比例成正比。因此，可以从各个位点的峰高出发，构建混合峰与单个贡献者DNA量之间的函数关系，并以此建立目标优化模型进行比例反演。在建模过程中，需综合考虑多个基因座的共性，并对重叠等位基因峰进行去卷积。优化目标可以设定为混合图谱的预测峰高与观测峰高之间的最小二乘误差，在约束条件下求解贡献者比例向量。由于比例求解为一个典型的非负且和为1的约束优化问题，因此可引入如粒子群优化（PSO）、差分进化（DE）等全局优化算法，提高求解精度和稳定性。该模型的输出比例可通过与已知比例样本的比较进行误差评估，如计算均方根误差（RMSE）或平均绝对误差（MAE）来量化识别性能。

解题思路：

2.1 问题背景与建模目标

在法医DNA分析中，STR（Short Tandem Repeat）图谱不仅体现了DNA片段在各个基因位点的等位基因分布结构，还通过峰高（height）间接反映了各贡献者在混合样本中的相对比例。每个个体在一个基因座上最多携带两个等位基因，因此若知道混合样本中贡献者人数 n，则可以推断该STR图谱在各个位点上的峰高是这 n 位个体在相应基因型基础上“贡献比例加权”后的结果。

本题的目标是：在已知贡献者人数的前提下，依据附件2中多个不同混合比例样本的STR峰高数据，建立数学模型以识别每位贡献者在样本中所占的混合比例，进而提升后续基因型推理的准确性。

2.2 数据结构理解与变量定义

根据附件2的格式，每个样本包含多个marker（基因座），每个marker下存在若干个峰值（等位基因），每个峰值具有 size（片段长度）和 height（峰高）属性。假设某一混合样本包含 n 位个体，记其混合比例为向量：

$\mathbf{p} = [p_1, p_2, \dots, p_n]^\top, \quad$ 其中 $\quad p_i \geq 0, \quad \sum_{i=1}^{n} p_i = 1$

设该样本在第 j 个marker处观测到 m_j 个等位基因峰值，其峰高为：

$\mathbf{H}j^{\text{obs}} = [h{j1}, h_{j2}, \dots, h_{jm_j}]^\top$

设参与混合的 n 位个体在该marker上的理论基因型为 ${\mathbf{G}{j}^{(1)}, \dots, \mathbf{G}{j}^{(n)}}$ ，其中 $\mathbf{G}{j}^{(i)}$ 表示第 i 位个体在第 j 个marker的等位基因组合，例如 $\mathbf{G}{j}^{(i)} = [a_{j1}^{(i)}, a_{j2}^{(i)}]$ 表示其拥有两个等位基因。则混合后在marker j 上的理论峰高模型为：

$\mathbf{H}j^{\text{pred}} = \sum{i=1}^{n} p_i \cdot \delta(\mathbf{G}_{j}^{(i)})$

其中 \delta(\mathbf{G}_{j}^{(i)}) 表示将该个体的等位基因转换为峰值向量（按等位基因size对齐，非其等位基因位置记为0），例如对于长度为 m_j 的全局等位基因列表，其贡献为 0-1 向量。

2.3 模型构建与目标函数设计

上述过程可形式化为一个约束最小化问题：最小化预测峰高与真实观测峰高之间的误差，记为：

$\min_{\mathbf{p}} \quad L(\mathbf{p}) = \sum_{j=1}^{M} \left| \mathbf{H}j^{\text{obs}} - \sum{i=1}^{n} p_i \cdot \delta(\mathbf{G}_{j}^{(i)}) \right|_2^2$

约束条件为：

$\sum_{i=1}^{n} p_i = 1, \quad p_i \geq 0 \quad \forall i$

其中 M 为总共参与计算的marker数目。该问题本质为带约束的非负最小二乘拟合问题。由于真实STR峰高受到扩增偏差、背景噪声等因素影响，为增强鲁棒性，可以加入正则项，例如：

$\min_{\mathbf{p}} \quad L(\mathbf{p}) + \lambda \cdot |\mathbf{p}|_2^2$

其中 $\lambda$ 为正则化系数，可用交叉验证选择。此正则项用于防止比例向量的过拟合震荡，增强泛化能力。

2.4 求解策略与优化算法

该优化问题具备连续性与凸性，常见的求解方法包括：

带约束的最小二乘法（如scipy.optimize中的 $lsq_linear$ ）；

投影梯度下降法，每轮梯度更新后投影至概率单纯形上；

智能优化算法，如粒子群优化（PSO）适用于非凸扰动问题建模场景。

若采用粒子群优化，其核心为构造适应度函数：

$F(\mathbf{p}) = -L(\mathbf{p})$

粒子更新过程包括：

$p_i^{(t+1)} = p_i^{(t)} + w \cdot v_i^{(t)} + c_1 \cdot r_1 \cdot (p_{i}^{\text{best}} - p_i^{(t)}) + c_2 \cdot r_2 \cdot (p^{\text{global}} - p_i^{(t)})$

其中 w 为惯性权重， $c_1, c_2$ 为加速系数， $r_1, r_2$ 为[0,1]间随机数，通过适应度函数最大化实现比例向量的收敛搜索。

2.5 模型评估指标

在比例预测完成后，需与真实混合比例进行比对，设真实比例为 $\mathbf{p}^{\text{true}}$ ，预测比例为 $\hat{\mathbf{p}}$ ，则可用以下指标评估模型效果：

均方根误差（RMSE）：

$\text{RMSE} = \sqrt{ \frac{1}{n} \sum_{i=1}^n \left( p_i^{\text{true}} - \hat{p}_i \right)^2 }$

平均绝对误差（MAE）：

$\text{MAE} = \frac{1}{n} \sum_{i=1}^n \left| p_i^{\text{true}} - \hat{p}_i \right|$

同时也可观察混合比例预测对后续基因型匹配任务的提升程度，作为间接验证标准。

2.6 方法优势与扩展

该模型利用STR图谱中峰高的定量特征，结合线性权重叠加思想，能有效将观察峰值反推至各贡献者比例，并具备良好的可解释性与扩展性。在后续问题中（如基因型识别），该比例估计模型可直接作为已知量代入，有助于提高推理准确率。

Python代码：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import re
from scipy.optimize import minimize

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# ---------- Step 1: 加载 Excel 文件 ----------
file_path = '附件2：不同混合比例的STR图谱数据.xlsx'
xls = pd.ExcelFile(file_path)
sheet_names = xls.sheet_names

# ---------- Step 2: 提取比例向量 ----------
def extract_proportion(sheet_name):
    try:
        part = sheet_name.split('-')[1].split(';')
        nums = [float(n) for n in part]
        total = sum(nums)
        return [x / total for x in nums]
    except:
        return []

# ---------- Step 3: 查找含2人混合样本 ----------
target_sheet = None
for name in sheet_names:
    if '-' in name and ';' in name:
        p = extract_proportion(name)
        if len(p) == 2:
            target_sheet = name
            true_p = p
            break

if target_sheet is None:
    print("未找到两人混合样本，请检查附件2命名格式")
    exit()

# ---------- Step 4: 加载子表 ----------
df_dict = pd.read_excel(file_path, sheet_name=target_sheet)
if '不同比例的数据集' not in df_dict:
    print("sheet中未找到 '不同比例的数据集' 子表")
    exit()

df = df_dict['不同比例的数据集']

# ---------- Step 5: 提取有效峰高 ----------
height_cols = [col for col in df.columns if isinstance(col, str) and col.startswith("Height")]
observed_heights = []

for _, row in df.iterrows():
    heights = [row[col] for col in height_cols if pd.notna(row[col])]
    if len(heights) >= 2:
        observed_heights.append(sorted(heights[:2]))

# ---------- Step 6: 定义优化目标 ----------
def proportion_loss(p, obs_heights):
    if len(p) != 2: return 1e6
    p = np.clip(p, 0.0001, 0.9999)
    loss = 0
    for obs in obs_heights:
        obs = np.array(obs) / np.sum(obs)
        pred = np.array([p[0], p[1]])
        loss += np.sum((obs - pred) ** 2)
    return loss

# ---------- Step 7: 优化求解 ----------
x0 = [0.5, 0.5]
bounds = [(0, 1), (0, 1)]
constraints = {'type': 'eq', 'fun': lambda x: np.sum(x) - 1}

res = minimize(proportion_loss, x0, args=(observed_heights,), method='SLSQP', bounds=bounds, constraints=constraints)
predicted_p = res.x

# ---------- Step 8: 可视化结果 ----------
if len(true_p) == 2 and len(predicted_p) == 2:
    df_plot = pd.DataFrame({
        '贡献者': ['贡献者1', '贡献者2'],
        '真实比例': true_p,
        '预测比例': predicted_p
    })

    bar_width = 0.35
    index = np.arange(2)

    plt.figure(figsize=(8, 6))
    plt.bar(index, df_plot['真实比例'], bar_width, label='真实比例')
    plt.bar(index + bar_width, df_plot['预测比例'], bar_width, label='预测比例')

    plt.xlabel('贡献者')
    plt.ylabel('比例值')
    plt.title(f'样本 {target_sheet} 的贡献者比例识别结果')
    plt.xticks(index + bar_width / 2, df_plot['贡献者'])
    plt.legend()
    plt.tight_layout()
    plt.savefig('比例.png')

    print("真实比例：", true_p)
    print("预测比例：", predicted_p)
else:
    print("比例结果长度不匹配，无法绘图。")