GP-Marginal.ipynb

https://github.com/pymc-devs/pymc-examples/blob/main/examples/gaussian_processes/GP-Marginal.ipynb

(gp_marginal)=

Marginal Likelihood Implementation

:::{post} June 4, 2023 :tags: gaussian process, time series :category: reference, intermediate :author: Bill Engels, Chris Fonnesbeck :::

The gp.Marginal class implements the more common case of GP regression: the observed data are the sum of a GP and Gaussian noise. gp.Marginal has a marginal_likelihood method, a conditional method, and a predict method. Given a mean and covariance function, the function $f(x)$ is modeled as,

f(x)∼GP(m(x), k(x,x′)) .f(x) \sim \mathcal{GP}(m(x),\, k(x, x')) \,.

The observations $y$ are the unknown function plus noise

ϵ∼N(0,Σ)y=f(x)+ϵ\begin{aligned} \epsilon &\sim N(0, \Sigma) \\ y &= f(x) + \epsilon \\ \end{aligned}

The .marginal_likelihood method

The unknown latent function can be analytically integrated out of the product of the GP prior probability with a normal likelihood. This quantity is called the marginal likelihood.

p(y∣x)=∫p(y∣f,x) p(f∣x) dfp(y \mid x) = \int p(y \mid f, x) \, p(f \mid x) \, df

The log of the marginal likelihood, $p(y \mid x)$, is

log⁡p(y∣x)=−12(y−mx)T(Kxx+Σ)−1(y−mx)−12log⁡(Kxx+Σ)−n2log⁡(2π)\log p(y \mid x) = -\frac{1}{2} (\mathbf{y} - \mathbf{m}_x)^{T} (\mathbf{K}_{xx} + \boldsymbol\Sigma)^{-1} (\mathbf{y} - \mathbf{m}_x) - \frac{1}{2}\log(\mathbf{K}_{xx} + \boldsymbol\Sigma) - \frac{n}{2}\log (2 \pi)

$\boldsymbol\Sigma$ is the covariance matrix of the Gaussian noise. Since the Gaussian noise doesn't need to be white to be conjugate, the marginal_likelihood method supports either using a white noise term when a scalar is provided, or a noise covariance function when a covariance function is provided.

The gp.marginal_likelihood method implements the quantity given above. Some sample code would be,

python

The .conditional distribution

The .conditional has an optional flag for pred_noise, which defaults to False. When pred_sigma=False, the conditional method produces the predictive distribution for the underlying function represented by the GP. When pred_sigma=True, the conditional method produces the predictive distribution for the GP plus noise. Using the same gp object defined above,

python

If using an additive GP model, the conditional distribution for individual components can be constructed by setting the optional argument given. For more information on building additive GPs, see the main documentation page. For an example, see the Mauna Loa CO$_2$ notebook.

Making predictions

The .predict method returns the conditional mean and variance of the gp given a point as NumPy arrays. The point can be the result of find_MAP or a sample from the trace. The .predict method can be used outside of a Model block. Like .conditional, .predict accepts given so it can produce predictions from components of additive GPs.

python

Example: Regression with white, Gaussian noise

The estimated values are close to their true values.

Using .conditional

The prediction also matches the results from gp.Latent very closely. What about predicting new data points? Here we only predicted $f_$, not $f_$ + noise, which is what we actually observe.

The conditional method of gp.Marginal contains the flag pred_noise whose default value is False. To draw from the posterior predictive distribution, we simply set this flag to True.

Notice that the posterior predictive density is wider than the conditional distribution of the noiseless function, and reflects the predictive distribution of the noisy data, which is marked as black dots. The light colored dots don't follow the spread of the predictive density exactly because they are a single draw from the posterior of the GP plus noise.

Using .predict

We can use the .predict method to return the mean and variance given a particular point.

Authors

Watermark