問題
將GP與不同的似然函數(shù)結合可以構建非常靈活的模型惯疙，然而這會使得模型的推理變得不可解（intractable）土砂，因為似然函數(shù)不再是高斯過程的共軛祝谚。所以我們需要變分推理來近似f的后驗分布或者MCMC方法從f后驗分布來采樣，從而預測在新的點的函數(shù)值们豌。
GPflow.models.GPMC和gpflow.train.HMC的結合使用就實現(xiàn)了MCMC方法涯捻。
模型
從完全的貝葉斯角度看，一般GP模型的數(shù)據(jù)生成過程可以表示為：
\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望迎，然后從高斯分布得到一組隱函數(shù)f障癌，然后再由一個連接函數(shù)g(\cdot)映射到觀測函數(shù)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)
如果我們將隱函數(shù)f和模型超參數(shù)\theta都考慮乘模型參數(shù)辩尊，(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})聯(lián)合采樣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是常數(shù)轿亮，只需要從分子部分采樣即可。
例子1
準備數(shù)據(jù)
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)))
結果如下圖蒿涎，黑色散點為要擬合的點哀托，藍色線為預測函數(shù)，淺藍色帶為5%和95%分位數(shù)位置劳秋。
下圖是由采樣的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中的其他參數(shù)的先驗分布。
至此补履，GPMC模塊也講完了添坊。
