高斯过程

import numpy as np

def kernel(X1, X2, l=1.0, sigma_f=1.0):
    '''
    各向同性平方指数内核.
    计算点X1与点X2的协方差矩阵.
    
    Args:
        X1: ndArray， m个点 (m x d).
        X2: ndArray， n个点 (n x d).

    返回:
        协方差矩阵 (m x n).
    '''
    sqdist = np.sum(X1**2, 1).reshape(-1, 1) + np.sum(X2**2, 1) - 2 * np.dot(X1, X2.T)
    return sigma_f**2 * np.exp(-0.5 / l**2 * sqdist)

# 编写绘图函数
import matplotlib.pyplot as plt

from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

def plot_gp(mu, cov, X, X_train=None, Y_train=None, samples=[]):
    X = X.ravel()
    mu = mu.ravel()
    uncertainty = 1.96 * np.sqrt(np.diag(cov))
    
    plt.fill_between(X, mu + uncertainty, mu - uncertainty, alpha=0.1)
    
    plt.plot(X, mu, label='均值')
    for i, sample in enumerate(samples):
        plt.plot(X, sample, lw=1, ls='--', label=f'样本 {i+1}')
    if X_train is not None:
        plt.plot(X_train, Y_train, 'rx')
    plt.legend(bbox_to_anchor=(1.04,0.5), loc="center left")

def plot_gp_2D(gx, gy, mu, X_train, Y_train, title, i):
    ax = plt.gcf().add_subplot(1, 2, i, projection='3d')
    ax.plot_surface(gx, gy, mu.reshape(gx.shape), cmap=cm.coolwarm, linewidth=0, alpha=0.2, antialiased=False)
    ax.scatter(X_train[:,0], X_train[:,1], Y_train, c=Y_train, cmap=cm.coolwarm)
    ax.set_title(title)

%matplotlib inline

# 有限个输入数据点
X = np.arange(-5, 5, 0.2).reshape(-1, 1)

# 先验的均值与方差
mu = np.zeros(X.shape)
cov = kernel(X, X)

# 从先验分布（多元高斯分布）中抽取样本点
samples = np.random.multivariate_normal(mu.ravel(), cov, 3)

# 画出GP的均值, 置信区间
plot_gp(mu, cov, X, samples=samples)

# 倒入计算逆矩阵的函数inv()
from numpy.linalg import inv

def posterior_predictive(X_s, X_train, Y_train, l=1.0, sigma_f=1.0, sigma_y=1e-8):
    '''
    计算后验分布的充分统计量 
    给定 m 个训练数据 X_train 与 Y_train 
    给定 n 个新输入数据 inputs X_s.
    
    Args:
        X_s: 新输入数据 (n x d).
        X_train: 训练输入数据 (m x d).
        Y_train: 训练输出数据 (m x 1).
        l: 核函数的长度参数.
        sigma_f: 核函数的纵向波动参数.
        sigma_y: 噪音参数.
    
    返回:
        后验均值向量 (n x d) 与协方差矩阵 (n x n).
    '''
    K = kernel(X_train, X_train, l, sigma_f) + sigma_y**2 * np.eye(len(X_train))
    K_s = kernel(X_train, X_s, l, sigma_f)
    K_ss = kernel(X_s, X_s, l, sigma_f) + 1e-8 * np.eye(len(X_s))
    K_inv = inv(K)
    
    # 公式 (4)
    mu_s = K_s.T.dot(K_inv).dot(Y_train)

    # 公式 (5)
    cov_s = K_ss - K_s.T.dot(K_inv).dot(K_s)
    
    return mu_s, cov_s

# 无噪音的5个输入数据
X_train = np.array([-4, -3, -2, -1, 1]).reshape(-1, 1)
# y=sin(x)
Y_train = np.sin(X_train)

# 计算后验预测分布的均值向量与协方差矩阵
mu_s, cov_s = posterior_predictive(X, X_train, Y_train)

# 从后验预测分布中抽取3个样本
samples = np.random.multivariate_normal(mu_s.ravel(), cov_s, 3)
plot_gp(mu_s, cov_s, X, X_train=X_train, Y_train=Y_train, samples=samples)

# 定义噪音参数
noise = 0.4

# 带有噪音的训练数据
X_train = np.arange(-3, 4, 1).reshape(-1, 1)
Y_train = np.sin(X_train) + noise * np.random.randn(*X_train.shape)

# 可以对比地看一下带噪音的训练数据与不带噪音的训练数据的区别
plt.figure()
plt.plot(X_train, np.sin(X_train) + noise * np.random.randn(*X_train.shape), lw=1, ls='-.', label='有噪音')
plt.plot(X_train, np.sin(X_train) + 0.0 * np.random.randn(*X_train.shape), lw=2.5, ls='-', label='无噪音')
plt.plot(X_train, np.sin(X_train) + 0.0 * np.random.randn(*X_train.shape), 'rx')
plt.legend(bbox_to_anchor=(1.04,0.5), loc="center left")
plt.show()

# 计算后验预测分布的均值向量以及方差矩阵
mu_s, cov_s = posterior_predictive(X, X_train, Y_train, sigma_y=noise)

# 从后验预测分布中抽取3个样本点
samples = np.random.multivariate_normal(mu_s.ravel(), cov_s, 3)
plot_gp(mu_s, cov_s, X, X_train=X_train, Y_train=Y_train, samples=samples)

