尧图网站设计 尧图网站设计YAOTU DESIGN
ARTICLE DETAIL

资讯详情

深耕网站设计与一线实操的经验洞察。

CARS特征选择算法全解析:近红外光谱降维与PLS建模实战

CARS特征选择算法全解析:近红外光谱降维与PLS建模实战 简介这是一份面向光谱数据分析与机器学习特征选择场景的CARS算法Python实现适用于近红外光谱波段筛选、化学计量建模以及食品、药品无损检测等研究。CARS采用自适应重加权迭代策略在每一轮评估中利用统计量或模型预测误差衡量各波长段对目标变量的贡献动态保留重要波段、降低其权重更新滞后并逐步剔除冗余和噪声特征从而减少数据维度、增强模型稳定性与可解释性。资源以RAR压缩包提供包内包含1个Python脚本整体约2KB代码精炼便于逐行理解算法细节。目前已有2104人学习下载是一份经常被参考的特征选择实现。脚本覆盖数据预处理、权重初始化、迭代特征评估与权重更新、模型训练验证、性能比较及停止条件判断等完整环节支持替换光谱数据后直接运行输出优选波段子集与权重可辅助回归或分类建模显著降低人工筛选与调参成本。1. 先搞清楚CARS在干什么再决定要不要下载它近红外光谱数据里一个样品动辄上千个波长点但和目标成分真正相关的往往只有几十个。高维、强相关、带噪声的矩阵直接丢给PLS或随机森林模型常常把噪声当成信号学进去验证集误差大得离谱。CARSCompetitive Adaptive Reweighted Sampling竞争性自适应重加权采样是光谱特征选择里出镜率很高的算法它不投影、不压缩坐标而是直接从原始波长里挑子集保留每个波长的物理含义。我这次拆的是CARS特征选择资源包里的CARS.py是能直接改路径跑的完整实现不是demo。文章会把它的迭代过程拆开给可复现的最小脚本梳理参数边界也把互信息特征选择和递归特征消除放在一起对比。手里有近红外光谱、需要降维建模的人看完可以直接上手改。2. 拆解CARS.py自适应重加权的迭代逻辑与七个关键参数2.1 先理解CARS的反直觉机制每轮都在淘汰而不是在选优不少做特征选择的人一上来就问“怎么衡量哪个波长最重要”CARS的思路恰恰反过来。它不是在每轮给所有波长打分然后挑前几名而是先默认所有波长都参与建模然后找出当前模型里最不重要的那批变量把它们砍掉下一轮在剩余变量上重新建模、重新排序。真正留下来的特征是在一轮轮淘汰中“幸存”的而不是一次排序选出来的。这个机制由三层结构组成。第一层是蒙特卡洛采样每一轮只随机抽取一部分样本去拟合PLS模型模拟不同样本组合下的系数变化避免某一组特殊样本主导整个选择结果。第二层是删减比例控制常见实现用指数衰减函数EDF决定每轮删多少变量前期变量冗余多删得快后期剩下的大多是有效信息删得慢防止一口气把重要波段误删。第三层是自适应重加权用当前PLS回归系数的绝对值作为每个波段的权重下一轮只对保留波段重新计算权重随迭代更新。三层叠加到一起变量数从一千多逐步收敛到几十个同时每一轮都留一组交叉验证误差用来跟踪效果。EDF常被当成玄学其实它只干一件事算出“第i轮应该删除当前保留变量的百分比”。常用形式是p_i a * exp(-k * i)其中a由初始保留比例和最低保留比例反推k控制衰减快慢。CARS.py里一般先把总变量数和目标最少变量数代入算出衰减斜率再在每一轮用这个比例乘以当前保留数得到本轮删除数量。这样能保证前几轮砍得多、后几轮砍得少不会出现迭代到一半变量已经被清空的情况。2.2 核心循环的代码骨架主循环里到底发生了什么下面这段骨架是我按CARS.py里最常见的实现思路整理的把主循环拆成五步采样、拟合、淘汰、交叉验证、记录最优。拿到原始文件后对照这段代码很容易定位每一块在做什么。import numpy as np from sklearn.cross_decomposition import PLSRegression from sklearn.model_selection import KFold from sklearn.metrics import mean_squared_error def cars_select(X, y, n_mc50, sample_ratio0.9, n_components10, drop_ratio_start0.1, min_keep5, random_state42): n_samples, n_vars X.shape keep np.arange(n_vars) # 当前保留的原始波长索引 best_rmssecv np.inf best_keep keep.copy() for mc in range(n_mc): # 1) 蒙特卡洛采样每轮只用 90% 样本拟合 PLS rng np.random.default_rng(random_state mc) idx rng.choice(n_samples, int(n_samples * sample_ratio), replaceFalse) Xs, ys X[idx], y[idx] # 2) 在当前保留变量上拟合 PLS系数绝对值作为权重 pls PLSRegression(n_componentsn_components, scaleTrue) pls.fit(Xs[:, keep], ys) coef np.abs(pls.coef_.ravel()) # 3) 按权重从小到大排序淘汰最不重要的 drop_ratio 比例 # 这里用固定比例演示真实 CARS.py 通常用 EDF 动态计算 n_drop max(1, int(len(keep) * drop_ratio_start)) order np.argsort(coef) kill keep[order[:n_drop]] keep_new np.setdiff1d(keep, kill) if len(keep_new) min_keep: break # 4) 5 折交叉验证粗略评估当前子集的 RMSECV kf KFold(n_splits5, shuffleTrue, random_staterandom_state) errs [] for tr_idx, va_idx in kf.split(X[:, keep_new]): m PLSRegression(n_componentsn_components, scaleTrue) m.fit(X[tr_idx][:, keep_new], y[tr_idx]) pred m.predict(X[va_idx][:, keep_new]).ravel() errs.append(mean_squared_error(y[va_idx], pred)) rmssecv np.sqrt(np.mean(errs)) # 5) 当前子集误差更低就记录下来作为最优特征子集 if rmssecv best_rmssecv: best_rmssecv rmssecv best_keep keep_new.copy() keep keep_new return best_keep, best_rmssecv第二步是整个算法的权重来源。PLS回归系数本身表示“这个波长的变化对目标值的影响方向和强度”系数绝对值大说明该波段与浓度或性质的协变关系强应该保留绝对值接近零的要么是纯噪声要么它的信息已经被其它波长重复表达删掉损失不大。第三步要特别强调真实CARS.py里很少用固定drop_ratio_start而是用EDF每轮动态计算删除比例。上面代码用固定比例是为了让逻辑好跟踪但你把固定0.1套到自己的数据上时前几轮可能砍得太狠把弱但有用的波段提前删掉。拿到脚本后先确认它有没有EDF参数有就用动态比例。第四步的KFold交叉验证是评估子集质量的关键。这里设5折是折中方案样本只有几十个的小数据集改3折更稳妥样本到几百个时可以上7折。折数越多RMSECV估计越准但耗时成倍上涨而且样本量小的时候折数太多反而让每折训练集过短误差估计不稳定。2.3 数据进CARS之前平滑、散射校正和矩阵形状检查CARS.py接收的X是二维数组形状是(样本数, 波长点数)y是一维浓度或性质数组。我经手的近红外数据里几乎没有拿原始吸光度直接跑就效果好的至少要先做两步处理。第一步是平滑去噪。最常用的是Savitzky-Golay平滑窗口一般取5~15个点多项式阶数1~3。窗口太小去不干净噪声太大会把相邻波段的细微差异抹平我一般从窗口7开始试用RMSECV来选。第二步是散射校正。固体粉末样品的颗粒度、装样松紧会造成光谱基线平移和斜率变化常用的有标准正态变换SNV和多元散射校正MSC。我的经验是固体样品先试SNV液体样品先试MSC两个都跑一遍取验证集误差小的路径。容易踩的一个坑是输入格式。CARS迭代里频繁做X[:, keep]这种切片要求X保持二维结构千万别在预处理环节把光谱压成一维。另外如果你的y是多列比如同时预测水分、蛋白、脂肪一定要分开跑CARS因为不同目标的响应波长不一样硬塞在一起选出的特征往往两头不讨好。调用前我会做一次格式检查from sklearn.preprocessing import StandardScaler assert X.ndim 2 and X.dtype in (np.float64, np.float32), X 必须是二维浮点矩阵 X np.nan_to_num(X, nannp.nanmean(X)) # 每个波长单独标准化把系数权重拉到同一尺度 scaler StandardScaler() X_scaled scaler.fit_transform(X) best_idx, best_rmssecv cars_select(X_scaled, y, n_mc50) print(选中的波长索引:, best_idx) print(最优 RMSECV:, best_rmssecv)这里有个细节值得说StandardScaler是对每一列波长做标准化它不会改变波段间的相关性但会把回归系数的尺度拉齐。如果光谱基线漂移明显换成SNV再跑通常比单纯StandardScaler更对症因为StandardScaler只按列处理SNV按样本行整体校正对颗粒散射的影响处理得更到位。实际项目里我习惯两条预处理路径都跑看最终交叉验证误差再定。2.4 七个关键参数怎么设从n_mc到random_stateCARS.py里最常需要调的参数有七个按影响程度排序如下表。参数常见范围设置建议n_mc50 ~ 500蒙特卡洛迭代次数初步探索取100精修取200sample_ratio0.7 ~ 0.95每轮采样比例样本量少于100时取0.9n_components3 ~ 20PLS主成分数用交叉验证确定不要拍脑袋EDF初值 / 删减比例0.05 ~ 0.2优先用EDF动态控制避免固定比例误删min_keep3 ~ 20保留变量下限防止迭代收敛到空集折数3 ~ 10交叉验证折数小样本用3大样本用5或7random_state任意整数固定随机种子否则结果不可复现参数这块是被忽略最多的部分。CARS的表现对n_components和random_state尤其敏感。n_components设得过大PLS会在高方差方向上过拟合CARS跟着把噪声波段当成高贡献波段留下来设得太小模型欠拟合重要波段的系数表达不出来选出的子集没有意义。我自己的习惯是先固定n_mc100对3到15个主成分做一轮网格搜索找到RMSECV最低点再把那个主成分数带回CARS正式运行。随机种子的影响也很大。因为蒙特卡洛采样本身就是随机的不固定种子的话同一份数据跑两次选出的波段重合率可能只有70%。初学的人会觉得算法不稳定其实只是种子没锁住。正式实验前先固定random_state每次运行的结果单独保存最后做交集分析这才是CARS的正确使用姿势。2.5 样本量特别小怎么办如果你的样本量只有三四十个按上面默认参数跑会明显感觉RMSECV曲线抖动厉害。我的处理方式是把sample_ratio提到0.95每轮只丢两三个样本同时把折数降到3增加每折训练集大小。另一个补救是加大n_mc到200以上让随机采样覆盖更多样本组合再用多次运行的波段频率做最终筛选。样本量小的时候不要迷信单次运行的最优子集复现性比单次误差更重要。3. 把CARS跑起来数据输入、调用方式与输出解读3.1 从csv到CARS.py的最小运行流程拿到CARS.rar后我做的第一件事是检查CARS.py末尾的入口。这类资源最常见的写法是文件末尾带一个if __name__ __main__块里面写死数据路径和处理参数你需要做的就是把路径换成自己的数据。如果你的光谱是csv用下面这种方式接入import numpy as np import pandas as pd # 约定第一列是样品编号最后一列是浓度中间全是波长 df pd.read_csv(near_ir_data.csv) wavelengths df.columns[1:-1].astype(float) X df.iloc[:, 1:-1].values y df.iloc[:, -1].values # 导入 CARS.py 中的选择函数 from cars import cars_select best_idx, score cars_select(X, y, n_mc100, n_components8, random_state42) print(保留波段序号:, best_idx) print(对应波长:, wavelengths[best_idx]) np.save(selected_idx.npy, best_idx)pandas读进来之后先检查有没有NaN和字符串列。光谱数据偶尔会有表头塞单位、列名带空格的情况pandas会把它们读成object类型直接传给CARS会报类型错误。我习惯在得到X之后立刻打印X.shape和X.dtype确认是float二维矩阵再往下走。还有一个容易被忽略的点波长列如果是整数型的nm值astype(float)没问题但如果仪器导出的波长带小数点后多位记得保留原始精度不要随意round否则最后输出波段位置会和光谱图对不上。3.2 输出文件里能看到什么CARS.py跑完通常输出三类东西选中的波段索引常见格式是npy或txt每轮RMSECV轨迹一般是数组或列表最终PLS模型系数用来反查哪些波段贡献最大。如果脚本只返回一组索引我强烈建议自己把RMSECV轨迹也存下来因为它是判断收敛的关键。import matplotlib.pyplot as plt fn np.load(cars_trace.npy) # 假设原始脚本把每轮 RMSECV 存成了 npy plt.plot(range(len(fn)), fn) plt.xlabel(iter) plt.ylabel(RMSECV) plt.savefig(cars_trace.png, dpi150)看这条曲线有个经验误差先下降后上升说明算法在“删冗余信息”和“开始删有效信息”之间找到了分界点最优子集就在拐点附近。如果从第一轮到最后一轮都在下降说明变量删得还不够应该加大n_mc或调低保留下限。如果曲线全程平坦多半是n_components设得太小模型本身没拟合起来CARS删什么都差不多这时候回头调主成分数比调CARS参数更有用。3.3 选完特征之后用交叉验证比较全光谱和子集模型CARS的价值不是“选出一组波长编号”而是给后续建模提供一个更干净的输入。我的固定步骤是CARS选出特征后用同样的PLS配置在选出的子集上建模比较全光谱模型和特征子集模型的验证集R²与RMSECV这一步能直观看到降维带来的收益。def eval_pls(X, y, n_components8): from sklearn.cross_decomposition import PLSRegression from sklearn.model_selection import cross_val_predict pls PLSRegression(n_componentsn_components, scaleTrue) y_pred cross_val_predict(pls, X, y, cv5) r2 1 - np.sum((y - y_pred) ** 2) / np.sum((y - np.mean(y)) ** 2) rmsecv np.sqrt(np.mean((y - y_pred) ** 2)) return r2, rmsecv r2_full, rmsecv_full eval_pls(X_scaled, y) r2_cars, rmsecv_cars eval_pls(X_scaled[:, best_idx], y) print(ffull spectrum: R2{r2_full:.3f}, RMSECV{rmsecv_full:.3f}) print(fcars selected: R2{r2_cars:.3f}, RMSECV{rmsecv_cars:.3f})这里的cross_val_predict用的是5折外样本预测比在训练集上重拟合更接近真实泛化表现。选出的特征如果能让RMSECV明显下降同时变量数从一千降到几十这步就算成功。如果误差持平而变量数大幅下降同样算成功因为模型更简单、抗噪声能力更强。如果误差明显变差优先怀疑三个地方折数太高导致训练样本过少PLS主成分数对子集来说偏小或者CARS的删减比例太激进把关键波段误删了。比较完之后我一般会把选中的波长序号和对应波长值存成csv顺便记录预处理参数、CARS参数和本次RMSECV。光谱分析项目经常要追溯“为什么选这组波段”有一份完整参数记录比事后翻代码快得多。4. CARS常见问题与避坑记录从运行报错到结果不可复现4.1 现象同一份数据跑两次选出的波段不一致这是我见过最多的“CARS不靠谱”反馈。原因不复杂CARS的蒙特卡洛采样是随机的每轮随机抽90%的样本抽到的组合不同PLS系数就有波动淘汰路径自然不同。解决分两步。第一固定random_state第二正式结果跑5到10次统计每个波段被选中的频率只留下出现频率高于60%的波段作为高重现波段。from collections import Counter freq Counter() for seed in range(10): idx, _ cars_select(X_scaled, y, random_stateseed) freq.update(idx.tolist()) high_freq [k for k, v in freq.items() if v 6] print(10 次运行中至少出现 6 次的波段数:, len(high_freq))固定种子之后结果可复现但新手要明白一个边界固定种子只是消除随机性不代表选出的波段就是唯一最优解。不同种子下模型的RMSECV可能有差异这正是CARS的随机采样特性决定的。如果把它当黑匣子用忽略种子和迭代次数的影响很容易得出误导性结论。4.2 现象特征越选越少RMSECV反而上涨这个现象通常由两个原因叠加。第一个是PLS主成分数设得太高模型在训练集上过拟合CARS把噪声波段当成高贡献波段留下真正有效的波段反而被挤掉。第二个是删减比例设得激进比如把固定删减比例调到0.2重要波段还没来得及被重新加权就已经被误删后续迭代里它再也没有机会回来。我的解决顺序是先用3到15个主成分做网格搜索找到交叉验证误差最低的n_components再把删减比例调到0.05~0.08或者改用EDF动态比例让删减过程温和一些。还有一个容易被忽略的机制CARS每轮的权重是基于当前保留变量重新计算的删掉一个变量之后其它变量的权重会重新分配。因此前面误删的影响会被后续迭代放大这也是为什么CARS不能只看最终输出必须配合RMSECV轨迹一起判断。如果误差曲线在某个位置开始掉头向上说明从那之后删掉的波段已经开始伤害模型了。4.3 现象n_mc500跑了一个小时还没结束CARS每轮都要做一次PLS拟合和交叉验证n_mc大、折数多、样本量大时总耗时非常可观。这不是死机是算法本身就是这样设计的。省时间的做法是分两阶段先用n_mc50粗跑看RMSECV轨迹有没有出现拐点确认拐点位置后再用n_mc200在拐点附近精修。粗跑阶段把折数降到3精修阶段再上5折能省一半以上的时间。另外n_components不要设到20以上。近红外数据通常是高维小样本20个主成分的协方差矩阵估计非常慢而且几乎必然过拟合。我的上限是15多数数据集在8~12之间就能达到平稳误差。4.4 现象CARS选出的波段落在一个很窄的区间看起来“太集中”近红外光谱的波长是连续的真正相关的官能团区域往往就是一两个小区间所以波段集中本身不是问题。反而要警惕的是选出的波段散布在完全无关的区域比如水分特征峰和蛋白质特征峰同时被保留那说明X和y的关联不够单一或者预处理没有把散射基线处理干净建议先做SNV或MSC再跑一次。另一个更常见的情况是保留的波段数量仍然偏多比如一千个波长选完还剩两百个。这通常是min_keep设得偏高或者RMSECV曲线的最低点比拐点滞后太多。我的做法是在RMSECV轨迹上找“拐点”而不是“最低点”作为子集规模。最低点往往已经过拟合拐点对应的子集更保守跨批次数据的重现性也更好。5. CARS与其它特征选择方法的边界互信息特征选择和递归特征消除5.1 三者逻辑差异模型系数、信息量、递归训练很多教程把特征选择讲成工具箱里随便拿一把其实CARS、互信息特征选择和递归特征消除背后的逻辑完全不同。CARS属于“模型系数加权迭代删减”核心依赖PLS回归系数判断波长重要程度选出的特征直接服务于后续的PLS建模。互信息特征选择从信息论角度出发计算每个波长与目标y之间的互信息衡量的是“知道这个波长后对y的不确定性减少多少”它不依赖具体预测模型所以不受模型过拟合的影响。递归特征消除RFE则是每轮训练一个模型、删掉权重最低的特征、在剩余特征上递归重复。它和CARS看起来像实际上是两回事RFE用全部样本训练CARS故意只抽部分样本制造多样性RFE对任何有coef_或feature_importances_的模型都通用CARS则天然为光谱线性模型设计。这三者不是对立关系而是不同风险偏好。CARS的优势是模型相关性强选出来的波段直接服务PLS互信息的优点是单变量独立性不受模型容错机制影响适合冷启动RFE的优势是通用性强不限定光谱场景。5.2 参数和适用场景对比方法核心实现主要参数适合场景CARSPLS系数 蒙特卡洛 EDFn_mc、n_components、删减比例近红外光谱加PLS建模互信息特征选择单变量互信息估计n_neighbors、选择阈值冷启动快速粗筛、检查单变量相关RFE递归训练删特征估计器、n_features_to_select中小维度表格数据、非线性模型互信息门槛最低sklearn一行就能算出来from sklearn.feature_selection import mutual_info_regression mi mutual_info_regression(X_scaled, y, n_neighbors5, random_state42) top_mi np.argsort(mi)[-50:]n_neighbors是连续变量互信息估计里的平滑参数影响非常大。5是折中值样本量上千可以试3噪声明显时试10。互信息分数不区分方向和交互作用也就是说两个变量单独看都不显著、合并起来却很强的组合单变量互信息会漏掉。这是所有单变量筛选方法的通病。递归特征消除的代码套路是这样的from sklearn.feature_selection import RFE from sklearn.linear_model import Lasso estimator Lasso(alpha0.01) rfe RFE(estimator, n_features_to_select30, step0.1) rfe.fit(X_scaled, y) print(RFE selected count:, np.sum(rfe.support_))RFE的step参数控制每轮删除的比例0.1表示每轮删10%当前变量。它比CARS好理解但计算量更大因为每轮都要在全量样本上重新训练Lasso。对上千波长的光谱数据RFE跑起来明显比CARS慢而且Lasso的alpha需要额外调。5.3 推荐组合CARS粗筛、互信息精排我最近做的一个近红外数据集就是这么处理的1700个波长先跑CARS选了70多个再用互信息从这70个里精排到30个最终PLS的RMSECV比单独用CARS或单独用互信息都低。原因是CARS保留波段间协同信息而互信息能把CARS里“模型勉强拟合出来”的冗余项再压掉一层。from sklearn.feature_selection import SelectKBest, mutual_info_regression selector SelectKBest(mutual_info_regression, k30) X_selected selector.fit_transform(X_scaled[:, best_idx], y) final_columns best_idx[selector.get_support()]这段代码的意思很直接在CARS选出的best_idx子集上再按互信息排序挑前30个波段。两段式比单跑CARS多一道工序换来的是子集更紧凑、跨批次重现性更高。代价是超参数变成两组调参成本高一截所以更适合要做长期稳定模型的课题或项目而不是一次性出结果的临时分析。如果做的是判别任务比如近红外真药假药分类CARS同样能用只是内部模型从PLSRegression换成带coef_的分类器比如PLSCanonical或逻辑回归。CARS.py里如果写死PLSRegression需要自己动手改一行接口差异不大。6. 验证CARS结果的三板斧RMSECV轨迹、多次运行交集、波段解释拿到CARS选出的波段别急着写进报告先用三个方法交叉验证。第一看RMSECV轨迹是否出现“先降后升”的拐点如果有最优子集在拐点附近而不是终点。第二用5到10个不同随机种子分别跑CARS统计每个波段被选中的频率高频波段大概率对应真实化学信息偶尔被选中的可能是噪声或随机扰动。第三把选出的波长对照近红外官能团归属表水分在1400~1450nm和1900~1950nm附近有强吸收蛋白质的N-H吸收在2050~2100nm附近。如果选出的波段完全落在这类已知区域之外就要检查预处理或参数是不是出了问题。一个很实用的小技巧是画“选中频率-波长”叠加图把多次运行的选择频率柱状图直接叠在平均光谱上。频率峰和光谱峰重合的位置往往就是最有解释力的波段比单独看一组索引直观得多。freq_arr np.zeros(X.shape[1]) for seed in range(10): idx, _ cars_select(X_scaled, y, random_stateseed) freq_arr[idx] 1 plt.bar(wavelengths, freq_arr, width2, alpha0.5) plt.plot(wavelengths, X_scaled.mean(axis0), colork, linewidth0.8) plt.xlabel(wavelength (nm)) plt.ylabel(selection frequency) plt.savefig(car_frequency.png, dpi150)三种验证都通过之后我会把选中的波段、预处理参数、CARS参数和验证指标整理成一份记录。从那以后每次拿到新的光谱数据凡是准备正式使用的特征子集我都会强制走一遍这三步验证先看轨迹、再看交集、最后对照波段含义缺一步都不敢直接拿去建模型。希望帮到你。本文还有配套的精品资源点击获取
返回列表