高斯过程回归与分类：不确定性建模实战

f(x) ~ GP(m(x), k(x, x'))

f = [f(x_1), ..., f(x_n)]ᵀ ~ N(μ, K)
μ_i = m(x_i)
K_ij = k(x_i, x_j)

import numpy as np

def sample_gp_prior(n_points=200, n_samples=5, length_scale=0.3, seed=0):
    """从 GP 先验采样：无需数据，纯由核决定函数形状"""
    rng = np.random.default_rng(seed)
    X = np.linspace(0, 1, n_points)[:, None]
    sq = (X - X.T) ** 2
    K = np.exp(-0.5 * sq / length_scale ** 2)      # RBF 核
    K += 1e-8 * np.eye(n_points)                   # 抖动保证数值稳定
    L = np.linalg.cholesky(K)
    return X.ravel(), (L @ rng.standard_normal((n_points, n_samples)))

X, samples = sample_gp_prior()
print("先验样本形状:", samples.shape)   # (200, 5)，5 条光滑随机函数

def zero_mean(X):
    return np.zeros(len(X))

def linear_mean(X):
    return 2.0 * X.ravel()      # 假设有上升趋势

def rbf(X1, X2, length_scale=1.0, variance=1.0):
    sq = (X1[:, None, :] - X2[None, :, :]) ** 2
    return variance * np.exp(-0.5 * sq.sum(-1) / length_scale ** 2)

def matern32(X1, X2, length_scale=1.0, variance=1.0):
    d = np.sqrt(((X1[:, None, :] - X2[None, :, :]) ** 2).sum(-1))
    r = np.sqrt(3) * d / length_scale
    return variance * (1 + r) * np.exp(-r)

X = np.array([[0.0], [0.5], [1.0]])
print("RBF 核矩阵:\n", rbf(X, X).round(3))
print("Matern 3/2 核矩阵:\n", matern32(X, X).round(3))

def combined_kernel(X1, X2, ls1=1.0, ls2=0.1):
    return rbf(X1, X2, ls1) + 0.5 * rbf(X1, X2, ls2)   # 长短两个尺度

μ* = k*ᵀ (K + σ_n²I)⁻¹ y
σ*² = k(x*, x*) - k*ᵀ (K + σ_n²I)⁻¹ k*

def gp_predict(X_train, y_train, X_test, length_scale=0.3, noise=0.05):
    K = rbf(X_train, X_train, length_scale) + noise ** 2 * np.eye(len(X_train))
    Ks = rbf(X_train, X_test, length_scale)
    Kss = rbf(X_test, X_test, length_scale)
    L = np.linalg.cholesky(K)
    alpha = np.linalg.solve(L.T, np.linalg.solve(L, y_train))
    mu = Ks.T @ alpha
    v = np.linalg.solve(L, Ks)
    cov = Kss - v.T @ v
    return mu, np.sqrt(np.clip(np.diag(cov), 0, None))

rng = np.random.default_rng(0)
X_tr = rng.uniform(0, 1, 15)[:, None]
y_tr = np.sin(2 * np.pi * X_tr.ravel()) + rng.normal(0, 0.1, 15)
X_te = np.linspace(-0.2, 1.2, 100)[:, None]
mu, std = gp_predict(X_tr, y_tr, X_te)
print("测试点预测前 3 个:", mu[:3].round(3))
print("对应标准差:", std[:3].round(3))

log p(y | X, θ) = -½ yᵀ (K_θ + σ_n²I)⁻¹ y
                  - ½ log|K_θ + σ_n²I|
                  - n/2 log(2π)

from scipy.optimize import minimize

def neg_log_marginal_likelihood(log_params, X, y):
    ls, var, noise = np.exp(log_params)          # 在对数空间优化，保证正
    K = var * rbf(X, X, ls) + noise ** 2 * np.eye(len(X))
    L = np.linalg.cholesky(K)
    alpha = np.linalg.solve(L.T, np.linalg.solve(L, y))
    return 0.5 * y @ alpha + np.log(np.diag(L)).sum() + 0.5 * len(X) * np.log(2 * np.pi)

res = minimize(neg_log_marginal_likelihood, np.log([0.3, 1.0, 0.1]),
               args=(X_tr, y_tr), method="L-BFGS-B")
print("最优超参数 (ls, var, noise):", np.exp(res.x).round(3))

1. 用牛顿法找后验众数 f̂
2. 在 f̂ 处算 Hessian，得到近似高斯 N(f̂, H⁻¹)
3. 对新点用该高斯做预测积分

from sklearn.gaussian_process import GaussianProcessClassifier
from sklearn.gaussian_process.kernels import RBF
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split

X, y = make_classification(n_samples=200, n_features=2, n_redundant=0,
                           n_informative=2, random_state=0)
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3, random_state=0)
gpc = GaussianProcessClassifier(kernel=RBF(length_scale=1.0), random_state=0)
gpc.fit(X_tr, y_tr)
print("测试精度:", round(gpc.score(X_te, y_te), 4))

# 概念示意：用诱导点近似核矩阵的低秩结构
import numpy as np

def inducing_approx(X, Z, length_scale=0.3):
    """Z: (m, d) 诱导点，返回近似核矩阵"""
    Kmm = rbf(Z, Z, length_scale)
    Knm = rbf(X, Z, length_scale)
    Kmm_inv = np.linalg.pinv(Kmm)
    return Knm @ Kmm_inv @ Knm.T      # 低秩近似，秩 <= m

Z = np.linspace(0, 1, 8)[:, None]     # 8 个诱导点
print("近似核矩阵形状:", inducing_approx(X_tr, Z).shape)

import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel, ConstantKernel, Matern
from sklearn.preprocessing import StandardScaler

rng = np.random.default_rng(42)
X = rng.uniform(0, 5, 60)[:, None]
y = np.sin(X.ravel()) + 0.1 * X.ravel() + rng.normal(0, 0.15, 60)

scaler = StandardScaler().fit(X)
X_s = scaler.transform(X)

# 核 = 常数尺度 × RBF + 白噪声（学习噪声水平）
kernel = ConstantKernel(1.0) * RBF(length_scale=1.0) + WhiteKernel(noise_level=0.1)
gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10,
                              normalize_y=True, random_state=0)