params = [
    (0.3, 1.0, 0.2),
    (3.0, 1.0, 0.2),
    (1.0, 0.3, 0.2),
    (1.0, 3.0, 0.2),
    (1.0, 1.0, 0.05),
    (1.0, 1.0, 1.5),
]

plt.figure(figsize=(12, 5))

for i, (l, sigma_f, sigma_y) in enumerate(params):
    mu_s, cov_s = posterior_predictive(X, X_train, Y_train, l=l, 
                                       sigma_f=sigma_f, 
                                       sigma_y=sigma_y)
    plt.subplot(3, 2, i + 1)
    plt.subplots_adjust(top=2)
    plt.title(f'l = {l}, sigma_f = {sigma_f}, sigma_y = {sigma_y}')
    plot_gp(mu_s, cov_s, X, X_train=X_train, Y_train=Y_train)

from numpy.linalg import cholesky, det, lstsq
from scipy.optimize import minimize

def nll_fn(X_train, Y_train, noise, naive=True):
    '''
    基于给定数据X_train和Y_train以及噪声水平
    返回一个可以计算负边际对数似然的函数
    
    Args:
        X_train: 训练输入 (m x d).
        Y_train: 训练输出 (m x 1).
        noise: Y_train的噪声水平.
        naive: 如果 True 那么使用公式(7)来实现
               如果 False 那么使用数值方法来实现
               
        
    返回:
        最小的目标对象
    '''
    def nll_naive(theta):
        # 使用公式(7)来实现 
        # 与下面的nll_stable的实现相比在数值上不稳定
        K = kernel(X_train, X_train, l=theta[0], sigma_f=theta[1]) + \
            noise**2 * np.eye(len(X_train))
        return 0.5 * np.log(det(K)) + \
               0.5 * Y_train.T.dot(inv(K).dot(Y_train)) + \
               0.5 * len(X_train) * np.log(2*np.pi)

    def nll_stable(theta):
        # 数值上更稳定，相比于nll_naive 
        K = kernel(X_train, X_train, l=theta[0], sigma_f=theta[1]) + \
            noise**2 * np.eye(len(X_train))
        L = cholesky(K)
        return np.sum(np.log(np.diagonal(L))) + \
               0.5 * Y_train.T.dot(lstsq(L.T, lstsq(L, Y_train)[0])[0]) + \
               0.5 * len(X_train) * np.log(2*np.pi)
    
    if naive:
        return nll_naive
    else:
        return nll_stable

# 求解可满足最小化目标函数的参数 l 及 sigma_f.
# 实际上，我们应该使用不同的初始化多次运行最小化，以避免局部最小化，
# 但是为了简单起见，此处将其跳过
res = minimize(nll_fn(X_train, Y_train, noise), [1, 1], 
               bounds=((1e-5, None), (1e-5, None)),
               method='L-BFGS-B')

# 将优化结果存储在全局变量中，以便我们以后可以将其与其他实现的结果进行比较
l_opt, sigma_f_opt = res.x

# 使用优化的核函数参数计算后验预测分布的参数，并绘制结果图
mu_s, cov_s = posterior_predictive(X, X_train, Y_train, l=l_opt, sigma_f=sigma_f_opt, sigma_y=noise)
plot_gp(mu_s, cov_s, X, X_train=X_train, Y_train=Y_train)

print(l_opt, sigma_f_opt)

0.9872536793237083 0.8613778055591963

#噪音参数
noise_2D = 0.1

rx, ry = np.arange(-5, 5, 0.3), np.arange(-5, 5, 0.3)
gx, gy = np.meshgrid(rx, rx)

X_2D = np.c_[gx.ravel(), gy.ravel()]

X_2D_train = np.random.uniform(-4, 4, (100, 2))
Y_2D_train = np.sin(0.5 * np.linalg.norm(X_2D_train, axis=1)) + \
             noise_2D * np.random.randn(len(X_2D_train))

plt.figure(figsize=(14,7))

mu_s, _ = posterior_predictive(X_2D, X_2D_train, Y_2D_train, sigma_y=noise_2D)
plot_gp_2D(gx, gy, mu_s, X_2D_train, Y_2D_train, 
           f'参数优化前: l={1.00} sigma_f={1.00}', 1)

res = minimize(nll_fn(X_2D_train, Y_2D_train, noise_2D), [1, 1], 
               bounds=((1e-5, None), (1e-5, None)),
               method='L-BFGS-B')

mu_s, _ = posterior_predictive(X_2D, X_2D_train, Y_2D_train, *res.x, sigma_y=noise_2D)
plot_gp_2D(gx, gy, mu_s, X_2D_train, Y_2D_train,
           f'参数优化后: l={res.x[0]:.2f} sigma_f={res.x[1]:.2f}', 2)

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ConstantKernel, RBF

rbf = ConstantKernel(1.0) * RBF(length_scale=1.0)
gpr = GaussianProcessRegressor(kernel=rbf, alpha=noise**2)

# 1D 训练样本
gpr.fit(X_train, Y_train)

# 计算后验预测分布的均值向量与协方差矩阵
mu_s, cov_s = gpr.predict(X, return_cov=True)

# 获得最优核函数参数
l = gpr.kernel_.k2.get_params()['length_scale']
sigma_f = np.sqrt(gpr.kernel_.k1.get_params()['constant_value'])

# 与前面手写的结果比对
assert(np.isclose(l_opt, l))
assert(np.isclose(sigma_f_opt, sigma_f))

# 绘制结果
plot_gp(mu_s, cov_s, X, X_train=X_train, Y_train=Y_train)