简介：2021年华为杯数学建模竞赛D题的完整建模解答资源，面向参赛学生与药物优化建模学习者，聚焦抗乳腺癌候选药物的特征筛选、回归预测与二分类评估。压缩包共83个文件，包括32个Python计算脚本、22个CSV结果数据、14张PNG可视化图表、5个Excel分析表格及详细README与代码说明，整体31.82MB，目录沿pycode结构展开，便于对照赛题流程逐模块快速检索。目前已有3478人学习下载。内容覆盖卡方检验、RFE特征选择、XGBoost回归、SVM分类、交叉验证与超参数调优等关键环节，并给出缺失值处理、异常值检测、标准化归一化等预处理代码。辅助脚本配有散点图、箱线图与ROC曲线绘制逻辑，可帮助读者复现赛题完整方案，理解从数据清洗、模型训练到结果解读的闭环，并迁移到生物信息与药物研发场景。
1. 抗乳腺癌药物优化建模在解什么题：小样本高维数据的组合拳
2021年华为杯D题“抗乳腺癌候选药物的优化建模”，题面材料并不复杂：一组ERα靶点的生物活性值pIC50，外加729个分子描述符，样本量却只有197个化合物；另一个附件里还要在几十个ADMET属性上做二分类预测。题面短，真正难的是数据结构——特征数量远超样本数量，任何不做筛选的模型都会把噪声当信号。这不是传统意义上的大数据竞赛，反而是一次非常典型的“小样本高维”数学建模训练。它适合两类读者：一类是准备华为杯或研究生数模竞赛的学生，需要一套可复用的解题框架；另一类是日常做QSAR或药物筛选的从业者，想把特征选择、回归、分类串成一条完整流水线。这篇笔记就是我按这个流程把题目拆开、把每个阶段的坑补齐后的解答思路，直接照着复现即可。
2. 特征选择：从729个分子描述符到二十个稳健特征的三层漏斗
2.1 第一层：缺失率、常量列与极端分布清洗
拿到729列特征矩阵，第一件事不是跑模型，而是做结构性体检。那个矩阵看起来“干净”，实际上暗藏很多无用列：某些描述符在绝大多数化合物上取值完全相同，有些整列都是同一个数，有些列的缺失率超过一半。这些列留在特征集里，对树模型来说是纯噪声，会让模型反复在无意义的边界上分裂，拉低泛化能力；对线性模型来说，它们的系数会被噪声主导，解释性立刻崩掉。
清洗我按三步做，顺序不能乱。第一步统计缺失率，阈值设在30%；超过的直接删除，剩余缺失用列中位数填充。这里我不太推荐kNN插补，因为样本只有197个，kNN插补会把少量高杠杆点的偏离扩散到邻近样本，后面模型反而更飘。第二步计算每列标准差，小于0.01的视为常量列，整列剔除。第三步看偏度，对偏度绝对值大于3的列做对数变换，压制重尾对后续相关系数和LASSO的干扰。三步做完，特征数一般会从729降到450到550之间，但有效信息几乎没丢。
import numpy as np import pandas as pd df = pd.read_csv("molecular_descriptors.csv", index_col=0) y = df["pIC50"] X_raw = df.drop(columns=["pIC50"]) # 第一步：缺失率超过30%整列删除 miss_ratio = X_raw.isnull().mean() X_raw = X_raw.drop(columns=miss_ratio[miss_ratio > 0.30].index) # 第二步：标准差小于0.01视为常量列 std_vals = X_raw.std() X_raw = X_raw.drop(columns=std_vals[std_vals < 0.01].index) # 第三步：中位数填补缺失，再处理重尾列 X_raw = X_raw.fillna(X_raw.median()) skew_vals = X_raw.skew() for col in skew_vals[skew_vals.abs() > 3].index[:60]: X_raw[col] = np.log1p(X_raw[col] - X_raw[col].min() + 1)三个阈值——缺失率30%、标准差0.01、偏度3——都是经验值，可以直接套用。唯独对数变换这步要小心：log1p只对正值友好，如果一列最小值是负数，直接平移再取对数，变换后分布依然难看。我的处理是先对整列做缩尾（上下1%分位截断），再平移取log，效果比硬变换好得多，这个坑我在第一版解答里踩过。
第一层结束后，你会面对一个信号更稠密、但维度仍然很高的特征矩阵。接下来进入第二层，开始砍相关性。
2.2 第二层：相关系数排序与VIF去共线性
729列分子描述符不是互相独立的，它们是从同一个分子结构体系里导出的不同视角。比如分子体积和分子表面积天然正相关，某些拓扑指数之间也存在明显的线性依赖。把这种高度相关的特征全部保留，树模型会在它们之间随机分配分裂点，重要性排序完全失真；线性回归则陷入多重共线性，系数估计变得极不稳定，换一个训练集结果就面目全非。
第二层的做法分两步。第一步做过滤式筛选：计算每个特征与pIC50的Pearson相关系数和Spearman秩相关系数，各取绝对值，按两者中的较大值排序，保留前150到200个。为什么要同时算两个？活性与描述符的关系不总是线性的。有些描述符只在秩相关上有信号，线性相关系数算出来只有0.1左右，只看Pearson会把它们一票否决，Spearman能接住这类弱单调关系。
第二步去共线性。对保留特征做相关矩阵，每次找到相关性最高的一对特征，剔除其中与pIC50相关系数较低者，循环迭代，直到所有特征两两相关系数绝对值小于0.8。之后再看VIF，把方差膨胀因子大于10的特征逐个剔除。这一步有点笨重，但对最终模型的稳定性贡献极大。
corr_with_y = X_clean.corrwith(y).abs() drop_set = set() while True: corr_mat = X_clean.drop(columns=drop_set).corr().abs() upper = corr_mat.where(np.triu(np.ones(corr_mat.shape), k=1).astype(bool)) max_pair = np.unravel_index(np.nanargmax(upper.values), upper.shape) if upper.values[max_pair] < 0.80: break col1, col2 = upper.columns[max_pair[0]], upper.columns[max_pair[1]] if corr_with_y[col1] >= corr_with_y[col2]: drop_set.add(col2) else: drop_set.add(col1) X_second = X_clean.drop(columns=drop_set)这里有个容易被忽略的问题：当两个特征相关系数极高时，剔除哪一个是会影响模型形态的。我的原则是保留与目标相关性更高的那列，而不是保留方差更大的那列。相关性高意味着它对活性的解释更直接，方差大只说明数值分布宽，不代表信息量。
相关系数阈值0.8对分子描述符来说已经偏紧，会让剩余特征偏少，但换来的是后续LASSO和随机森林能各自稳定收敛。有些队伍用0.9甚至0.95，特征数更多、训练集拟合更好，但对噪声也更敏感。如果你目标是拿一个稳定提交结果，宁可用0.8偏紧一点，模型坏在复杂上的概率小很多。
第二层跑完，特征数一般降到60到100个。接下来第三步，用模型内嵌方法把特征压到20至40个。
2.3 第三层：LASSO与随机森林重要性交叉取交集
第二层之后的特征数量，对线性模型尚可接受，对树模型仍然偏多。样本只有197个，50个特征意味着每棵树的分裂点选择会非常稀疏，过拟合风险直线上升。所以第三层要做嵌入式筛选，让LASSO和随机森林各投一次票，再取交集。
from sklearn.linear_model import LassoCV from sklearn.ensemble import RandomForestRegressor X_mid = X_second # LASSO用内部5折交叉验证选惩罚强度，系数非零者入选 lasso = LassoCV(cv=5, random_state=42, max_iter=20000) lasso.fit(X_mid, y) coef_series = pd.Series(lasso.coef_, index=X_mid.columns) lasso_feats = coef_series[np.abs(coef_series) > 1e-6].index.tolist() # 随机森林500棵树，max_features用sqrt，降低无关特征干扰 rf = RandomForestRegressor( n_estimators=500, max_features="sqrt", random_state=42, n_jobs=-1 ) rf.fit(X_mid, y) imp = pd.Series(rf.feature_importances_, index=X_mid.columns) rf_feats = imp.sort_values(ascending=False).head(40).index.tolist() # 交集过小就回退，否则用交集+LASSO独有特征补足 inter = list(set(lasso_feats) & set(rf_feats)) extra = [f for f in lasso_feats if f not in rf_feats][:10] final_feats = list(dict.fromkeys(inter + extra)) print(f"LASSO {len(lasso_feats)} 个，RF {len(rf_feats)} 个，" f"交集 {len(inter)} 个，最终 {len(final_feats)} 个")有人会问为什么不用PCA之类的降维方法。PCA的问题在于映射后的主成分是数学上的抽象组合，丢失了“哪个描述符真正重要”的可解释性。华为杯这类数模竞赛的评分非常看重论文里“为什么选这些特征”的说理过程，PCA给不出这个答案，LASSO加随机森林则能在保留原始特征ID的同时完成筛选。
另一个实操细节：如果交集太小，比如只有三四个，我不硬凑。交集小说明两种模型对数据结构的判断差异很大，通常是因为特征集合里存在多套互相竞争的信号。此时正确的做法是回到第二层，把相关系数阈值从0.8放宽到0.9，多放一些弱相关特征进来，再跑第三层。这个回环是正常流程，不是返工。最终特征数控制在20到40个之间是一个平衡点：少于20个容易欠拟合，多于40个解释成本和过拟合风险同时上来。
3. 做厚pIC50回归模型：基线、早停与五折交叉验证的完整管线
3.1 划分、标准化与基线回归：先测水位再谈优化
特征定稿后，我先立一个不可动摇的骨架：固定随机种子的训练测试划分，外加五折交叉验证。为什么固定随机种子？小样本场景里一次划分的运气能让R²上下波动0.1，随机种子天天变，你连“哪个特征集更好”都判断不了。我习惯把random_state固定在42，20%的样本留作验证集，剩下80%做训练和交叉验证。
然后跑基线回归。基线的任务是测水位——知道最简单的模型能到什么水平，再去堆复杂度。一般我用Ridge和ElasticNet各跑一次。R²落在0.5到0.7之间都算正常水位，说明特征筛选保留下来的信息是靠谱的；如果基线R²连0.4都不到，问题基本出在特征选择阶段，先回去找特征的问题，别急着换模型。
from sklearn.model_selection import train_test_split, KFold, cross_val_score from sklearn.preprocessing import StandardScaler from sklearn.linear_model import Ridge from sklearn.metrics import r2_score X_feat = X_second[final_feats] X_train, X_val, y_train, y_val = train_test_split( X_feat, y, test_size=0.2, random_state=42 ) # 标准化只fit训练集，验证集只transform scaler = StandardScaler().fit(X_train) X_tr_s = scaler.transform(X_train) X_va_s = scaler.transform(X_val) ridge = Ridge(alpha=1.0) ridge.fit(X_tr_s, y_train) pred = ridge.predict(X_va_s) print(f"Ridge 验证集 R2={r2_score(y_val, pred):.4f}") kfold = KFold(5, shuffle=True, random_state=42) cv_ridge = cross_val_score(Ridge(alpha=1.0), X_tr_s, y_train, cv=kfold, scoring="r2") print(f"Ridge 5折平均 R2={cv_ridge.mean():.4f} (+/- {cv_ridge.std():.4f})")这段代码里最关键的一行是StandardScaler().fit(X_train)。标准化只能从训练集学均值和方差，验证集和测试集只允许transform，不允许fit。有队伍在特征筛选阶段把全量数据扔进去，标准化的统计量里混进了验证集信息，验证集R²虚高，真实提交时直接翻车。这个顺序是铁律，任何时候都不要破。
基线跑完后，我还会顺手看一个东西：Ridge系数随alpha变化的稳定性。把alpha从0.001扫到100，观察每个特征的系数是否出现符号翻转。符号翻转说明该特征在不同正则强度下作用方向不确定，这种特征即使重要性分数高，也不稳定，我一般会从最终特征集里划掉。这步用np.logspace网格扫一遍很快，价值却很大，相当于给特征集做了一次压力测试。
3.2 树模型的关键超参与早停策略：别让模型赢在训练集
基线确认信息结构没坏，接下来上树模型。XGBoost和LightGBM在这个数据量级下都很靠谱，我优先用XGBoost，因为它的内置正则和早停接口更透明，调参行为好解释。树模型在小样本上最大的敌人是过拟合，控制手段就三个：树深度、学习率、迭代轮数。
我习惯把n_estimators设大，靠早停决定停在哪里，而不是拍脑袋设个固定的500或1000。早停的验证集必须固定，不能拿交叉验证的每一折做早停——如果早停用的数据和评估用的数据是同一批，停止点就是过拟合点，结果天然乐观。
import xgboost as xgb xgb_model = xgb.XGBRegressor( n_estimators=3000, max_depth=3, learning_rate=0.02, subsample=0.8, colsample_bytree=0.8, reg_lambda=1.0, early_stopping_rounds=100, random_state=42, ) # 在固定划分的验证集上早停，记录最佳迭代轮数 xgb_model.fit( X_train, y_train, eval_set=[(X_val, y_val)], verbose=False ) print(f"最佳迭代轮数: {xgb_model.best_iteration}")这里几个参数值得展开。max_depth=3对应197个样本基本是上限，再深就过拟合。learning_rate=0.02配合3000棵候选树，学习慢但稳，不会因为单棵树权重太大而跳变。subsample=0.8和colsample_bytree=0.8是双重抽样，每棵树只看80%的样本和80%的特征，显著降低树之间的相关性，是控制方差的经典配置。reg_lambda是L2正则项，对树模型的实际约束有限，但保留它无损。
我在第一版代码里踩过一个坑：把eval_set直接传进cross_val_score的fit_params里，结果早停信号被同一个固定验证集单方面决定，五折交叉验证的意义被架空了。正确的做法是分两步——先用固定划分的验证集做早停，再用五折交叉验证只评估性能，两者不要混用。
如果XGBoost五折平均R²还在0.6以下，我会换SVR试试。带RBF核的SVR在100到200个样本的回归上经常有超预期表现，原因在于它只关注支持向量附近的样本，不会被全体样本的噪声拉偏。跑SVR之前必须标准化，核函数对尺度极其敏感，这一点和Ridge一样。SVR的好处是超参数少，调C（惩罚系数）和gamma（核宽度）就够，不容易过拟合。我在D题里用SVR做过一个独立模型，后来集成时就靠它兜底。
4. 把ADMET属性变成分类器：不平衡标签下的建模与评估细节
4.1 每个属性独立建模：阈值怎么定，标签怎么切
华为杯D题的第二部分是ADMET属性的二分类预测。ADMET是吸收、分布、代谢、排泄、毒性五个词组的缩写，每个属性都有一批化合物的实验测量值。题目要求在这些属性上做二分类，常见的做法是依据某个阈值把实验数值切成“良好/不良”两个档位，再训练分类器。关键是我第一版方案在这里走了弯路：把所有ADMET属性拼成一个多标签矩阵，试图用一个模型一次预测所有属性。
结果F1低得惨不忍睹。原因很直白：不同属性的正负样本比例差异悬殊，有的属性95%是“良好”，只有5%是“不良”。拼在一起后模型梯度被多数类主导，少数属性完全被淹没，等于白训。正确的做法是每个属性独立建模，每个属性单独定阈值、单独做特征筛选，因为同一个分子描述符对吸收和毒性的解释力完全不是一回事。
from sklearn.model_selection import StratifiedKFold from sklearn.ensemble import GradientBoostingClassifier from sklearn.metrics import f1_score, roc_auc_score results = {} for prop in ADMET_properties: X_prop = admet_feat[prop_cols[prop]] # 各属性自己的特征子集 y_prop = admet_label[prop] # 各属性自己的标签列 model = GradientBoostingClassifier( n_estimators=300, max_depth=3, random_state=42 ) skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42) f1s, aucs = [], [] for tr_idx, va_idx in skf.split(X_prop, y_prop): X_tr, X_va = X_prop.iloc[tr_idx], X_prop.iloc[va_idx] y_tr, y_va = y_prop.iloc[tr_idx], y_prop.iloc[va_idx] model.fit(X_tr, y_tr) proba = model.predict_proba(X_va)[:, 1] y_pred = (proba >= 0.5).astype(int) f1s.append(f1_score(y_va, y_pred)) aucs.append(roc_auc_score(y_va, proba)) results[prop] = [np.mean(f1s), np.mean(aucs)]这段代码的核心是StratifiedKFold。ADMET的标签分布严重不平衡时，普通KFold可能把某个属性的全部正例切到验证集里，训练过程没见过正例，模型退化成只会猜多数类的“复读机”。StratifiedKFold按类别比例分层抽样，保证每个折的类别比例和全量一致，这是不平衡分类交叉验证的基本前提。如果你的代码里还在用KFold处理二分类不平衡数据，赶快换掉。
阈值怎么切也是门学问。有些属性是连续的毒性值，有些是离散的分级。我一般先看属性的分布直方图，找到自然的拐点或题目给定的药理学标准；没有标准时用中位数切分，保证两个类别各占约50%——这是最保守的做法，不会让类别失衡到模型学不动。如果某个属性按中位数切之后某一类仍低于20%，那就是天然不平衡属性，靠后续的指标选择和阈值调整来救。
4.2 评价指标：别让准确率骗了你，阈值要单独调
第二个容易踩的坑是评价指标。某个属性90%是“良好”、10%是“不良”时，一个“全猜良好”的模型准确率直接90%，看起来很美，实际毫无价值。竞赛评审看的是你在少数类上的表现，F1和AUC才是主角。F1兼顾精确率和召回率，AUC看排序质量，两个都稳才算真稳。
predict_proba输出的是0到1的概率，默认二值化阈值是0.5。但对正例只有10%的属性，0.5通常不是最优切割点。我的做法是在验证集上遍历阈值从0.1到0.9、步长0.05，选择F1最高的阈值作为该属性的最终判定线。这个“每属性一个阈值”的细节，往往比调模型超参涨分还快。
best_thr, best_f1 = 0.5, 0.0 for thr in np.arange(0.1, 0.9, 0.05): pred_t = (proba_total >= thr).astype(int) score = f1_score(y_total, pred_t) if score > best_f1: best_f1, best_thr = score, thr print(f"最优阈值 {best_thr:.2f}，对应F1 {best_f1:.4f}")如果还想再进一步，可以做概率校准——用IsotonicRegression或Platt scaling把模型输出概率映射到真实概率分布。校准后的概率跨属性可比，你可以给化合物做综合风险评分，比如在论文里写“该化合物ADMET综合风险排名前5%”，比单看某个属性的原始概率有说服力得多。
赛道里还有一个业务导向值得写进论文：药品筛选宁可多筛掉几个候选物，也不要放过一个毒性化合物。也就是说，对毒性属性，阈值可以偏向高敏感度方向——宁愿误杀也不漏杀。这个业务逻辑会让你的阈值选择有据可依，评审看到的不再是一堆冷冰冰的0和1。
5. 华为杯D题建模避坑指南：四个让纸面分数翻车的暗坑
5.1 特征筛选在标准化之后做：数据泄露的经典剧情
现象：验证集R²高达0.92，提交后实际结果掉到0.55，评委判定模型存在严重数据泄露。百口莫辩。
原因：把全量数据喂给同一个StandardScaler做标准化，再在标准化后的全量矩阵上做特征筛选。验证集和测试集的统计量在特征选择和缩放阶段就已经暴露给模型了，验证集分数是假繁荣，真实场景中没有这个先知信息。
解决：严格按“先划分训练集和验证集→只对训练集fit标准化器→在训练集上做特征选择→用同一套规则transform验证集”的顺序执行。所有判断特征是否保留的依据，只能来自训练集。顺便说一句，特征选择也要在标准化之后做，因为LASSO这类惩罚模型对特征尺度敏感，先标准化再筛才公平。
5.2 回归特征集直接拿去做ADMET分类
现象：ADMET分类F1普遍低于0.55，调了n_estimators、learning_rate，还是上不去。
原因：回归任务筛选出的特征是冲着最大化pIC50解释力去的，ADMET各属性与活性的分子机制并不一致。一个在回归里强相关的描述符，在毒性分类里可能纯粹是噪声，用同一套特征集等于给两个任务穿同一件尺码不合的衣服。
解决：ADMET部分单独走一遍特征筛选流程，筛选依据换成对应属性的标签，而不是pIC50。哪怕最终两组特征有交集，也必须独立筛。这一步多花半小时，换来的是每个任务各自可解释，论文写起来也更扎实。
5.3 只看R²不看残差分布
现象：R²0.78看着不错，但残差图上有一条明显的二次曲线趋势——低预测值系统性偏高，高预测值系统性偏低。
原因：回归模型在取值两端拟合不足，个别高杠杆点拉高了整体R²，掩盖了局部系统性偏差。只看R²就像只看考试成绩不看错题分布，错过的全是信息。
解决：每次训练完都画残差图，横轴是预测值，纵轴是残差，看点是否围绕0均匀分布。如果两端有系统性弯折，说明模型偏置未修干净，后面用分位数校准修正分布。把这个图放进论文里，比口头说“模型表现良好”有说服力一万倍。
5.4 只交“最终答案”，忽略解答包里的复现性要求
现象：预测结果表交上去了，代码是一坨充满绝对路径和未固定随机种子的脚本，裁判换台机器重跑，结果对不上，过程也讲不清。
原因：“解答2.zip”这类提交包的评审逻辑是三维的：结果分、方法分、复现分。你可以用一个黑匣子模型拿到不错的结果分，但论文里不解释特征筛选逻辑、代码不固定随机种子，方法和复现两块分直接丢掉，得不偿失。
解决：打包前做三件事。第一，代码里所有随机种子固定，用np.random.seed(42)、random.seed(42)、模型参数里的random_state=42三管齐下。第二，参数表导出成YAML文件，谁来看都能说清每个值的来历。第三，写一个run.sh，按顺序跑通数据清洗→特征筛选→回归→分类→输出结果表，保证换台机器也能复现。这个习惯要不了多少时间，关键时刻能救人一命。
6. 让结果更稳的最后一公里：集成、分位数校准与可复现提交
6.1 模型集成与分位数校准
走到这个阶段，单模型的性能基本到顶了，继续堆迭代轮数只是浪费电。最实际的提升来自集成——把Ridge、XGBoost和SVR的预测值按验证集R²加权平均。比如Ridge验证R²0.62、XGBoost0.71、SVR0.68，权重就按三者R²占总和的比例分配。这个做法在验证集上通常比最好的单模型再高0.02到0.05，稳定且真实。分类侧同理，对概率做加权平均而不是硬投票，因为投票会丢失置信度信息，少数类的概率加权更稳。
分位数校准是我最想推荐的一个技巧，看起来有点玄学，原理却很简单：小样本回归模型的预测值分布通常方差被压缩，低值偏高、高值偏低。把验证集预测值排序取分位数，用训练集真实pIC50的对应分位数替换，能让预测分布整体对齐真实分布。这一步修的是系统性偏置，不是单点误差，所以不会给验证集涨分，但能让提交结果的分布看起来更真实可信，评审对这类的误差异普遍更宽容。
6.2 参数留痕与提交前检查
最后再说一个老生常谈但真没人爱做的事：固定随机种子、导出参数表、写复现脚本。我第一次打华为杯时，最大的教训就是模型文件命名混乱，三天后复盘时已经分不清哪个是加过特征筛选的、哪个是没加过的。现在我做任何数据集的第一版，都会强制自己在项目里建三个文件：config.yaml存全部参数、run.sh串完整流程、results/目录按日期归档结果。这个习惯救了我很多次，也让复盘变得异常轻松。2021年D题这道题，考的不只是算法能力，更是对数据处理规范性的理解——特征要不要标准化、筛选放在哪一步、早停和交叉验证怎么不打架、不平衡分类指标怎么选，每一个细节都在考你是否真的理解模型的行为边界。希望这份流程和踩坑记录能帮你在下次数模竞赛里少走几步弯路，让结果配得上你的努力。
本文还有配套的精品资源，点击获取
路由协议验证)
