问题

将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%分位数位置。

MCMC得到拟合函数（蓝色曲线）

下图是由采样的variance值所画的直方图，代表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模块也讲完了。