gaussianprocesses

• predict_heteroscedastic_gp() no longer returns noise_variance; use observation_noise_variance.
• CRAN preparation (#40).
• citation("gaussianprocesses") gives a citation from inst/CITATION.
• DESCRIPTION has no BugReports field: GitLab serves the issue list only to logged-in users, and R's incoming check flags the public work items page. The README links the tracker.
• CI runs R 4.6.1, adds a manual cran-check job that installs every suggested package and checks the built tarball with its tests, and a manual text-checks job for spelling and URLs. RELEASING.md describes CRAN submission and resubmission.
• Stability tiers and the versioning policy for 1.0.
• Every exported function is stable, experimental, or deprecated. Its help page says which in a new "Stability" section, and the reference index is grouped by tier. ?gaussianprocesses sets out the policy. From 1.0.0, stable interfaces change incompatibly only in a major release, after a deprecation period, and experimental interfaces may change in a minor release, with a NEWS entry.
• Experimental:
• the sparse approximations;
• latent models with their likelihoods and scores;
• the heteroscedastic method;
• state-space inference (optimize_time_series_gp() and fit_time_series_gp(method = "state_space"));
• the spectral-mixture and changepoint kernels;
• derivative observations (fit_gp(derivative = )) and predict_gradient_gp();
• hyperparameter uncertainty;
• the numerical and calibration diagnostics;
• the benchmark and simulation-scenario helpers.

Everything else is stable.

• Schema 4 is the model format of version 1.0, and every 1.x release will read models saved by any earlier version. The tests now read models saved by the 0.1.0 release and a set of schema-4 models, and check that they give the results they were saved with.
• CONTRIBUTING.md records the versioning and deprecation policy and two decisions for 1.0: the kernel registry stays internal, and the package stays pure R, with the measurements behind both.
• Validation against independent references.
• inst/validation/references.R defines 35 references, at least one for every inference method: other implementations (DiceKriging, kernlab, nlme, gplite, and LAPACK's pivoted Cholesky decomposition), brute-force computations (dense linear algebra, refits, numerical integration and differentiation, optim()), and closed forms. Each has a tolerance and, for other packages, the version it was set with. A test runs every reference, skipping those whose package is not installed, and checks that every exported function belongs to a validated method or is listed as performing no inference.
• The DiceKriging and kernlab comparison scripts in inst/benchmarks and the gplite, nlme, and DiceKriging comparison tests are replaced by references. CI now installs DiceKriging and kernlab, so every reference runs there.
• A new validation article tabulates the references as it is built and reproduces three analyses on data shipped with R: Rasmussen and Williams's Mauna Loa CO2 kernel on datasets::co2, a classifier for versicolor and virginica irises against logistic regression, and a Poisson model of datasets::discoveries against a quadratic Poisson regression.
• State-space (Kalman) inference for time series.
• fit_time_series_gp(method = "state_space") computes the exact GP posterior and log marginal likelihood in O(n) time by Kalman filtering and Rauch-Tung-Striebel smoothing, for Matérn-1/2, 3/2, and 5/2 kernels, their sums, and scaled versions of them. The kernel registry supplies their stochastic differential equations, and the transition matrices are in closed form. Results agree with exact inference to a relative 1e-8 on irregular and tied inputs. Other kernels raise an error of class gaussianprocesses_state_space_error.
• forecast_gp(), predict(), logLik(), fitted(), residuals(), summary(), log_marginal_likelihood(), and rolling_origin_gp(method = "state_space") accept the new models, and sample_gp_posterior() draws by forward filtering, backward sampling.
• New optimize_time_series_gp() estimates the hyperparameters of a time-series GP with either method; the state-space likelihood is optimized with finite-difference gradients.
• inst/benchmarks/state-space-scaling.R times the method from 1000 to 100 000 observations, and the time-series vignette compares state-space and exact forecasts.
• Inducing-point selection and optimization for sparse Gaussian processes.
• select_inducing_points(method = "variance", kernel = ...) performs greedy conditional-variance selection, the pivot order of a pivoted Cholesky decomposition of K(X, X), in O(n m^2) time without forming K(X, X). The sparse fitters accept selection = "variance".
• optimize_sparse_gp(optimize_inducing = TRUE) optimizes the inducing inputs jointly with the hyperparameters, with analytical gradients from the kernel input derivatives. It is recommended for VFE; with FITC it overfits.
• Sparse fits detect inducing points that are nearly identical under the kernel, record them in $inducing_collapse, and warn (class gaussianprocesses_inducing_collapse_warning).
• The research note on sparse approximations reports the gap to the exact marginal likelihood that optimizing the inducing points closes on the simulation scenarios, reproduces FITC's noise underestimation with free inducing points, and times the greedy selection.
• Variational sparse Gaussian processes, gradients, and sparse sampling.
• fit_sparse_gp() gains method = "vfe", the collapsed variational bound of Titsias (2009), which lies below the exact log marginal likelihood and equals it when every input is an inducing point. VFE predicts with the deterministic training conditional. FITC stays the default.
• Sparse models report the trace term tr(K_ff - Q_ff) / sigma^2, a diagnostic of how well the inducing points cover the data, and record their method; models now use schema 4, and older sparse models are read as FITC.
• New optimize_sparse_gp() optimizes kernel parameters and noise for VFE (the default) or FITC with fixed inducing points, using the optimizer of optimize_gp() and analytical gradients computed in O(n m^2) time without n x n matrices.
• sample_gp_posterior() draws from FITC and VFE posteriors, and predict_sparse_gp() predicts marginal variances in blocks (block_size), with memory linear in the number of prediction inputs.
• A research note compares VFE and FITC on the failure modes reported by Bauer, van der Wilk, and Rasmussen (2016), and reports time and memory for up to 50 000 observations and 500 inducing points.
• Poisson count models with latent Gaussian processes.
• Latent models accept mean coefficients estimated by maximum likelihood (estimate_coefficients("ml"), the default of linear_mean()): fit_latent_gp() maximizes the approximate marginal likelihood over them with its analytical gradient, and optimize_latent_gp() estimates them jointly with the hyperparameters as mean.<name>. A constant mean estimated this way is the log level of a count model. coef() returns the coefficients of latent models.
• predict_latent_gp() returns, for Poisson likelihoods, rate_interval, the exact interval of the lognormal rate, and prediction_interval, equal-tailed intervals for new counts from the Poisson-lognormal predictive distribution, which are conservative because counts are discrete; for Gaussian likelihoods it returns the prediction_interval of the exact model. sample_gp_posterior() returns rate draws.
• The Poisson-lognormal distribution function is computed by quadrature in the latent value or, where the Poisson probabilities change too sharply for that, through the gamma representation of the Poisson distribution function; it agrees with numerical integration to 1e-9 for rates from 0.3 to 5000.
• New gp_count_scores() scores count forecasts with the log score, Dawid-Sebastiani score, Pearson dispersion (an overdispersion check), coverage of the count intervals, and randomized PIT values; gp_holdout_scores() returns them for latent Poisson models.
• A new article analyses datasets::discoveries and reports a simulation: rate intervals cover the true rate at their nominal level with known hyperparameters, count intervals are conservative, randomized PIT values are uniform, and the dispersion check detects overdispersion.
• The Poisson fits agree with gplite.
• Binary classification with latent Gaussian processes.
• New gp_classification_scores() scores probability forecasts with the log loss, Brier score, ROC AUC, reliability bins, and expected calibration error. gp_holdout_scores() returns these scores for latent models with a Bernoulli likelihood, with the log loss from predictive log probabilities, and the scores of the Gaussian predictive for a Gaussian likelihood.
• predict_latent_gp() returns probability_interval, the interval of the class probability itself, and sample_gp_posterior() returns probability draws for Bernoulli likelihoods.
• optimize_latent_gp() warns (class gaussianprocesses_bound_warning) when a parameter ends at a bound of the search and records it in $optimization$at_bound. With separable classes the signal variance grows to a bound or to a very large value.
• A new research note simulates 300 data sets from a known GP classifier. With the true hyperparameters, the Laplace predictive probabilities are slightly underconfident (pooled reliability deviations of up to 7 standard errors), while their Brier score exceeds the oracle's by only 0.015; the note explains why and documents the plug-in and MacKay shortcuts.
• The Laplace fits agree with the gplite package (now suggested) to 1e-6 in the mode, approximate marginal likelihood, and latent predictions.
• Latent Gaussian-process models with the Laplace approximation, for responses observed through a likelihood such as binary labels or counts.
• New fit_latent_gp(x, y, kernel, likelihood, ...) finds the posterior mode by Newton's method on the well-conditioned matrix B = I + W^(1/2) K W^(1/2) (Rasmussen and Williams, Algorithm 3.1), with the package's jitter policy, a step-halving line search, and a convergence test on the Newton decrement. $laplace reports the iterations, step halvings, decrement, and gradient norm, and a warning is raised if the iterations do not converge. Likelihoods that are not log-concave are refused. The Newton update is written so that it stays accurate for counts in the thousands, where the usual form loses all digits.
• log_marginal_likelihood() and log_marginal_likelihood_gradient() return the Laplace approximation and its analytical gradient, including the implicit dependence of the mode on the hyperparameters, for kernel and likelihood parameters.
• New optimize_latent_gp() optimizes the hyperparameters with the optimizer of optimize_gp(), which now shares its setup, multi-start, and diagnostics code, and warm-starts each Newton iteration from the previous mode.
• New predict_latent_gp() gives the latent posterior (Algorithm 3.2) and the predictive distribution of new responses through the likelihood, including class probabilities. sample_gp_posterior() draws latent values from the Gaussian approximation and observations from the likelihood. predict(), logLik(), nobs(), fitted(), and summary() have methods for latent models.
• Means may have fixed coefficients or a Gaussian coefficient prior.
• With gaussian_likelihood() the approximation is exact, and the model equals fit_gp() to 1e-10. ?fit_latent_gp reports the error of the approximation against tensor-product quadrature for up to three observations.
• Likelihoods for non-Gaussian responses, the foundation for latent Gaussian-process models.
• New gaussian_likelihood(), bernoulli_likelihood() (logit and probit links), and poisson_likelihood() (log link, with exposure) build likelihood specifications: plain data with a print method, like kernel specifications. Binary responses may be 0/1, logical, or a two-level factor whose second level is the event, as in glm().
• New evaluate_likelihood() returns the log density and its first three derivatives in the latent values. They stay finite and accurate for latent values up to 40 in magnitude: the probit derivatives use the inverse Mills ratio in log space, and its continued fraction in the far tail.
• New likelihood_predictive() gives the predictive mean, variance, class probability, and log density of new responses when the latent value is Gaussian, in closed form where one exists and otherwise by adaptive Gauss-Hermite quadrature centred at the mode of the integrand.
• Gauss-Hermite rules come from the Golub-Welsch eigenvalue method, with weights from Christoffel sums that keep the smallest weights accurate.
• Non-Gaussian models will be fitted through one entry point that takes a likelihood; fit_gp() remains the exact Gaussian model (see the API audit).
• Proper scoring rules, hyperparameter uncertainty, and multi-start diagnostics.
• New gp_scores() scores Gaussian predictions with the log score, the CRPS in closed form, interval scores and coverage across a grid of levels, PIT values, SMSE, and MSLL. gp_holdout_scores() scores any model on held-out data, and gp_loo_scores() on its exact leave-one-out predictions.
• New gp_hyperparameter_uncertainty() computes the observed information at the estimate from differences of the analytical gradient on the optimizer scale, with standard errors, correlations, intervals, and a flag for poorly identified hyperparameters. It detects the variance-length-scale ridge of a Matérn process observed within one length scale, and stays silent when the data span twenty length scales. New gp_profile_likelihood() profiles one hyperparameter.
• optimize_gp() records the distinct optima its starts reached in $optimization$optima, and summary() reports them.
• A new research note checks the standard errors in a 200-replicate simulation: for the noise variance they match the spread of the estimates, and for the variance and length scale they are 13% and 20% too small at 60 observations, with 95% intervals covering about 90% of the time.
• New spectral_mixture_kernel() (Wilson and Adams, 2013): a sum of components w exp(-2 pi^2 tau^2 v) cos(2 pi tau mu), one weight, frequency, and spectral variance per component and input dimension. A component with frequency 0 is an RBF kernel.
• Each component is a kernel registry entry, so hyperparameter gradients, input derivatives, and derivative observations work for it. Frequencies are real-valued (identity coordinate).
• New