高斯過程回歸 （GPR）

首先定義一個(gè)輸入空間X，定義一個(gè)函數(shù)f，它將X上的點(diǎn)映射到空間F。

F上的每個(gè)點(diǎn)都是一個(gè)隨機(jī)變量，GPR假設(shè)F上的點(diǎn)服從高斯過程，即對(duì)于任意有限個(gè)點(diǎn)f_1, ..., f_n，他們的聯(lián)合分布都是一個(gè)高斯分布。其中均值由均值函數(shù)定義，協(xié)方差矩陣由協(xié)方差函數(shù)定義。

GPR還假設(shè)了一個(gè)一維的高斯噪音\epsilon，它也是服從高斯分布的一維隨機(jī)變量。f+\epsilon即表示觀測到的樣本點(diǎn)y，所有樣本點(diǎn)的集合稱為Y。由于f無法實(shí)際觀測到，因此我們把f稱為隱函數(shù)，y稱為觀測函數(shù)。

以上就是GPR的數(shù)據(jù)生成過程。

那么對(duì)于給定一個(gè)新的輸入x_\ast，如何預(yù)測對(duì)應(yīng)的f_\ast或 y_\ast呢？

首先明確我們要求的是f_\ast的后驗(yàn)分布

p(f_\ast | \bm{x}_\ast, \bm{X}, \bm{y}) = \mathcal{N} (\mu, \sigma^2) (1)

其中\mu=\bm{k_\ast}^T (K+\sigma_n^2 I)^{-1} \bm{y}， \sigma^2 = k_{\ast\ast} - \bm{k_\ast}^T (K + \sigma^2 I)^{-1}\bm{k_\ast}。推理過程可以參考GPML的17頁。

注意(1)中包含了GPR模型的超參數(shù)，我們采用什么原則優(yōu)化超參數(shù)呢？

這里要說一下，GPML書中寫道，"The posterior in eq. (2.5) combines the likelihood and the prior, and captures everything we know about the parameters."，也就是說后驗(yàn)包含了我們知道的所有信息，包括先驗(yàn)和數(shù)據(jù)。
根據(jù)GPML的第五章model selection的109頁，有三個(gè)級(jí)別的參數(shù)優(yōu)化策略。第一層是最大化模型隱函數(shù)f的后驗(yàn)概率，第二層是最大化模型超參數(shù)\theta的后驗(yàn)概率，最后是將模型類型也當(dāng)做變量，最大化模型H的后驗(yàn)概率。

實(shí)際上我們采用的是第二層策略 (5.5)，最大化超參數(shù)\theta的后驗(yàn)概率。但是由于分母無法計(jì)算且是關(guān)于\theta的常數(shù)，因此最大化分子部分。分子部分是邊際似然函數(shù)\log p(\bm{y}|X)和超參數(shù)先驗(yàn)分布p(\theta)的乘積。

\log p(\bm{y}|X)=-\frac{1}{2}\bm{y}^T(K+\sigma_n^2I)^{-1}\bm{y} - \frac{1}{2}\log|K+\sigma^2I|-\frac{n}{2}\log2\pi (2)

這個(gè)策略又叫第二類型最大似然估計(jì)（ML-II）。

例子1

import gpflow
import numpy as np
import matplotlib
%matplotlib inline
matplotlib.rcParams['figure.figsize'] = (12, 6)
plt = matplotlib.pyplot

準(zhǔn)備數(shù)據(jù)

N = 12
X = np.random.rand(N,1)
Y = np.sin(12*X) + 0.66*np.cos(25*X) + np.random.randn(N,1)*0.1 + 3
plt.plot(X, Y, 'kx', mew=2)

構(gòu)建模型

k = gpflow.kernels.Matern52(1, lengthscales=0.3)
m = gpflow.models.GPR(X, Y, kern=k)
m.likelihood.variance = 0.01

k = gpflow.kernels.Matern52(1, lengthscales=0.3)
meanf = gpflow.mean_functions.Linear(1.0, 0.0)
m = gpflow.models.GPR(X, Y, k, meanf)
m.likelihood.variance = 0.01

優(yōu)化模型

