python - 高斯过程后验(Python)

我使用下面的代码创建并采样了一个平均值为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)

最佳答案

只是为了进一步跟进。我已经在这里用juypter格式编写了代码：
https://github.com/SpaceMeerkat/Scariff
相关阅读材料如下：
https://spacemeerkat.wordpress.com/
以防万一有人想通过这种材料，并成为卡住像我一样。

关于python - 高斯过程后验(Python)，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/47517230/