PyMC 高斯过程实现类完全指南：Latent、Marginal、HSGP、Kron 与 TP 的选型与使用

PyMC 高斯过程实现类完全指南Latent、Marginal、HSGP、Kron 与 TP 的选型与使用【免费下载链接】pymcBayesian Modeling and Probabilistic Programming in Python项目地址: https://gitcode.com/GitHub_Trending/py/pymc本指南围绕 PyMC 的pymc.gp模块展开系统讲解其八种 GP 实现类——HSGP、HSGPPeriodic、Latent、LatentKron、Marginal、MarginalKron、MarginalApprox与TP——的数学定位、适用场景、核心方法与参数细节。读完你将能够根据数据分布假设与计算规模正确选择 GP 实现掌握prior/marginal_likelihood/conditional/predict的完整调用链并利用prior_linearized以线性模型的方式加速大规模推理。一、认识 PyMC 的 GP 实现家族在 PyMC 中pymc.gp子模块提供了从精确 GP 到各种近似方案的一整套实现。pymc/gp/__init__.py中集中导出了全部实现类它们可以按三条主线划分按是否显式包含噪声建模Latent与HSGP不假设高斯噪声可搭配任意似然如泊松、伯努利Marginal、MarginalKron、MarginalApprox则将 GP 先验与加性高斯噪声合并为边缘似然只适用于正态回归。按精确度与计算代价Latent/Marginal是精确实现HSGP/HSGPPeriodic是固定基向量Hilbert 空间近似MarginalApprox是基于诱导点的稀疏近似LatentKron/MarginalKron利用 Kronecker 结构加速网格数据。按分布族TP是以 Student-t 过程替换高斯过程的变体需要显式指定自由度nu。下表概括了八种实现的核心定位依据pymc/gp/gp.py与pymc/gp/hsgp_approx.py中的类文档实现类数学形式适用场景关键方法Latent隐函数 GP 先验无噪声假设非正态似然分类、计数回归、任意模型组件prior、conditionalMarginalGP 先验 加性高斯噪声正态回归、观测数据平滑marginal_likelihood、conditional、predictMarginalApprox诱导点稀疏近似DTC/FITC/VFE大数据量正态回归marginal_likelihood、conditional、predictLatentKronKronecker 乘积核的隐 GP多维规则网格、无噪声似然prior、conditionalMarginalKronKronecker 乘积核 高斯噪声多维规则网格回归marginal_likelihood、conditional、predictHSGPHilbert 空间基函数近似大样本平滑回归/分类、线性化加速prior、prior_linearized、conditionalHSGPPeriodic周期核的基函数近似1 维周期性数据季节、昼夜prior、prior_linearized、conditionalTPStudent-t 过程需指定nu对异常值更鲁棒的回归prior、conditional所有实现类都继承自pymc/gp/gp.py中的Base基类该基类定义了统一的接口prior(name, X, ...)、marginal_likelihood(name, X, ...)、conditional(name, Xnew, ...)与predict(Xnew, ...)并实现了__add__操作符——两个同类 GP 可以通过相加此时均值函数与协方差函数分别相加见gp.py第 49-55 行。二、Latent无噪声假设的隐函数 GPLatent是 GP 的直接实现不包含任何加性噪声假设因而得名“Latent”——底层函数值本身被视为潜变量。其数学形式为f(x) ~ GP(μ(x), k(x, x))mean_func均值函数默认pm.gp.mean.Zero()零均值。cov_func协方差函数默认pm.gp.cov.Constant(0.0)。典型用法节选自gp.py中Latent.prior的示例import numpy as np import pymc as pm # 一维列向量输入 X np.linspace(0, 1, 10)[:, None] with pm.Model() as model: cov_func pm.gp.cov.ExpQuad(1, ls0.1) gp pm.gp.Latent(cov_funccov_func) f gp.prior(f, XX) # 采样后在新点处构建条件分布 Xnew np.linspace(-1, 2, 50)[:, None] with model: fcond gp.conditional(fcond, XnewXnew)prior的参数细节Latent.prior的完整签名为gp.py第 159 行prior(name, X, n_outputs1, reparameterizeTrue, jitter1e-6, **kwargs)X函数输入值。一维输入必须是形状为(n, 1)的列向量。n_outputs输出 GP 的数量默认 1。例如gp.prior(f, XX, n_outputs3, dims(n_gps, x_dim))要求len(n_gps) 3且len(x_dim) X.shape[0]。reparameterize默认True通过协方差矩阵的 Cholesky 因子旋转随机变量实现重参数化f mu cholesky(cov) vv为标准正态见_build_prior。设为False时直接使用pm.MvNormal。jitter默认1e-6这是模块级常量JITTER_DEFAULT定义于pymc/gp/util.py第 27 行添加到协方差矩阵对角线以保障数值稳定性。conditional的条件分布Latent.conditional(name, Xnew, givenNone, jitter1e-6, **kwargs)基于训练点上的f值构造新点f*的条件高斯分布gp.py第 231 行。其底层_build_conditionalgp.py第 216-229 行执行标准的条件化推导Kxx cov_total(X); Kxs self.cov_func(X, Xnew); Kss self.cov_func(Xnew) L cholesky(stabilize(Kxx, jitter)) A solve_lower(L, Kxs) mu mean_func(Xnew) A^T (L^{-1} (f - mean_func(X))) cov Kss - A^T A即条件均值通过核矩阵求解获得条件协方差为舒尔补Schur complement。given参数允许传入{X: ..., f: ..., gp: ...}覆盖已存储的训练信息这在把Latent用作较大模型组件时非常有用。三、MarginalGP 先验与高斯噪声的合并Marginal实现的是“GP 先验 加性高斯噪声”的联合模型适用于正态回归。其边缘似然是对 GP 先验与正态似然乘积的积分y | X, θ ~ ∫ p(y | f, X, θ) p(f | X, θ) df典型用法节选自gp.py第 425-446 行X np.linspace(0, 1, 10)[:, None] with pm.Model() as model: cov_func pm.gp.cov.ExpQuad(1, ls0.1) gp pm.gp.Marginal(cov_funccov_func) sigma pm.HalfCauchy(sigma, beta3) y_ gp.marginal_likelihood(y, XX, yy, sigmasigma) Xnew np.linspace(-1, 2, 50)[:, None] with model: fcond gp.conditional(fcond, XnewXnew)marginal_likelihood参数说明签名见gp.py第 456 行marginal_likelihood(name, X, y, sigma, jitter1e-6, is_observedTrue, **kwargs)y观测数据形状(n,)是 GP 函数与高斯噪声之和。sigma高斯噪声标准差。可以是标量/随机变量也可以是协方差对象当传入非BaseCovariance的标量时内部会包装为pm.gp.cov.WhiteNoise(sigma)见gp.py第 497 行从而支持异方差或结构化噪声。is_observed默认True将y设为模型中的观测变量该参数已被标记为弃用。底层_build_marginal_likelihoodgp.py第 449-454 行只是简单地把核协方差与噪声协方差相加cov Kxx Knx再经stabilize(cov, jitter)加固。边际化带来的好处是采样时无需为每个观测点引入潜变量fHMC 的采样空间大幅缩小。predict与conditionalpredict(Xnew, pointNone, diagFalse, pred_noiseFalse, givenNone, jitter1e-6, modelNone)返回预测均值和可选方差的 NumPy 数组diagTrue时只返回对角方差计算更省pred_noiseTrue时预测中包含噪声方差适合与观测比较。conditional(name, Xnew, pred_noiseFalse, givenNone, jitter1e-6, **kwargs)构造预测点上的随机变量。pred_noiseTrue时条件分布包含噪声项用于模拟新的观测y否则仅给出隐函数f的分布。测试tests/gp/test_gp.py中的testMarginal系列用例展示了marginal_likelihood、conditional与predict的组合验证方式。四、MarginalApprox诱导点稀疏近似DTC / FITC / VFE当训练样本量很大时精确Marginal需要对n × n核矩阵做 Cholesky 分解O(n³)代价过高。MarginalApprox通过一小撮诱导点inducing pointsXu对全数据做低秩近似将复杂度降至O(n m²)m为诱导点个数。可用近似方式gp.py第 677-681 行DTCDeterministic Training Conditional确定性训练条件近似忽略条件方差修正。FITCFully Independent Training Conditional完全独立训练条件近似对每个点单独修正方差对噪声建模更细。VFEVariational Free Energy变分自由能近似来自 Titsias 的变分诱导变量方法是默认选项approxVFE。示例节选自gp.py第 694-720 行X np.linspace(0, 1, 10)[:, None] Xu np.linspace(0, 1, 5)[:, None] # 远小于 X 的诱导点集 with pm.Model() as model: cov_func pm.gp.cov.ExpQuad(1, ls0.1) gp pm.gp.MarginalApprox(cov_funccov_func, approxFITC) sigma pm.HalfCauchy(sigma, beta3) y_ gp.marginal_likelihood(y, XX, XuXu, yy, sigmasigma) Xnew np.linspace(-1, 2, 50)[:, None] with model: fcond gp.conditional(fcond, XnewXnew)三种近似的底层差异MarginalApprox的marginal_likelihood需要额外传入诱导点Xu。其底层_build_marginal_likelihood_loglikgp.py第 750-777 行用代码直观地区分了三种近似对 FITC每个点的噪声项被修正为Lamd clip(Kffd - Qffd, 0, inf) sigma²即保留训练点自身的方差细节对 VFELamd sigma²且额外计算了一个迹trace修正项(1/(2σ²)) * (Σ Kffd - Σ Qffd)这是变分下界推导的自然产物对 DTC同样Lamd sigma²但不含迹修正实现最简单。参数校验方面构造函数只接受(FITC, VFE, DTC)之一否则抛出NotImplementedErrorgp.py第 736-740 行__add__还要求相加的两个 GP 必须使用相同的approxgp.py第 742-748 行。诱导点的选取可以参考pymc/gp/util.py中提供的kmeans_inducing_points(n_inducing, X, **kmeans_kwargs)工具函数用 K-means 自动生成诱导点位置。tests/gp/test_gp.py中的TestMarginalApprox会针对 DTC/FITC/VFE 分别断言边缘似然、条件分布与predict的结果一致性。五、HSGPHilbert 空间基函数近似HSGPHilbert Space Gaussian Process是一种降秩 GP 近似它使用一组固定的基向量拉普拉斯特征函数其系数是平稳核功率谱密度的随机函数hsgp_approx.py第 171-183 行。与Latent一样它不假设高斯噪声可搭配任意似然与Latent的不同在于它只在有限边界域[-L, L]内逼近核函数因此适用范围严格限于训练与预测都在该边界内的场景。构造参数HSGP(m[25, 25], c4.0, cov_funccov_func) # 或 L[...]完整签名为hsgp_approx.py第 258-268 行HSGP(m, LNone, cNone, drop_firstFalse, parametrizationnoncentered, *, mean_funcZero(), cov_funcCovariance)m每个活跃维度的基向量个数列表。len(m)必须等于协方差函数的活跃维度数cov_func.n_dims否则抛出ValueError。总基向量数m* prod(m)保存在gp.n_basis_vectors属性中。示例中m[25, 25]意味着共25 × 25 625个基向量。L/c边界条件。二者必须且只能提供一个hsgp_approx.py第 281-282 行。若提供c则L c * S其中S是数据半宽(max(X) - min(X))/2官方建议c 1.2低于此值会发出警告hsgp_approx.py第 287-288 行。边界L的设定与set_boundary、calc_eigenvalues、calc_eigenvectors三个模块级函数直接相关hsgp_approx.py第 33-74 行特征值由(π·S_arr/(2L))²解析给出特征向量则是正弦基函数。drop_first默认False。第一个基向量往往“很平”、与截距项高度相似当模型已有截距时丢弃它可能改善采样。该参数已被标记为未来版本将弃用DeprecationWarning。parametrizationnoncentered默认或centered。非中心化时系数beta ~ Normal(0, 1)再乘以功率谱密度平方根中心化时beta ~ Normal(0, sqrt_psd)。当 GP 信号强于噪声时中心化参数化可能采样效率更高。cov_func必须是实现power_spectral_density方法的平稳核如ExpQuad、Matern52。注意Periodic核不适用它有专门的HSGPPeriodic。基向量个数与边界的启发式选择pymc/gp/hsgp_approx.py提供了approx_hsgp_hyperparams(x_range, lengthscale_range, cov_func)工具函数hsgp_approx.py第 97-168 行基于 Ruitort-Mayol 等人的建议自动给出最小m与cx_range应覆盖训练与预测的全部范围。例如训练数据在[0, 10]而要在[7, 15]预测就应传x_range[0, 15]。lengthscale_range依据你对长度尺度的先验来定。例如认为长度尺度 95% 先验质量落在[1, 5]就传lengthscale_range[1, 5]。支持的cov_func取值为expquad、matern52、matern32大小写不敏感其余值抛出ValueError。返回的(m, c)中c越大越能容纳大长度尺度m越大越能刻画小长度尺度但计算代价随之上升。注意这些建议基于一维 GP。prior / prior_linearized / conditionalHSGP 的核心是固定基结构因此它额外提供了prior_linearized(X)方法hsgp_approx.py第 329-425 行返回(phi, sqrt_psd)——即拉普拉斯特征函数矩阵形状(n, m*)与功率谱密度平方根向量。这允许绕过 GP 接口、按线性回归的方式建模with pm.Model() as model: eta pm.Exponential(eta, lam1.0) ell pm.InverseGamma(ell, mu5.0, sigma5.0) cov_func eta**2 * pm.gp.cov.ExpQuad(1, lsell) gp pm.gp.HSGP(m[200], L[10], cov_funccov_func) X pm.Data(X, X) phi, sqrt_psd gp.prior_linearized(XX) beta pm.Normal(beta, sizegp.n_basis_vectors) f pm.Deterministic(f, phi (beta * sqrt_psd))线性化形式的优势在于预测时只需pm.set_data({X: x_new})更新数据再用pm.sample_posterior_predictive(idata, var_names[f])生成后验预测整个流程与线性模型无异多个 GP 可以共享同一组基phi只各自维护自己的系数实现计算加速。实现细节上prior_linearized会固定训练数据的中点_X_center预测时用训练中点而不是测试中点做中心化确保基的一致性hsgp_approx.py第 403-407 行。prior(name, X, dimsNone, hsgp_coeffs_dimsNone)与conditional(name, Xnew, dimsNone)则是对线性化形式的封装前者按parametrization生成系数{name}_hsgp_coeffs并构造f mean_func(X) phi (beta * sqrt_psd)hsgp_approx.py第 456-473 行后者基于已存系数构造新点的确定性输出hsgp_approx.py第 476-514 行。在测试层面tests/gp/test_hsgp_approx.py的TestHSGP会校验m/L/c/parametrization的各种非法组合均抛出ValueError第 149-171 行用c2与L[12]两种方式等价构造 HSGP第 175-179 行通过假设检验KS 检验验证 HSGP 先验与未近似的pm.gp.Latent先验在分布上不可区分第 205-224 行条件分布同理第 234-251 行。六、HSGPPeriodic周期核的基函数近似HSGPPeriodic是针对Periodic协方差函数的近似。严格来说它并非 Hilbert 空间近似而是同一篇论文Ruitort-Mayol et al., 2022 附录 B中基于随机谐振子的级数展开但 API 与HSGP保持一致可作为Latent的即插即用替代hsgp_approx.py第 517-530 行。with pm.Model() as model: scale pm.HalfNormal(scale, 10) cov_func pm.gp.cov.Periodic(1, period1, ls0.1) gp pm.gp.HSGPPeriodic(m25, scalescale, cov_funccov_func) f gp.prior(f, XX) Xnew np.linspace(-1, 2, 50)[:, None] with model: fcond gp.conditional(fcond, XnewXnew)参数说明hsgp_approx.py第 535-543 行m基向量个数必须为正整数。该近似仅实现于一维情形cov_func必须是Periodic实例且n_dims 1否则抛出ValueErrorhsgp_approx.py第 586-602 行。scaleGP 效应的标准差方差平方根默认1.0可用随机变量。cov_funcPeriodic核。注意方差控制通过scale完成不要在核里重复放幅度参数。底层的calc_basis_periodichsgp_approx.py第 77-94 行用角频率w0 2π/period构造余弦基cos(j·w0·X)与正弦基sin(j·w0·X)j 0..m-1。由于正弦分量的第一个特征函数恒为零实际使用的系数个数为2m - 1见prior中size(m*2-1)的写法hsgp_approx.py第 717-730 行。线性化接口prior_linearized返回((phi_cos, phi_sin), psd)同样支持pm.set_data式的线性模型预测。tests/gp/test_hsgp_approx.py的TestHSGPPeriodic校验了m非法值、非Periodic核、多维周期核等错误路径第 261-281 行并同样用假设检验确认其先验与Latent先验一致第 293-311 行。七、LatentKron 与 MarginalKronKronecker 结构加速当输入覆盖多维规则网格时完整核矩阵大小N × NN prod(n_d)会迅速失控。若核函数本身可分解为各维度核的 Kronecker 乘积则可以利用KroneckerNormal分布与kron_dot等工具避免显式构造巨大协方差矩阵。LatentKronLatentKron是不含噪声假设的 Kronecker 结构 GPgp.py第 910-923 行构造时传入各维度独立的协方差函数列表X1 np.linspace(0, 1, 10)[:, None] X2 np.linspace(0, 2, 5)[:, None] Xs [X1, X2] with pm.Model() as model: cov_func1 pm.gp.cov.ExpQuad(1, ls0.1) cov_func2 pm.gp.cov.ExpQuad(1, ls0.3) gp pm.gp.LatentKron(cov_funcs[cov_func1, cov_func2]) f gp.prior(f, XsXs) Xnew1 np.linspace(-1, 2, 10)[:, None] Xnew2 np.linspace(0, 3, 10)[:, None] Xnew np.concatenate((Xnew1, Xnew2), axis1) # 非完整网格也可以 # Xnew pm.math.cartesian(Xnew1, Xnew2) # 完整网格同样支持 with model: fcond gp.conditional(fcond, XnewXnew)cov_funcs协方差函数列表内部会被包装为pm.gp.cov.Krongp.py第 962-968 行。prior的Xs是各维度输入列表其总协方差定义在完整网格cartesian(*Xs)上每个X[i]必须能无报错地传入对应的cov_funcs[i]。内部实现上_build_prior对各维度核分别做 Cholesky 分解再用kron_dot将旋转变量映射到全网格gp.py第 974-982 行避免构造N × N矩阵。LatentKron不支持加法__add__抛TypeError。MarginalKronMarginalKron是 Kronecker 结构 高斯噪声的回归版本接口融合了Marginal与Kronkron_gp pm.gp.MarginalKron(mean_funcself.mean, cov_funcsself.cov_funcs) f kron_gp.marginal_likelihood(f, self.Xs, self.y, sigmaself.sigma)其marginal_likelihood需要Xs各维度输入列表、y按完整网格展平的观测形状(N,)与sigmapredict支持diagTrue只返回对角方差。测试tests/gp/test_gp.py中的TestMarginalKron第 467 行起专门对比了MarginalKron与普通Marginal在相同数据上的边缘似然和预测结果验证 Kronecker 实现与精确实现在数值上等价。八、TPStudent-t 过程TP是Latent的学生 t 版本对异常值更鲁棒代价是必须指定自由度参数nu且不支持加法gp.py第 272-301 行f(X) ~ TP(μ(X), k(X, X), ν)with pm.Model() as model: cov_func pm.gp.cov.ExpQuad(1, ls0.1) gp pm.gp.TP(cov_funccov_func, nu3) # 或 scale_funccov_func, nunu f gp.prior(f, XX)nu自由度必填未提供时抛出ValueErrorgp.py第 303-305 行。scale_func协方差尺度函数旧参数名cov_func仍可用但会触发FutureWarning提示改用scale_funcgp.py第 306-312 行。底层实现中_build_prior使用pm.StudentT重参数化时或pm.MvStudentT非重参数化时构造旋转变量_build_conditionalgp.py第 360-372 行则把自由度更新为nu2 nu n条件协方差按(nu beta - 2)/(nu2 - 2)缩放最终以pm.MvStudentT返回。九、选型决策指南综合以上分析选择实现类时可以按三步决策确定似然/噪声假设需要高斯噪声边缘化 →Marginal/MarginalKron/MarginalApprox非正态似然分类、计数、异质噪声或 GP 作为复杂模型组件 →Latent/HSGP/HSGPPeriodic。评估数据规模与结构样本量大 → 优先HSGP固定基m*个基向量或MarginalApprox诱导点m个诱导点多维规则网格 →LatentKron/MarginalKron周期数据 →HSGPPeriodic。需要鲁棒性数据含离群点 →TP指定nu。对于中等规模数据的精确回归Marginal是最直接的起点当模型复杂化或数据增长后HSGP的prior_linearized将 GP 降维为“基 × 系数”的线性形式配合pm.set_data与pm.sample_posterior_predictive即可在保持贝叶斯推断框架的同时获得接近线性模型的扩展性。十、进一步阅读GP 模块的 API 汇总docs/source/api/gp.rst及其下的docs/source/api/gp/implementations.rst、docs/source/api/gp/cov.rst、docs/source/api/gp/mean.rst、docs/source/api/gp/util.rst。全部实现源码pymc/gp/gp.pyTP、Latent、Marginal、MarginalApprox、LatentKron、MarginalKron与pymc/gp/hsgp_approx.pyHSGP、HSGPPeriodic及基函数计算工具。协方差函数库含power_spectral_density实现pymc/gp/cov.py。工具函数JITTER_DEFAULT、stabilize、kmeans_inducing_points、plot_gp_dist等pymc/gp/util.py。配套测试tests/gp/test_gp.py、tests/gp/test_hsgp_approx.py提供了各实现与精确 GP 一致性验证、参数校验与预测逻辑的完整示例。【免费下载链接】pymcBayesian Modeling and Probabilistic Programming in Python项目地址: https://gitcode.com/GitHub_Trending/py/pymc创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考