Uncertainty of estimated hyperparameters

Approximates the uncertainty of a model's hyperparameters from the curvature of the log marginal likelihood at its optimum, and flags hyperparameters that the data identify poorly.

Usage

gp_hyperparameter_uncertainty(
  model,
  step = 1e-04,
  interval_level = 0.95,
  condition_threshold = 1e+06,
  log_scale_threshold = log(2)
)

Arguments

model: A fitted gaussianprocesses_model, normally from optimize_gp().
step: Step of the central differences of the gradient, on the optimizer's coordinate scale.
interval_level: Level of the approximate intervals.
condition_threshold: Condition number of the observed information above which the hyperparameters are flagged as poorly identified.
log_scale_threshold: Standard deviation, on the log scale, of the least determined combination of log-scale parameters above which they are flagged as poorly identified. The default, $\log 2$, flags a combination known only to within a factor of two at one standard deviation.

Value

An object of class gaussianprocesses_hyperparameter_uncertainty with

parameters: Data frame: parameter, estimate, coordinate ("log" or "identity"), coordinate_estimate, standard_error on the coordinate scale, and lower and upper on the natural scale.
covariance, correlation: On the coordinate scale; NULL when the information is not positive definite.
information, eigenvalues, condition_number, weakest_log_scale_sd, poorly_identified, weakest_direction: The observed information and its diagnostics.
gradient: The gradient at the estimate.

Details

The hyperparameters are those optimize_gp() estimated (all kernel parameters and the noise for other models), on its coordinate scale $\eta$: logarithms of positive parameters and real parameters themselves. The observed information $I = -\nabla^2 \ell(\hat\eta)$ is computed by central differences of the analytical gradient (log_marginal_likelihood_gradient()), one refit per displaced coordinate. Its inverse is the covariance of the Laplace approximation with a flat prior on the coordinate scale: $$\widehat{\mathrm{Cov}}(\hat\eta) \approx I^{-1}.$$ Standard errors and correlations are reported on the coordinate scale, and approximate intervals $\hat\eta \pm z\,\mathrm{se}$ are mapped back to the natural scale, so that intervals for positive parameters stay positive. Estimated mean coefficients are profiled out (the likelihood is the profile, restricted, or marginal likelihood of the mean), and the approximation conditions on them.

The approximation needs a maximum: gradient reports the gradient at the estimate, which should be close to zero. Some combinations of hyperparameters can be nearly unidentifiable: under infill asymptotics, only $\sigma^2 / \ell^{2\nu}$ is consistently estimable for a Matérn kernel of smoothness $\nu$ (Zhang, 2004). The likelihood is then nearly flat along a ridge: the information is nearly singular in absolute terms, though not always in relative ones, because the curvature across the ridge can be small too. poorly_identified is therefore TRUE when the information is not positive definite, when its condition number exceeds condition_threshold, or when the standard deviation along the least determined combination of the log-scale parameters, weakest_log_scale_sd, exceeds log_scale_threshold. weakest_direction is the eigenvector of the smallest eigenvalue, the direction along the ridge; for a Matérn kernel under infill it is close to $(2\nu, 1) / \sqrt{4\nu^2 + 1}$ in (log variance, log length scale). Real-valued parameters, such as changepoint locations, have no common scale with the others and enter only the first two criteria. gp_profile_likelihood() shows the likelihood along one parameter without the quadratic approximation.

Stability

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

References

Zhang, H. (2004). Inconsistent estimation and asymptotically equal interpolations in model-based geostatistics. Journal of the American Statistical Association, 99(465), 250–261.

Examples

set.seed(1)
x <- seq(0, 10, length.out = 60)
y <- sin(x) + rnorm(60, sd = 0.2)
model <- optimize_gp(x, y, rbf_kernel(), noise_variance = 0.1, n_starts = 1)

uncertainty <- gp_hyperparameter_uncertainty(model)
uncertainty
#> Hyperparameter uncertainty (Laplace approximation on the optimizer scale)
#>       parameter estimate coordinate standard_error lower 95% upper 95%
#>        variance 1.892900        log          1.010   0.26020  13.77000
#>    length_scale 2.135400        log          0.229   1.36200   3.34800
#>  noise_variance 0.031965        log          0.196   0.02176   0.04695
#>   condition number of the information: 78.1
#>   weakest log-scale combination: standard deviation 1.03
#>   poorly identified along: 0.98 variance + 0.19 length_scale (coordinate scale)
uncertainty$correlation
#>                  variance length_scale noise_variance
#> variance        1.0000000    0.8545129     -0.0938447
#> length_scale    0.8545129    1.0000000     -0.1067184
#> noise_variance -0.0938447   -0.1067184      1.0000000