机器学习中的高斯过程简介-好文

 

 

 

 

import matplotlib.pyplot as plt
import numpy as np
from itertools import cycle

color_cycle = cycle('kbryg')
n_variable = 20
n_sample = 5

plt.figure(figsize=(500, 300))
sigma_s = np.eye(n_variable) # 协方差阵为对角矩阵

for _ in range(n_sample):
    point = np.random.multivariate_normal(np.zeros(n_variable), sigma_s)
    plt.plot(np.arange(n_variable), point, color=next(color_cycle))
    plt.scatter(np.arange(n_variable), point)

plt.rcParams['font.sans-serif'] = ['SimHei']  # 正确显示中文
plt.rcParams['axes.unicode_minus'] = False  # 正确显示负号
plt.title('不使用核函数的情况', fontsize=24)
plt.xticks(np.arange(n_variable))
plt.xlabel('y', fontsize=20)
plt.show()

 

核函数来估计协方差，能很好反映各个变量y yy间的关系。因为核函数充分考虑了样本间的相似性(大多用d dd度量)。
例如采用squared exponential核函数后，曲线变得平滑了：
 

xs = np.linspace(0, 1, n_variable)  # 构建y对应的特征
sigma_s = np.exp(-(np.expand_dims(xs, axis=0) - np.expand_dims(xs, axis=1)) ** 2 / 2) 
 # 依据特征xs构建协方差阵

 

import numpy as np
from scipy.optimize import minimize
from scipy.spatial.distance import pdist, cdist, squareform


class Kernel:
    # RBF核函数
    def __init__(self, theta):
        self.theta = theta

    def __call__(self, X, Y=None):
        if Y is None:
            dists = pdist(X / self.theta, metric='sqeuclidean')
            K = np.exp(-0.5 * dists)
            K = squareform(K)
            np.fill_diagonal(K, 1)
        else:
            dists = cdist(X / self.theta, Y / self.theta, metric='sqeuclidean')
            K = np.exp(-0.5 * dists)
        return K


class GPR:
    def __init__(self):
        self.K = None
        self.X = None
        self.y = None

    def log_marginal_likelihood(self, theta):
        # 计算对数边缘似然函数
        K = Kernel(theta)
        sigma = K(self.X)
        log_likelihood = np.log(np.linalg.det(sigma)) + \
                         self.y @ np.linalg.inv(sigma) @ self.y + \
                         sigma.shape[0] * np.log(2 * np.pi)
        return - 0.5 * log_likelihood

    def fit(self, X, y):
        # 训练，获取核函数最优参数
        self.X = X
        self.y = y

        def obj_func(theta):
            return - self.log_marginal_likelihood(theta)

        theta_opt = minimize(obj_func, np.array([1.0]), method='BFGS')
        self.K = Kernel(theta_opt.x[0])

    def predict(self, X_pred):
        # 预测，获取y*的条件概率分布
        K_pred = self.K(X_pred)
        K_train = self.K(self.X)
        K_pred_train = self.K(self.X, X_pred)
        K_inv = np.linalg.inv(K_train)
        mu = K_pred_train.T @ K_inv @ self.y
        sigma = K_pred - K_pred_train.T @ K_inv @ K_pred_train
        return mu, np.diagonal(sigma)


if __name__ == '__main__':
    import matplotlib.pyplot as plt

    gpr = GPR()

    # 生成样本数据
    coefs = [6, -2.5, -2.4, -0.1, 0.2, 0.03]


    def f(x):
        total = 0
        for exp, coef in enumerate(coefs):
            total += coef * (x ** exp)
        return total


    xs = np.linspace(-5.0, 3.5, 100)
    ys = f(xs)

    X_train = np.array([-4, -1.5, 0, 1.5, 2.5, 2.7])
    y_train = f(X_train)

    X_train = X_train.reshape(-1, 1)
    X_pred = np.linspace(-8, 7, 80).reshape((-1, 1))

    gpr.fit(X_train, y_train)
    y_pred, y_std = gpr.predict(X_pred)

    plt.plot(xs, ys, color='k', linewidth=2, label='True')
    plt.scatter(X_train, y_train, color='b', marker='*', linewidths=3, label='Train_data')
    plt.plot(X_pred, y_pred, color='r', label='Pred')

    plt.fill_between(X_pred.reshape(1, -1)[0], y_pred - y_std, y_pred + y_std, color='darkorange',
                     alpha=0.2)
    plt.legend()
    plt.show()

 

参考资料
http://bridg.land/posts/gaussian-processes-1
http://www.datalearner.com/blog/1051459170229238
https://en.wikipedia.org/wiki/Gaussian_process
--------------------- 
作者：slx_share 
来源：CSDN 
原文：https://blog.csdn.net/slx_share/article/details/83573339 
版权声明：本文为博主原创文章，转载请附上博文链接！