Predict from a latent Gaussian-process model

Computes the approximate posterior of the latent function at new inputs and the predictive distribution of new responses through the likelihood: class probabilities for Bernoulli likelihoods and predictive means and variances for every likelihood.

Usage

predict_latent_gp(
  model,
  x,
  interval_level = 0.95,
  include_covariance = FALSE,
  variance_tolerance = sqrt(.Machine$double.eps),
  exposure = NULL,
  quadrature_order = 40L,
  block_size = NULL
)

Arguments

model

A fitted gaussianprocesses_latent_model from fit_latent_gp() or optimize_latent_gp().

x

Numeric vector or matrix of prediction inputs.

interval_level

Probability level of the latent credible intervals.

include_covariance

If TRUE, also return the full latent posterior covariance matrix.

variance_tolerance

Relative tolerance used when deciding whether a small negative posterior variance is floating-point error or a material numerical failure.

exposure

Optional positive exposure of the new responses for Poisson likelihoods: one value or one per prediction input. It defaults to 1.

quadrature_order

Number of Gauss-Hermite nodes for predictive integrals without a closed form (likelihood_predictive()).

block_size

Number of prediction inputs processed together when include_covariance = FALSE, as in predict_gp().

Value

An object of class gaussianprocesses_latent_prediction with the shared latent fields of predict_gp(), x, mean, latent_variance, latent_sd, interval_level, and latent_interval (and latent_covariance with include_covariance = TRUE), and the predictive distribution of a new response: response_mean, response_variance, probability \(\Pr(y_* = 1)\) and probability_interval for Bernoulli likelihoods, rate_interval, prediction_interval, and exposure for Poisson likelihoods, and prediction_interval for Gaussian likelihoods. A non-Gaussian response has no observation-noise variance, so the observation fields of predict_gp() are replaced by these; prediction_interval has its meaning there, an interval for a new response.

Details

Under the Laplace approximation the latent values at new inputs are Gaussian (Rasmussen and Williams, 2006, Algorithm 3.2), with $$\mathrm{E}[f_*] = m(x_*) + k_*^\top \nabla \log p(y \mid \hat f), \qquad \mathrm{Var}[f_*] = k(x_*, x_*) - v^\top v,$$ where \(k_* = k(X, x_*)\), \(v = R^{-\top} W^{1/2} k_*\), and \(R^\top R = B\). The predictive distribution of a new response integrates the likelihood against this Gaussian with likelihood_predictive().

Class probabilities

For a Bernoulli likelihood, probability is \(\Pr(y_* = 1 \mid y)\) for the response coded 1: the second level of a factor response (model$response_levels), or TRUE. It averages the link over the latent posterior, $$\Pr(y_* = 1 \mid y) = \int \pi(f_*) N(f_* \mid \mu_*, \sigma_*^2) df_*,$$ which for the probit link is \(\Phi(\mu_* / \sqrt{1 + \sigma_*^2})\) exactly and for the logistic link is computed by Gauss-Hermite quadrature. It is not \(\pi(\mu_*)\): the plug-in probability ignores the latent uncertainty and is too confident, by up to 0.04 at \(\sigma_*^2 = 1\) and 0.11 at \(\sigma_*^2 = 4\). MacKay's approximation \(\sigma(\kappa \mu_*)\) with \(\kappa = (1 + \pi \sigma_*^2 / 8)^{-1/2}\) is closer but still off by up to 0.005 at \(\sigma_*^2 = 1\) and 0.01 at \(\sigma_*^2 = 4\), so it is not used.

probability_interval describes the uncertainty of the class probability \(\pi(f_*)\) itself: because the link is increasing, the interval_level quantiles of \(\pi(f_*)\) are the link applied to the latent interval, exactly under the Gaussian approximation. Joint uncertainty, for example of a difference of probabilities, comes from the probability draws of sample_gp_posterior().

Multi-class (softmax) classification is not supported.

Counts

For a Poisson likelihood the rate \(\lambda_* = E_* e^{f_*}\) is lognormal: its mean is \(E_* e^{\mu_* + \sigma_*^2 / 2}\) (response_mean, also the mean of the count) and its quantiles are \(E_* e^{\mu_* \pm z \sigma_*}\) exactly, because the exponential is increasing (rate_interval). A new count has the Poisson-lognormal predictive distribution; prediction_interval gives its equal-tailed quantiles at interval_level, computed by quadrature (gp_count_scores()). Counts are discrete, so these intervals cover at least interval_level: their coverage is conservative.

Stability

Experimental: this interface may change in a minor release, and every change is listed in NEWS. See gaussianprocesses-package for the policy.

See also

fit_latent_gp(), sample_gp_posterior().

Examples

x <- seq(-3, 3, length.out = 30)
y <- as.numeric(sin(x) + 0.3 * cos(5 * x) > 0)
model <- fit_latent_gp(x, y, rbf_kernel(variance = 4), bernoulli_likelihood())

prediction <- predict_latent_gp(model, c(-2, 0, 2, 6))

# Far from the data the latent variance grows and the probability moves
# towards 1/2, unlike the logistic function of the latent mean.
cbind(prediction$probability, plogis(prediction$mean))
#>           [,1]       [,2]
#> [1,] 0.0832661 0.04920686
#> [2,] 0.6497020 0.67195584
#> [3,] 0.8880144 0.92368837
#> [4,] 0.4972552 0.49546867