问题
将GP与不同的似然函数结合可以构建非常灵活的模型，然而这会使得模型的推理变得不可解（intractable），因为似然函数不再是高斯过程的共轭。所以我们需要变分推理来近似f的后验分布或者MCMC方法从f后验分布来采样，从而预测在新的点的函数值。
GPflow.models.GPMC和gpflow.train.HMC的结合使用就实现了MCMC方法。
模型
从完全的贝叶斯角度看，一般GP模型的数据生成过程可以表示为：
\theta \sim p(\theta)
f \sim \mathcal {GP}\Big(m(x; \theta), k(x, x'; \theta)\Big)
y_i \sim p\Big(y | g(f(x_i)\Big)
首先从\theta的先验分布采样一个\theta，然后从高斯分布得到一组隐函数f，然后再由一个连接函数g(\cdot)映射到观测函数y。
模型推理
首先明确我们要求的是f_\ast的后验分布
p(f_\ast | \bm{x}_\ast, \bm{X}, \bm{y}) = \int p(f_\ast | \bm{x}_\ast, \bm{X}, \bm{f}) p(\bm{f} | \bm{X}, \bm{y}) ~d{\bm{f}} (1)
如果我们将隐函数f和模型超参数\theta都考虑乘模型参数，(1)可以改写为
p(f_\ast | \bm{x}_\ast, \bm{X}, \bm{y}) = \int p(f_\ast | \bm{x}_\ast, \bm{X}, \bm{f}, \theta) p(\bm{f}, \theta | \bm{X}, \bm{y}) ~d{\bm{f}} ~d \theta (2)
我们只需使用Hamiltonian Monte Carlo (HMC) 从p(\bm{f}, \theta | \bm{X}, \bm{y})联合采样f和\theta即可。利用贝叶斯法则，上式改写为p(\bm{f}, \theta | \bm{X}, \bm{y}) = \frac {p(\bm{y} | \bm{X}, \bm{f}, \theta) p(\bm{f}) p(\theta)} {p(\bm{y}|\bm{X})}。分子部分对于f和\theta是常数，只需要从分子部分采样即可。
例子1
准备数据
import gpflow
from gpflow.test_util import notebook_niter
import numpy as np
import matplotlib
%matplotlib inline
matplotlib.rcParams['figure.figsize'] = (12, 6)
plt = matplotlib.pyplot
X = np.linspace(-3,3,20)
Y = np.random.exponential(np.sin(X)**2)
with gpflow.defer_build():
k = gpflow.kernels.Matern32(1, ARD=False) + gpflow.kernels.Bias(1)
l = gpflow.likelihoods.Exponential()
m = gpflow.models.GPMC(X[:,None], Y[:,None], k, l)
m.kern.kernels[0].lengthscales.prior = gpflow.priors.Gamma(1., 1.)
m.kern.kernels[0].variance.prior = gpflow.priors.Gamma(1., 1.)
m.kern.kernels[1].variance.prior = gpflow.priors.Gamma(1., 1.)
用AdamOptimizer先找一个较好的初始点
m.compile()
o = gpflow.train.AdamOptimizer(0.01)
o.minimize(m, maxiter=notebook_niter(15)) # start near MAP
采样
s = gpflow.train.HMC()
samples = s.sample(m, notebook_niter(500), epsilon=0.12, lmax=20, lmin=5, thin=5, logprobs=False)#, verbose=True)
求平均
xtest = np.linspace(-4,4,100)[:,None]
f_samples = []
for i, s in samples.iterrows():
f = m.predict_f_samples(xtest, 5, initialize=False, feed_dict=m.sample_feed_dict(s))
f_samples.append(f)
f_samples = np.vstack(f_samples)
画图
rate_samples = np.exp(f_samples[:, :, 0])
line, = plt.plot(xtest, np.mean(rate_samples, 0), lw=2)
plt.fill_between(xtest[:,0],
np.percentile(rate_samples, 5, axis=0),
np.percentile(rate_samples, 95, axis=0),
color=line.get_color(), alpha = 0.2)
plt.plot(X, Y, 'kx', mew=2)
plt.ylim(-0.1, np.max(np.percentile(rate_samples, 95, axis=0)))
结果如下图，黑色散点为要拟合的点，蓝色线为预测函数，浅蓝色带为5%和95%分位数位置。
下图是由采样的variance值所画的直方图，代表variance的后验分布。
代码解读
p(\bm{y} | \bm{X}, \bm{f}, \theta)对应GPMC._build_likelihood()。
p(\bm{f})对应GPMC.__init__()中的self.V.prior = Gaussian(0., 1.)。
p(\theta)对应kernel和likelihood中的其他参数的先验分布。
至此，GPMC模块也讲完了。
