高斯后验过程（Python）答案

【问题标题】：Gaussian Process Posterior (Python)高斯后验过程（Python）
【发布时间】：2018-05-11 01:31:21
【问题描述】：

我使用以下代码创建并采样了均值 = 0 的联合高斯先验：

import numpy as np
import matplotlib.pyplot as plt 
from math import pi 
from scipy.spatial.distance import cdist
import scipy.stats as sts

x_prior = np.linspace(-10,10,101)
x_prior = x_prior.reshape(-1,1)
mu = np.zeros(x_prior.shape)

#defining the Kernel for the covariance function

def sec(a,b, length_scale , sigma) : 
    K = sigma * np.exp(-1/(2*length_scale) * cdist(a,b)**2)
    return K 

#defining the Gaussian Process prior

def GP(a , b, mu , kernel , length_scale, sigma , samples ) :
    f = np.random.multivariate_normal(mu.flatten(), kernel(a ,b , length_scale , sigma ) , samples)
    return f

prior = GP(x_prior ,x_prior, mu , sec , 100, 1 , 5)

plt.figure()
plt.grid()
plt.title('samples from the Gaussian prior')
plt.plot(x_prior , prior.T)
plt.show()

然后，当添加一些“观察到的”数据时，我希望计算这些点的后验，但这就是我卡住的地方。

这是我引入新数据的代码：

x_train = np.array([-10,-8,5,-1,2])
x_train = x_train.reshape(-1,1)
def straight_line(m , x , c):
    y = 5*x + c
    return y
ytrain = straight_line(5 , x_train , 0)

据我了解，在给定与观测数据相关的先前和新 x 值的情况下，您计算新数据的条件分布。

然后，您是否希望通过对平均值进行某种更改以包含新的 y 值来更新多元变量成为后验？

我已使用以下资源进行尝试：

http://katbailey.github.io/post/gaussian-processes-for-dummies/ https://www.robots.ox.ac.uk/~mebden/reports/GPtutorial.pdf

但我真的很想了解每个阶段会发生什么，以及为什么，所以当我得到一个后验（我做不到）时，我确切地知道我是如何到达那里的。

以下是我一直在尝试实施但至今无济于事的一些解决方案：

K_train = sec(x_train , x_train , 1,1)
K_prior = sec(x_prior , x_prior , 1,1)
K_pt =  sec(x_prior , x_train , 1,1)
K_tp = sec(x_train , x_prior ,  1,1)  ## = k_tp transpose
prior = sts.multivariate_normal(mu.flatten(), K_prior) 
#mean_test = np.dot(K_p , np.linalg.inv(K_prior))
mean_function = np.dot(np.dot(K_tp ,np.linalg.inv(K_prior).T) , prior )
covariance_function = K_train - np.dot(np.dot(K_tp ,np.linalg.inv(K_prior).T) , K_pt)

【问题讨论】：

标签： python machine-learning process gaussian sampling

【解决方案1】：

只是为了进一步跟进。我在这里把我的代码写成了 Juypiter 格式：

https://github.com/SpaceMeerkat/Scariff

这里有相关的通读材料：

https://spacemeerkat.wordpress.com/

以防万一有人想通过这种材料工作并像我一样陷入困境。

【讨论】：

【解决方案2】：

只是为看到此内容的任何人提供的更新。我找到了阅读本文的解决方案：

https://arxiv.org/pdf/1711.10834.pdf

以及以下代码：

mean_function = np.dot(np.dot(K_pt ,np.linalg.inv(K_train)), ytrain) 

covariance_function = K_prior - np.dot(np.dot(K_pt ,np.linalg.inv(K_train)) , K_tp) 

f = np.random.multivariate_normal(mean_function[:,0],covariance_function , 100)

其中 f 是您从中采样的后关节高斯分布

【讨论】：

【解决方案3】：

观察数据 x(1:N), y(1:N) 后，使用高斯过程对新点 x* 的预测（克里金法）具有以下形式：

下面的代码显示了上述贝叶斯更新方程的实现，以在给定先验和观察数据的情况下计算后验（这里蓝色星代表训练数据点，红线表示相应的预测，GP 和绿色带是置信区间):

x_train = np.linspace(-10, 4, 10).reshape(-1,1)
y_train = np.random.random(10)
x_p = np.linspace(-10, 4, 50).reshape(-1,1) 
K_train = kernel(x_train, x_train, length_scale=2, sigma=1)
K_pt = kernel(x_p, x_train, length_scale=2, sigma=1)
K_tp = kernel(x_train, x_p, length_scale=2, sigma=1)
K_prior = kernel(x_p, x_p, length_scale=2, sigma=1)
# compute posterior
mean_function = np.dot(np.dot(K_pt ,np.linalg.inv(K_train)), y_train) 
covariance_function = K_prior - np.dot(np.dot(K_pt ,np.linalg.inv(K_train)) , K_tp) 
plt.plot(x_train, y_train, '*')
plt.plot(x_p, mean_function)
plt.fill_between(x_p.ravel(), mean_function-3*np.sqrt(np.diag(covariance_function)), mean_function+3*np.sqrt(np.diag(covariance_function)), color='g', alpha=.2)

预测后验分布如下图所示：

【讨论】：