Treatments of mean-function coefficients

A parametric mean $m(x) = h(x)^\top \beta$ has coefficients $\beta$ that are fixed, estimated, or marginalized. Pass the result of these functions as the coefficients argument of linear_mean(), polynomial_mean(), basis_mean(), or as the value of constant_mean(); a numeric vector fixes the coefficients instead.

Usage

estimate_coefficients(method = c("ml", "reml"))

coefficient_prior(mean = 0, covariance = NULL)

Arguments

method: "ml" or "reml": the likelihood that log_marginal_likelihood() reports and optimize_gp() maximizes.
mean: Prior mean $b$ of the coefficients: one value, recycled, or one per coefficient.
covariance: Prior covariance $B$: NULL for the vague limit $B^{-1} \to 0$; a positive number or vector, for a diagonal covariance; or a symmetric positive-definite matrix.

Value

An object of class gaussianprocesses_coefficients.

Details

Write $H$ for the matrix of basis rows $h(x_i)^\top$, $C = K + \Sigma$ for the observation covariance, and $p$ for the number of coefficients.

Estimated. The coefficients are the generalized least-squares estimate $\hat\beta = (H^\top C^{-1} H)^{-1} H^\top C^{-1} y$, which maximizes the likelihood for given hyperparameters. Predictions plug in $\hat\beta$ and do not include its uncertainty. With method = "ml" the likelihood is the profile likelihood, the Gaussian likelihood at $\hat\beta$. With method = "reml" it is the restricted likelihood $$\ell_R = -\tfrac12 y^\top P y - \tfrac12 \log|C| - \tfrac12 \log|H^\top C^{-1} H| - \tfrac{n - p}{2} \log 2\pi,$$ with $P = C^{-1} - C^{-1} H (H^\top C^{-1} H)^{-1} H^\top C^{-1}$, which accounts for the $p$ degrees of freedom used by $\hat\beta$ and so estimates variances with less bias.

Marginalized. With the prior $\beta \sim N(b, B)$, the model is exactly a Gaussian process with mean $h(x)^\top b$ and covariance $k(x, x') + h(x)^\top B h(x')$ (Rasmussen and Williams, 2006, section 2.7). Predictions use the posterior mean $\bar\beta = (B^{-1} + H^\top C^{-1} H)^{-1}(H^\top C^{-1} y + B^{-1} b)$ and add the coefficient uncertainty $R_* (B^{-1} + H^\top C^{-1} H)^{-1} R_*^\top$, $R_* = H_* - K_{*X} C^{-1} H$, to the predictive covariance; posterior draws include it too. The likelihood is the marginal likelihood of that Gaussian process. In the vague limit $\bar\beta = \hat\beta$ and the likelihood, up to a constant, is the restricted likelihood above.

Rank-deficient bases make estimated coefficients, and those with a vague prior, unidentifiable, and raise an error of class gaussianprocesses_rank_deficiency_error. A proper prior identifies them.

Stability

Stable: from version 1.0.0 this interface changes incompatibly only in a major release, after a deprecation period. Results and options that concern an experimental model class, kernel, or argument follow that interface's tier. See gaussianprocesses-package for the policy.

References

Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press. Section 2.7.

Examples

x <- seq(0, 4, length.out = 25)
y <- 1 + 0.5 * x + sin(2 * x)

estimated <- fit_gp(x, y, rbf_kernel(length_scale = 0.5),
  noise_variance = 0.01, mean = linear_mean())
coef(estimated)
#> intercept        x1 
#> 0.9194838 0.6207013 

# A vague prior gives the same coefficients and adds their uncertainty to
# predictions away from the data.
marginalized <- fit_gp(x, y, rbf_kernel(length_scale = 0.5),
  noise_variance = 0.01, mean = linear_mean(coefficient_prior()))
predict_gp(estimated, 8)$latent_sd
#> [1] 1
predict_gp(marginalized, 8)$latent_sd
#> [1] 2.22785