Fit a GP to ordered time-series observations

Usage

fit_time_series_gp(
  time,
  y,
  kernel,
  unit = "auto",
  method = c("exact", "state_space"),
  ...
)

Arguments

time: Numeric, Date, or POSIXt vector in strictly increasing order.
y: Numeric response vector.
kernel: Gaussian-process kernel specification.
unit: Time unit passed to gp_time_index().
method: "exact" conditions a Gaussian prior on the observations with fit_gp(). "state_space" computes the same posterior and log marginal likelihood in \(O(n)\) time by Kalman filtering and smoothing; it needs a Matérn-1/2, 3/2, or 5/2 kernel, a sum of such kernels, or scaled versions of them (see State-space inference).
...: Additional arguments passed to fit_gp(). With method = "state_space" only noise_variance (positive) and mean (with fixed coefficients) apply.

Value

A fitted model containing explicit time-index metadata used by forecast_gp(): a gaussianprocesses_model for method = "exact" and a gaussianprocesses_state_space_model for method = "state_space".

Details

This function was called fit_gp_time() in version 0.1.0.

State-space inference

On one input, a Matérn process with \(\nu = p + 1/2\) is the first component of the state \(x = (f, f', \ldots, f^{(p)})\) of the linear stochastic differential equation \(dx = F x\,dt + L\,dW\), with \(F\) the companion matrix of \((s + \lambda)^{p+1}\), \(\lambda = \sqrt{2\nu}/\ell\), and white noise \(W\) of spectral density \(q_c\) (Hartikainen and Särkkä, 2010). The stationary state covariance \(P_\infty\) solves \(F P_\infty + P_\infty F^\top + L q_c L^\top = 0\), and \(k(\tau) = H e^{F\tau} P_\infty H^\top\) for \(\tau \ge 0\) with \(H = (1, 0, \ldots, 0)\). Sums of kernels stack their states, and scaled kernels scale \(P_\infty\) and \(q_c\). Between consecutive times the state moves by \(A = e^{F\Delta t}\), which has a closed form because \(F + \lambda I\) is nilpotent, with process noise \(P_\infty - A P_\infty A^\top\).

The Kalman filter (with Joseph-form updates) gives the log marginal likelihood by the prediction-error decomposition, the Rauch-Tung-Striebel smoother gives the posterior at the training times, and forecast_gp() interleaves the requested times as unobserved steps. The results equal those of the exact GP up to rounding, at \(O(n)\) cost instead of \(O(n^3)\), and no numerical jitter is used. sample_gp_posterior() draws by forward filtering, backward sampling. Kernels without an exact finite-dimensional state (RBF, rational quadratic, periodic, linear, and products) raise an error of class gaussianprocesses_state_space_error.

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. method = "state_space" is experimental.

References

Hartikainen, J., and Särkkä, S. (2010). Kalman filtering and smoothing solutions to temporal Gaussian process regression models. IEEE International Workshop on Machine Learning for Signal Processing, 379–384.

Särkkä, S., and Solin, A. (2019). Applied Stochastic Differential Equations. Cambridge University Press.

Examples

time <- as.Date("2026-01-01") + 0:23
y <- sin(2 * pi * (0:23) / 12)

model <- fit_time_series_gp(
  time,
  y,
  kernel = time_series_kernel(period = 12),
  noise_variance = 0.05
)
model$time[c("origin", "unit")]
#> $origin
#> [1] "2026-01-01"
#> 
#> $unit
#> [1] "days"
#> 

# The same posterior by Kalman smoothing, for a Matern kernel.
kernel <- matern32_kernel(length_scale = 3)
exact <- fit_time_series_gp(time, y, kernel, noise_variance = 0.05)
kalman <- fit_time_series_gp(time, y, kernel, noise_variance = 0.05,
  method = "state_space")
c(logLik(exact), logLik(kalman))
#> [1] -10.59189 -10.59189