gpflow.train.ScipyOptimizer().minimize(m)
plot(m)
m.as_pandas_table()

在新的點(diǎn)上預(yù)測函數(shù)值

xx = np.linspace(-0.1, 1.1, 100).reshape(100, 1)
mean, var = m.predict_y(xx)

源代碼解讀

在執(zhí)行gpflow.train.ScipyOptimizer().minimize(m)時(shí)發(fā)生了什么？

首先模型初始化、計(jì)算目標(biāo)函數(shù)。GPR._build_likelihood()計(jì)算的是gaussian log marginal likelihood，所以GPR的目標(biāo)函數(shù)實(shí)際上是對(duì)數(shù)邊際似然+對(duì)數(shù)參數(shù)先驗(yàn)，見下面代碼：

  def build_objective(self):
        likelihood = self._build_likelihood()
        priors = []
        for param in self.parameters:
            unconstrained = param.unconstrained_tensor
            constrained = param._build_constrained(unconstrained)
            priors.append(param._build_prior(unconstrained, constrained))
        prior = self._build_prior(priors)
        return self._build_objective(likelihood, prior)

    def _build_objective(self, likelihood_tensor, prior_tensor):
        func = tf.add(likelihood_tensor, prior_tensor, name='nonneg_objective')
        return tf.negative(func, name='objective')  # 最小化的是負(fù)對(duì)數(shù)似然函數(shù)

其次構(gòu)建tensorflow計(jì)算圖，然后再利用L-BFGS-B優(yōu)化目標(biāo)函數(shù)。

在執(zhí)行m.predict_y(xx)時(shí)發(fā)生了什么？

調(diào)用GPR._build_predict()，根據(jù)GPML Regression一章的Algorithm2.1計(jì)算預(yù)測函數(shù)值。

例子2

例子1利用最大化\theta的似然函數(shù)（也就是f的邊際似然函數(shù)）來求得\theta，然后通過解析解得到預(yù)測函數(shù)值。我們也可以用MCMC方法來預(yù)測函數(shù)值。這里我們要求的是p(\theta | \bm{X}, \bm{y})，也可以寫成p(f_\ast | \bm{x}_\ast, \bm{X}, \bm{y}) = \int p(f_\ast | \bm{x}_\ast, \bm{X}, \bm{y}, \theta) p(\theta) ~d\theta。我們用HMC從先驗(yàn)p(\theta)采樣N個(gè)點(diǎn)，然后用\sum_i^N p(\bm{y} | \bm{X}, \theta) * p(\theta_i)來近似。

代碼

設(shè)置超參數(shù)的先驗(yàn)

m.clear()
m.kern.lengthscales.prior = gpflow.priors.Gamma(1., 1.)
m.kern.variance.prior = gpflow.priors.Gamma(1., 1.)
m.likelihood.variance.prior = gpflow.priors.Gamma(1., 1.)
m.mean_function.A.prior = gpflow.priors.Gaussian(0., 10.)
m.mean_function.b.prior = gpflow.priors.Gaussian(0., 10.)
m.compile()
m.as_pandas_table()

采樣

sampler = gpflow.train.HMC()
samples = sampler.sample(m, num_samples=500, epsilon=0.05, lmin=10, lmax=20, logprobs=False)

求平均

#plot the function posterior
xx = np.linspace(-0.1, 1.1, 100)[:,None]
plt.figure(figsize=(12, 6))
for i, s in samples.iloc[::20].iterrows():
    f = m.predict_f_samples(xx, 1, initialize=False, feed_dict=m.sample_feed_dict(s))
    plt.plot(xx, f[0,:,:], 'C0', lw=2, alpha=0.1)
    
plt.plot(X, Y, 'kx', mew=2)
_ = plt.xlim(xx.min(), xx.max())
_ = plt.ylim(0, 6)

至此，GPflow中的GPR就講完了。

參考文獻(xiàn)

[1] GPflow-master/doc/source/notebooks/regression.py
[2] Rasmussen, Carl Edward. "Gaussian processes in machine learning." Advanced lectures on machine learning. Springer, Berlin, Heidelberg, 2004. 63-71.