gp.fit(X_s, y)
print("学到的核:", gp.kernel_)

X_test = scaler.transform(np.linspace(-1, 6, 200)[:, None])
mu, std = gp.predict(X_test, return_std=True)
lower, upper = mu - 1.96 * std, mu + 1.96 * std
print("预测区间宽度（端点）:", round((upper - lower)[0], 3))
print("预测区间宽度（中部）:", round((upper - lower)[100], 3))

# 从简单到复杂依次尝试
kernels = {
    "RBF": RBF(length_scale=1.0),
    "Matern32": Matern(length_scale=1.0, nu=1.5),
    "Matern52": Matern(length_scale=1.0, nu=2.5),
    "RBF+noise": RBF(length_scale=1.0) + WhiteKernel(noise_level=0.1),
}
for name, k in kernels.items():
    g = GaussianProcessRegressor(kernel=k, normalize_y=True, random_state=0)
    g.fit(X_s, y)
    print(f"{name}: log-marginal-likelihood = {g.log_marginal_likelihood_value_:.2f}")

1. 用 GP 代理目标函数（给出均值与方差）
2. 用采集函数（EI / UCB / PI）决定下一个评估点
   - 均值高 → 利用（exploit）
   - 方差大 → 探索（explore）
3. 评估该点，更新 GP，重复

import numpy as np
from scipy.stats import norm

def expected_improvement(mu, std, best, xi=0.01):
    """EI 采集函数：平衡利用与探索"""
    z = (mu - best - xi) / (std + 1e-9)
    return (mu - best - xi) * norm.cdf(z) + std * norm.pdf(z)

mu = np.array([1.0, 1.5, 0.5])
std = np.array([0.1, 0.3, 0.8])
print("EI:", expected_improvement(mu, std, best=1.2).round(4))

样本量 < 1 万 且 需要不确定性？
  ├── 是 → GP（RBF/Matern + WhiteKernel）
  └── 否
        ├── 只需预测值 → 梯度提升树 / 神经网络
        ├── 超参搜索 → GP + 贝叶斯优化
        ├── 时序预测 → GP + 周期核（或 ARIMA）
        └── 大数据回归 → 稀疏 GP / 深度核 / 换模型