制剂 CQA 预测模型开发教程13支持向量回归与高斯过程回归——小样本建模与不确定度量化版本声明块工具/软件Python 3.10.11 scikit-learn 1.7.2 numpy 2.2.6 scipy 1.15.3数据NIR Shootout 2002 片剂数据集655 片 × 650 波长本文目标读完后你能说清 SVR 的 C 与 ε 各管什么、GPR 的核函数三项各管什么并能自己跑出预测值 ± 不确定度的预测带、核对它的区间覆盖率是否名副其实。一句话结论在 NIR Shootout 2002 的 assay 数据SNV 预处理、校正集 155 片上支持向量回归SVRrbf 核、C100、epsilon1.0取得 R²cal 0.9426、R²test 0.9144、RMSEP 4.6074、RPD 3.42是第 12 篇六模型对比中唯一在测试集上超过偏最小二乘回归PLS的模型高斯过程回归GPR核ConstantKernel(1.0) * RBF(length_scale50.0) WhiteKernel(1.0)、normalize_yTrue、alpha1e-6取得 R²test 0.9064、RMSEP 4.8181其预测标准差均值为 5.0287而 95% 置信区间的实测覆盖率为 0.9761——高于名义的 0.95说明区间略过宽。〇、本篇要解决的认知问题SVR 里的C与epsilon都是超参数它们各自控制模型的哪一部分为什么 ε 的量纲必须与 CQA 相同核技巧到底技巧在哪为什么 rbf 核能在只有 155 个校正样本时把非线性关系建起来高斯过程回归的核写成ConstantKernel * RBF WhiteKernel这三项分别管什么少一项会怎样normalize_yTrue到底改了什么为什么 assay 这种 y 波动很大的 CQA 必须打开它GPR 用return_stdTrue给出的预测标准差能直接当成这个预测值可信程度来报吗怎么核对一、机制解析1.1 SVR先容忍 ε再用 C 惩罚越界SVR支持向量回归Support Vector Regression与普通最小二乘最大的差别是损失函数。最小二乘对每一个残差都算平方并求和模型必须对每个点都讨好SVR 则用ε-不敏感损失ε-insensitive lossL_ε(r) max(0, |r| − ε) # r y − ŷ |r| ≤ ε → 损失为 0完全不管 |r| ε → 损失随超出部分线性增长用图表示就是一条ε 管道y │ ← 上边界ŷ ε │ ● ● │ ●────────● ← 管道内部损失 0这些点不进模型 │ ● ● │ ← 下边界ŷ − ε │ ● ○ ← 管道外的点成为支持向量被 C 惩罚 └────────────────────────► x由此得到三句话ε 决定多准才不算错ε 是管道半宽量纲必须与 y 相同。落到管道内的样本损失为 0不会出现在最终模型里——这正是 SVR 能稀疏的原因。C 决定越界得多疼C 是惩罚系数等价于正则化的倒数。C 越大越不允许样本落在管道外模型被迫更贴合训练数据容易过拟合C 越小模型越平偏差换方差。只有管道外的样本支持向量决定最终解因此 SVR 的解对噪声有一定免疫——这是它在 155 个小样本上的天然优势。把这两个参数放到本系列的数据上看量级关系立刻清楚参数SPEC 实测取值与 assay 尺度的关系直觉判断epsilon1.0 mg校正集 assay 的 SD 是 21.9803 mgε ≈ 0.045 个 SD管道很窄只忽略很小的误差C100.0相对 γ 的默认量级gammascale属于强惩罚允许模型有较强的局部弯曲经验法则ε 的第一个候选值应当取仪器重复性噪声的水平例如同一片药重复测 3 次的 RMSE而不是凭感觉填。ε 填成 y 的 SD 量级这里 20 左右管道会吞掉全部残差模型退化成一条水平线。1.2 核技巧不显式升维只算内积SVR 的原始形式是线性的在 650 维波长空间里找一个线性函数。核技巧要做的事是把样本映射到更高维空间 φ(x)在那里线性可分——但绝不真的去算 φ(x)只算内积K(x_i, x_j) 〈φ(x_i), φ(x_j)〉 核函数直接给出内积跳过 φ rbfK(x_i, x_j) exp(−γ·‖x_i − x_j‖²) γ 控制相关半径rbf 核的物理含义很直观两个光谱越接近它们的核值越接近 1越远则衰减到 0。于是模型可以写成新样本的预测值 各支持向量的相似度加权和——这正是非线性的来源。gammascalescikit-learn 默认取1 / (n_features × X.var())。本数据 SNV 后 X 的方差约为 1、n_features 650故 γ 约为 1/650 量级。不要手动把 γ 设成 1 之类的整数好看值那会让核的衰减速度提高几十倍模型退化成局部插值。1.3 GPR把函数本身当成随机变量高斯过程回归GPRGaussian Process Regression换了一个完全不同的视角它不问哪个函数最好而假设所有可能的函数服从一个高斯过程先验f(x) ~ GP( m(x), k(x, x) ) ↑ 均值函数 ↑ 协方差函数核观测模型是y f(x) εε 是测量噪声。给定训练集后用贝叶斯公式求后验于是对任一新样本 x*预测结果不是一个数而是一个正态分布p(y* | x*, 训练集) N( μ(x*), σ²(x*) ) ↑ 预测均值 ↑ 预测方差 ← 这就是 GPR 的独特价值对比第 12 篇的 SVRSVR 只给点估计GPR 顺手把不确定度也给了而且这个 σ² 不是事后拟合出来的是贝叶斯推导的直接产物。1.4 核函数三项分解谁管幅度、谁管形状、谁管噪声SPEC 实测的核长这样ConstantKernel(1.0) * RBF(length_scale50.0) WhiteKernel(1.0) └──────┬──────┘ └────────┬────────┘ └─────┬─────┘ ①信号幅度 ②信号形状 ③观测噪声组成数学角色管什么去掉它的后果ConstantKernel乘性常数 c信号方差幅度²决定预测值能摆多大核的幅度被钉死为 1无法匹配 y 的波动尺度RBF(length_scale50.0)平方指数核空间中的相关长度越大越平滑没有它就没有平滑假设GPR 退化为纯噪声模型WhiteKernel对角噪声项 δ(ij)独立同分布观测噪声mg² 量纲模型会强行穿过每个训练点区间变得极窄且过度自信三者是加法乘法的关系不是随便组合乘号把幅度与形状绑成一个信号核加号把信号与噪声分开。这是 GPR 核设计里最常用的标准写法。重要复现口径SPEC 给出的length_scale50.0是核的初始值不是最终拟合值。scikit-learn 的GaussianProcessRegressor默认optimizerfmin_l_bfgs_b拟合时会继续优化核超参数。本系列实测scikit-learn 1.7.2优化后的核为6.86**2 * RBF(length_scale35.4) WhiteKernel(noise_level0.0365)。如果为了冻结参数把optimizer设为None模型会退化得完全不能用——这正是本篇第一个反直觉点详见第三节。1.5normalize_yTrue让初值与 α 抖动项不再依赖 y 的量纲normalize_yTrue会在拟合前把 y 减去均值、除以标准差拟合完再变换回来。它改变的是优化问题的尺度normalize_yFalse核超参数与 alpha 的初值必须猜到 y 的量级上 normalize_yTrue y 被压到均值 0、方差 1alpha1e-6 成为与 CQA 无关的无量纲抖动在本数据上这一点尤其关键assay 的校正集 SD 是 21.9803 mg而 weight 只有 5.5842 mg。同一组alpha1e-6、length_scale50.0的初值在 weight 上和在 assay 上代表完全不同的正则强度。打开normalize_yTrue后超参数初值只与归一化后的尺度有关同一套配置才有跨 CQA 复用的可能。alpha1e-6是加在核矩阵对角线上的数值抖动项jitter作用是让K αI可逆不是噪声水平——噪声由WhiteKernel负责。这两个容易搞混。1.6 只有 155 个校正样本时该选谁把 SPEC 实测的两组结果并列assaySNV模型关键配置R²calRPDcalR²testRMSEPRPDtest是否给不确定度SVRrbf, C100, ε1.00.94264.190.91444.60743.42否GPRrbf 核三项normalize_yTrue——0.90644.8181—是σ 均值 5.0287PLS第 12 篇对照nLV50.96165.120.89795.03003.13否结论可以写得很具体只要点估计精度选 SVR它的 R²test 比 GPR 高 0.0080比 PLS 高 0.0165是三者最优。要不确定度选 GPR它多付出的代价是 R²test 低约 0.008换来的是每个预测值都附一个方差。对第 19 篇要讲的实时放行RTRT场景这个预测值离规格线还有几个 σ比R² 高 0.008重要得多。注意 SVR 的 R²cal0.9426反而低于 PLS0.9616这不是缺点而是健康模型的特征——校正集不追求满分泛化不退化。二、完整代码与逐行剖析2.1 SVR复现 R²test 0.9144第 13 篇 · 2.1SVR 在 assay 上的实现SNV 预处理 rbf 核importnumpyasnpimportpandasaspdfromscipy.ioimportloadmatfromsklearn.svmimportSVRfromsklearn.metricsimportr2_score,root_mean_squared_errordefload_shootout(path):解包 MATLAB DataSet Object返回三集光谱、三个 CQA 与波长轴。mloadmat(path)defunpack(name):# 外层是 (1,1) 结构体必须 [data][0, 0] 逐层解包arrm[name][data][0,0]# 原始 dtype 为大端 f8显式转 float64否则矩阵运算类型不匹配returnnp.asarray(arr,dtypenp.float64)defwave_of(name):# 波长轴不在顶层而是嵌在每个 DataSet Object 的 axisscale 字段里foriteminm[name][axisscale][0,0].ravel():arrnp.asarray(item)ifarr.size1andarr.dtype.kindiniuf:# 挑出数值型刻度数组returnarr.astype(np.float64).ravel()raiseValueError(未找到波长刻度)return{Xcal:unpack(calibrate_1),Xval:unpack(validate_1),Xtest:unpack(test_1),Ycal:unpack(calibrate_Y),Yval:unpack(validate_Y),Ytest:unpack(test_Y),wave:wave_of(calibrate_1),# 600 → 1898 nm}defsnv(X):逐样本标准正态变量变换铁律 10不许用全局均值/方差标准化。return(X-X.mean(axis1,keepdimsTrue))/X.std(axis1,ddof1,keepdimsTrue)defreport(tag,y_true,y_pred):统一的单集评价RMSE 走铁律 1 的 root_mean_squared_error。y_truenp.asarray(y_true,float).ravel()y_prednp.asarray(y_pred,float).ravel()rmseroot_mean_squared_error(y_true,y_pred)return{set:tag,R2:r2_score(y_true,y_pred),RMSE:rmse,RPD:float(np.std(y_true,ddof1)/rmse),bias:float(np.mean(y_pred-y_true))}if__name____main__:dload_shootout(nir_shootout_2002.mat)Xcal,Xtestsnv(d[Xcal]),snv(d[Xtest])ycal,ytestd[Ycal][:,2],d[Ytest][:,2]# 第 3 列 assaymg# gamma 不显式指定 → 取默认 scale 1/(n_features · X.var())与本系列实测口径一致svrSVR(kernelrbf,C100.0,epsilon1.0)svr.fit(Xcal,ycal)rows[report(cal,ycal,svr.predict(Xcal)),report(test,ytest,svr.predict(Xtest))]print(pd.DataFrame(rows).round(4).to_string(indexFalse))print(支持向量数,int(svr.support_.size),/,ycal.size)# 稀疏性直接可查实测输出与 SPEC 5.10 逐位一致setR²RMSE (mg)RPDbias (mg)cal0.94265.25094.19−0.2895test0.91444.60743.420.5442一个必须如实写出来的细节本模型的支持向量数是129 / 155。也就是说epsilon1.0 mg这条管道对 assay 的噪声水平而言偏窄绝大多数样本都落在管道外、都成了支持向量SVR 的稀疏性在此几乎没体现出来。这不影响它的预测精度R²test 0.9144 是实测的但意味着两件事①模型不能靠支持向量少来宣称轻量②如果把 ε 放大支持向量数会显著下降同时精度会先稳后降——ε 是稀疏度与拟合紧度之间的旋钮不是越小越好。判断方式很简单svr.n_support_.sum()越接近样本数模型越接近逐点加权而非稀疏核展开。2.2 GPR用return_stdTrue画预测带并核对覆盖率第 13 篇 · 2.2GPR 预测带 95% 区间覆盖率核对 接 2.1 的主流程继续执行沿用其中的 Xcal / ycal / Xtest / ytest 与 reportimportnumpyasnpimportmatplotlib.pyplotaspltfromsklearn.gaussian_processimportGaussianProcessRegressorfromsklearn.gaussian_process.kernelsimportConstantKernel,RBF,WhiteKernel# ① 核幅度 × 平滑形状 观测噪声三项缺一不可见 1.4kernelConstantKernel(1.0)*RBF(length_scale50.0)WhiteKernel(1.0)gprGaussianProcessRegressor(kernelkernel,normalize_yTrue,# 关键先把 y 归一到零均值单位方差超参数初值才不依赖 CQA 量纲alpha1e-6,# 核矩阵对角抖动仅保证可逆不承担噪声建模random_state42,# 铁律 9固定随机种子)# optimizer 用默认值 fmin_l_bfgs_b必须让优化器继续优化核超参数gpr.fit(Xcal,ycal)# ② return_stdTrue → 同时拿到预测均值与预测标准差形状同为 (n,)mu_cal,sd_calgpr.predict(Xcal,return_stdTrue)mu_test,sd_testgpr.predict(Xtest,return_stdTrue)# ③ 用 1.96σ 构造 95% 区间核对名义覆盖率 0.95 是否名副其实lomu_test-1.96*sd_test himu_test1.96*sd_test coveragefloat(np.mean((ytestlo)(ytesthi)))# 真实值落进区间的比例print(f优化后的核{gpr.kernel_})print(fR2test {r2_score(ytest,mu_test):.4f})print(fRMSEP {root_mean_squared_error(ytest,mu_test):.4f})print(f预测标准差均值 {sd_test.mean():.4f}mg)print(f95% 区间实测覆盖率 {coverage:.4f}名义 0.95)实测输出与 SPEC 5.11 一致优化后的核6.86**2 * RBF(length_scale35.4) WhiteKernel(noise_level0.0365) R2test 0.9064 RMSEP 4.8181 预测标准差均值 5.0287 mg 95% 区间实测覆盖率 0.9761 名义 0.952.3 把预测带画出来第 13 篇 · 2.3预测带可视化沿用 2.2 的 mu_test / sd_test / lo / hi / ytest# 横轴按参考值排序避免测试集无序导致的锯齿形来回穿越ordernp.argsort(ytest)fig,axplt.subplots(figsize(8,4.5))ax.fill_between(np.arange(ytest.size),lo[order],hi[order],alpha0.25,label95% 预测区间)ax.plot(np.arange(ytest.size),mu_test[order],lw1.2,labelGPR 预测均值)ax.scatter(np.arange(ytest.size),ytest[order],s4,ck,labelassay 参考值)ax.set_xlabel(测试样本按参考值排序);ax.set_ylabel(assay / mg)ax.legend();plt.tight_layout();plt.show()看这张图只需要确认两件事① 黑色的参考点是否绝大多数落在带内本系列实测覆盖率 0.9761应看到极少跑出带外的点② 带的宽度是否随样本变化——若整条带几乎等宽说明WhiteKernel的噪声项主导了方差、RBF项几乎没有随输入变化此时核结构设计需要检讨。2.4 怎么读 0.9761 这个覆盖率覆盖率 0.9761 0.95 的含义必须读准覆盖率含义工程解读≈ 0.95区间宽度与真实误差匹配理想0.9761本系列实测区间偏宽把不该包住的样本也包住了保守但安全不会漏报风险代价是区间偏大、判别力偏低明显 0.95区间偏窄过度自信危险真实值频繁落在区间外不可用于放行判定对制药场景略过宽是可接受的方向过窄才是事故。0.9761 说明这套核超参数下的 GPR 倾向于保守与预测标准差均值 5.0287 mg 相对 RMSEP 4.8181 偏大是同一件事——σ 的均值大于实际误差的均方根区间自然包得更多。2.5 反直觉的默认值项直觉实际行为GPR.optimizer“我都把核写死了它就该用我给的数”默认fmin_l_bfgs_b会继续优化SPEC 的length_scale50.0只是初值。设成None会让核停在未优化的初值上模型直接崩掉GPR.alpha是噪声水平是加在对角线上的数值抖动真正的观测噪声应由WhiteKernel表示normalize_y归一化只是预处理它决定了alpha与核超参数初值的有效尺度关掉后同一套初值在 weightSD 5.5842与 assaySD 21.9803上含义完全不同SVR.gamma越大越好默认scale已按特征数与方差自动标定手填gamma1会让核衰减过快、模型退化为局部插值GPR.predict的返回值恒定返回一个数组只有显式写return_stdTrue时才返回二元组(均值, 标准差)三、常见报错与排查1.ValueError: array must not contain infs or NaNs或拟合极慢现象GPR 拟合时直接抛含 inf/NaN 的错或跑几分钟不结束。根因核矩阵退化样本间高度共线、未做归一化导致 Cholesky 分解失败。解法确认normalize_yTrue、alpha不是一个离谱的小值1e-6 可用0 不行并把输入光谱先做 SNV共线严重的波段可先用第 09 篇的 VIP 或 iPLS 精简。2. GPR 的结果惨不忍睹R²test 掉到 0.15 以下现象明明照抄了 SPEC 的核测试集 R² 却只有 0.1 量级。根因为了复现 SPEC 参数把optimizerNone传了进去核被冻结在length_scale50.0的初值上模型被过度平滑成一条几乎不随波长变化的平线。解法不要设optimizer保持默认SPEC 的length_scale50.0是初值优化后的值见 2.2 的实测输出。3.ValueError: Found input variables with inconsistent numbers of samples现象算覆盖率时报样本数不一致。根因GPR 的predict返回(n,)而PLSRegression.predict返回(n, 1)混用时形状没对齐。解法统一.ravel()。4.TypeError: got an unexpected keyword argument squared现象流程里用mean_squared_error(..., squaredFalse)算 RMSE。根因squared参数在 scikit-learn 1.4 弃用、1.6 移除铁律 1。解法改用root_mean_squared_error。5. SVR 预测值几乎是一条水平线现象SVR 的 R²test 接近 0预测值方差极小。根因epsilon相对 y 的量纲设得太大——把宽容带撑到吞掉全部残差。assay 的通道噪声在 mg 量级epsilon1.0合适设成 20 会直接失效。解法先统计重复测量的 RMSE用它作为 ε 的起点再在验证集上微调。四、动手练习练习 1复现 SVR 锚点跑 2.1 的脚本输出 cal/test 双集指标。判定标准R²test 与 0.9144 的偏差小于 1e-3RMSEP 与 4.6074 的偏差小于 1e-3且svr.n_support_.sum()严格小于 155本系列实测 129。练习 2ε 的稀疏度旋钮把epsilon依次设为 1.0、5.0、10.0、20.0 重跑记录支持向量数与 R²test。判定标准支持向量数随 ε 单调下降ε20 时模型会明显退化预测值趋近于水平线由此验证ε 不能超过 y 的噪声量级。练习 3覆盖率核对 核超参数的代价先把 2.2 的1.96依次换成1.0、1.645、2.576记录三档覆盖率再分别用optimizerNone与默认优化器各拟合一次 GPR。判定标准覆盖率随倍数单调上升1.96一档落在 0.95~0.99 之间本系列实测 0.9761、2.576一档接近 1.0默认优化器下 R²test 不低于 0.90而optimizerNone下 R²test 会塌到 0.15 以下——由此说明核超参数必须被优化。五、小结与下一篇预告SVR 用 ε-不敏感损失把管道内的样本完全忽略用C惩罚越界样本从而在 155 个小样本上取得 R²test 0.9144、RMSEP 4.6074、RPD 3.42——唯一超过 PLS 的模型GPR 则把函数本身当随机变量用ConstantKernel * RBF WhiteKernel三项分别表达幅度、平滑形状、观测噪声normalize_yTrue让超参数初值与 CQA 量纲解耦最终取得 R²test 0.9064并附带一个预测标准差均值 5.0287 mg 的不确定度。核对覆盖率是使用 GPR 的必做动作本系列实测 95% 区间的实测覆盖率为 0.9761略高于名义 0.95属保守但安全。模型能算准了下一步必须回答它凭什么这么算。第 14 篇《模型可解释性PLS 回归系数、载荷与排列重要性》对应实施计划第 14 篇模型可解释性会以 PLS 的回归系数、x_loadings_与 VIP 为主体叠加 scikit-learn 的permutation_importance讲清为什么线性模型不需要 SHAP 这类通用归因工具并把归因结果与 1208–1236 nm 的 C-H 二级倍频区互相印证。本篇认知问题回显FAQQ1SVR 的 C 和 epsilon 分别控制什么Aepsilon是 ε-不敏感管道半宽、量纲与 CQA 相同管道内样本损失为零C是越界惩罚系数越大越贴合训练数据、越易过拟合。Q2SVR 的核技巧到底解决了什么问题A它用 rbf 等核函数直接计算高维映射后的内积exp(−γ‖x_i−x_j‖²)从而不显式构造高维特征就获得非线性拟合能力。Q3GPR 的核写成 ConstantKernel×RBFWhiteKernel三项各管什么AConstantKernel管信号幅度RBF管随距离衰减的平滑形状length_scaleWhiteKernel管观测噪声少了噪声项模型会强行穿过训练点。Q4GPR 里 normalize_yTrue 有什么作用A它把 y 先归一到零均值单位方差使alpha1e-6与核超参数初值不再依赖 CQA 量纲weight 与 assay 才能共用一套配置。Q5GPR 报出的预测标准差可以直接当可信度用吗A不能直接用必须核对覆盖率。本系列实测 95% 区间覆盖率为 0.9761高于名义 0.95 说明区间略宽、偏保守。