机器学习中的高斯过程：MQL5中的回归模型·进阶篇
（2/3）·当模型不再只吐一个点位，而是给出概率区间，你才真正看清预测的可信度边界
不少交易者把回归模型当成确定性信号源，看到数值就直接下单，却忽略预测本身可能误差极大。高斯过程的价值恰在于它顺手把不确定性也算了出来，而这层信息常被当作噪声丢掉。接上篇对核函数与先验的铺垫，本篇继续深挖贝叶斯推导下的完整回归链路。
似然把先验和报价连起来
高斯过程里，似然 p(y|f, X) 说的是：已知先验里的隐函数 f 和输入集 X，看到真实标签 y 的概率有多大。它是贝叶斯推断的枢纽，把函数先验和盘面真实观测接上。完整数据集下，似然走的是多元正态分布。 公式里 f(x) 是隐函数在观测点 X 上的取值向量，σy² 是噪声方差，‖y−f(x)‖² 是观测和隐函数之差的欧氏范数平方。这意味着 f(x) 越贴着 y，似然就越高——放在 MT5 做回归预测时，就是模型曲线越贴合已走出的 K 线，拟合可信度越大。 外汇和贵金属波动带杠杆，用这类概率模型只代表后市倾向，不预示必然，实盘前务必用历史数据回测验证。
「贝叶斯更新后得到的可信函数集」
后验分布是给定观测数据后，用似然对先验做贝叶斯更新得到的结果，它描述隐函数 f(x) 在新测试点 X* 上的取值分布。 均值 μ* 是更新后的预测值，协方差矩阵 Σ* 衡量预测不确定性——这两个量直接决定你敢不敢跟这一笔信号。 先验给出候选函数集，似然按数据给权重，后验给出最可信的函数及其误差范围。做外汇或贵金属回测时，Σ* 放大往往意味着行情结构突变，此时盲信 μ* 概率上要吃亏，属高风险场景。
◍ 高斯似然下怎么解回归
把观测噪声设成高斯分布后，高斯过程回归可以走纯线性代数解析解，不用上 MCMC 或复杂数值优化，计算量直接降一档。它既能吃带噪数据做回归，也能在无噪场景做严格插值，两条路的公式要分清楚。 无噪时就是插值：训练集 D={(x_n,y_n)}，且 y_n=f(x_n) 精确无误差。训练点 X 与测试点 X* 的函数值联合先验是多元正态，观测等于隐函数真值（似然为 δ 函数），于是 f* 的后验直接由条件分布公式给出——后验均值 μ_f* 是最优预测，后验协方差 Σ_f* 量化测试点不确定性，训练点处对角元为 0，采样函数必穿点。 带噪时 y=f(X)+ε，ε∼N(0,σ²)。这时后验要把 K(X,X) 对角加上 σ² 再算，模型不再硬穿每个点，而是在噪声允许带宽内平滑拟合。σ 和核参数一起优化。 一个容易混的点：加噪后 f* 和 y* 的预测分布不一样。μ_y*=μ_f*（噪声期望 0），但 Σ_y*=Σ_f*+σ²I。要新观测值区间就用 y* 那套，只要隐函数区间就用 f* 那套。 图例 5 用 5 个训练样本 + RBF 核画了三条无噪后验采样曲线，全过训练点；图例 6 同设定带噪，灰虚线是 f* 的 2σ 区间（约 95% 概率），远离数据处明显变宽。注意那灰带是 f* 的区间，不是 y* 的。外汇与贵金属行情噪声结构常非平稳，直接套固定 σ² 有高估近区、低估远区风险。
乔列斯基分解省下的算力不是小数
高斯过程拟合得好不好，核函数参数和噪声方差这组超参数占了大头。直接算 (K+σ²I)⁻¹·y 在样本量 n 变大后基本跑不动——n×n 矩阵求逆的复杂度是 O(n³)，n 上到几千就卡死。 工程上绕开求逆的办法是乔列斯基分解：把 K+σ²I 拆成下三角矩阵 L 和其转置 Lᵀ 的乘积，再解 L·z=y 得到 z=L⁻¹y。负对数边缘似然（NLML）三项随之改写：数据项变成 ½·zᵀz，行列式项变成对角元对数之和 ∑log(L_ii)，常数项仍是 n/2·log(2π)。 这套计算在 NegativeLogMarginalLikelihood 里封装好了，协方差矩阵 K 由 ComputeKernelMatrix 对所有选定核（RBF、线性、周期等）求和得到。优化目标 NLML 交给 OptimizeGP 调用 ALGLIB 的 BLEIC 算法去最小化。 别把正态当圣经：代码里 jitter=1e-6 不是装饰，K 接近奇异时没这行直接 Cholesky 报错返回 DBL_MAX，优化器会废。开 MT5 把 jitter 调到 1e-8 试一次，看分解失败率是否上升。
class=class="str">"cmt">//+------------------------------------------------------------------+ class=class="str">"cmt">//| Calculate the negative log-likelihood | class=class="str">"cmt">//+------------------------------------------------------------------+ class="type">class="kw">double NegativeLogMarginalLikelihood(class="kw">const matrix &x, class="kw">const matrix &y, class="type">int &kernel_list[], KernelParams ¶ms) { class="type">int n = (class="type">int)x.Rows(); class=class="str">"cmt">//--- Calculate the K covariance matrix matrix K = ComputeKernelMatrix(x, x, kernel_list, params); class=class="str">"cmt">//--- Add white noise variance if KERNEL_WHITE is selected class="type">class="kw">double white_noise_variance = class="num">0.0; for (class="type">int i = class="num">0; i < ArraySize(kernel_list); i++) { if (kernel_list[i] == KERNEL_WHITE) { white_noise_variance = params.white.sigma * params.white.sigma; break; } } class=class="str">"cmt">// Adding a little jitter prevents problems, class=class="str">"cmt">// if the K matrix is singular or close to singular class=class="str">"cmt">// due to the peculiarities of the calculations or the chosen kernel class="type">class="kw">double jitter = class="num">1e-6; K += matrix::Identity(n, n) * (white_noise_variance + jitter); class=class="str">"cmt">// ---Perform the Cholesky decomposition: K + sigma^class="num">2*I = L * L^T matrix L; if (!K.Cholesky(L)) { Print("Error: Cholesky decomposition failed"); class="kw">return DBL_MAX; } class=class="str">"cmt">//--- Solve the linear system L * z = y to find z = L^(-class="num">1) * y vector z = L.LstSq(y.Col(class="num">0)); if (z.Size() == class="num">0) { Print("Error: Unable to solve the system L * z = y"); class="kw">return DBL_MAX; } class=class="str">"cmt">//--- Calculate the first term in the NLML equation: class=class="str">"cmt">// class="num">1/class="num">2 * y^T * (K + sigma^class="num">2*I)^-class="num">1 * y = class="num">1/class="num">2 * z^T * z class="type">class="kw">double data_term = class="num">0.5 * z @ z; class=class="str">"cmt">// Scalar product z^T * z class=class="str">"cmt">//--- Calculate the second term: class="num">1/class="num">2 * log|K + sigma^class="num">2*I| = sum(log(L_ii)) vector diag = L.Diag(); class="type">class="kw">double log_det = class="num">0.0; for (class="type">int i = class="num">0; i < n; i++) log_det += MathLog(diag[i]); class=class="str">"cmt">// sum the logarithms of the L diagonal elements class=class="str">"cmt">//--- Calculate the third term: n/class="num">2 * log(2π) class="type">class="kw">double const_term = class="num">0.5 * n * MathLog(class="num">2 * M_PI); class="kw">return data_term + log_det + const_term; } class=class="str">"cmt">//+------------------------------------------------------------------+ class=class="str">"cmt">//| Calculate the covariance matrix | class=class="str">"cmt">//+------------------------------------------------------------------+ matrix ComputeKernelMatrix(class="kw">const matrix &X1, class="kw">const matrix &X2, class="type">int &kernel_list[], KernelParams ¶ms) {
「高斯过程超参的自动寻优落地」
把多个核拼进协方差矩阵后，下一步是让模型自己挑超参。下面这段直接调 Alglib 的 BLEIC 求解器，对 RBF、线性、周期等核的参数做有界非线性最小化，目标函数封装在 CNDimensional_GP 里（负对数边际似然 NLML）。 初始化阶段由 InitializeKernelParams 按 kernel_list 里的核类型算出参数总数，并填好初值 w、归一化尺度 s、上下界 bndl/bndu。BLEIC 的停止条件里 epsf 设了 0.0001、epso 与 epsi 都是 0.00001，diffstep 取 0.0001 做数值导数——这几个数直接决定收敛速度和过拟合倾向，外汇与贵金属行情噪声大，参数边界最好收窄再跑。 跑完用 MinBLEICResults 把最优 w 取回，Print 出 TerminationType；若返回 -8 说明内部完整性检查撞到无穷或 NaN，多半是初值尺度 s 没设对。开 MT5 把这段塞进 EA，接上你自己的 x_train/y_train，就能让小布式管线替你扫一遍核组合。
matrix K = matrix::Zeros(X1.Rows(), X2.Rows()); for (class="type">int i = class="num">0; i < ArraySize(kernel_list); i++) { class="kw">switch (kernel_list[i]) { case KERNEL_RBF: K += RBF_kernel(X1, X2, params.rbf.sigma_f, params.rbf.length); break; case KERNEL_LINEAR: K += Linear_kernel(X1, X2, params.linear.sigma_l); break; case KERNEL_PERIODIC: K += Periodic_kernel(X1, X2, params.periodic.sigma_f, params.periodic.length, params.periodic.period); break; case KERNEL_WHITE: class=class="str">"cmt">// WhiteKernel is added separately as sigma^class="num">2 * I break; } } class="kw">return K; } class=class="str">"cmt">//+------------------------------------------------------------------+ class=class="str">"cmt">//| Gaussian process hyperparameter optimization | class=class="str">"cmt">//+------------------------------------------------------------------+ KernelParams OptimizeGP(class="type">int &kernel_list[], matrix &x_train, matrix &y_train) { class="type">class="kw">double w[]; class=class="str">"cmt">// Array of initial hyperparameter values class="type">class="kw">double s[]; class=class="str">"cmt">// Array of scales for hyperparameters used for normalization in optimization class="type">class="kw">double bndl[], bndu[]; class=class="str">"cmt">// Arrays of lower and upper bounds for each hyperparameter CObject Obj; CNDimensional_GP ffunc(kernel_list, x_train, y_train);class=class="str">"cmt">// Object of a class implementing the NLML target function CNDimensional_Rep frep; class=class="str">"cmt">// Object for storing the optimization report class=class="str">"cmt">// Initialize parameters InitializeKernelParams(kernel_list, w, s, bndl, bndu); class=class="str">"cmt">/* The total number of parameters(num_params) is calculated depending on the kernel types in kernel_list. The initial values of the w parameters, scale and bounds for the parameters are set */ CMinBLEICStateShell state; CMinBLEICReportShell rep; class=class="str">"cmt">// Object for storing optimization results. class="type">class="kw">double epsg = class="num">0; class=class="str">"cmt">// Gradient precision(class="num">0 means gradient stopping is disabled) class="type">class="kw">double epsf = class="num">0.0001; class=class="str">"cmt">// Precision by function value class="type">class="kw">double epsw = class="num">0; class=class="str">"cmt">// accuracy by parameters class="type">class="kw">double diffstep = class="num">0.0001; class=class="str">"cmt">// Step for numerical calculation of derivatives. class="type">class="kw">double epso = class="num">0.00001; class=class="str">"cmt">// Parameters for external and internal convergence conditions in BLEIC class="type">class="kw">double epsi = class="num">0.00001; CAlglib::MinBLEICCreateF(w, diffstep, state); class=class="str">"cmt">// Create a BLEIC optimization object with w initial parameters and diffstep for numerical derivatives CAlglib::MinBLEICSetBC(state, bndl, bndu); class=class="str">"cmt">// Set parameter bounds(bndl, bndu). CAlglib::MinBLEICSetScale(state, s); class=class="str">"cmt">// Sets the scale of(s) parameters. CAlglib::MinBLEICSetInnerCond(state, epsg, epsf, epsw); CAlglib::MinBLEICSetOuterCond(state, epso, epsi); CAlglib::MinBLEICOptimize(state, ffunc, frep, class="num">0, Obj); CAlglib::MinBLEICResults(state, w, rep); Print("TerminationType =", rep.GetTerminationType()); class=class="str">"cmt">// Optimization completion code /* * -class="num">8 internal integrity control detected infinite or