Initialise a spectral-mixture kernel from the data's spectrum

Builds a spectral_mixture_kernel() whose components sit at the main peaks of the empirical spectrum of one-dimensional data, a starting point for optimize_gp().

Usage

initialize_spectral_mixture(x, y, n_components = 2L)

Arguments

x

Numeric vector of inputs, regularly or irregularly spaced.

y

Numeric vector of responses.

n_components

Number of components \(Q\).

Value

A spectral-mixture kernel specification with n_components components.

Details

The spectrum of the centred responses is the periodogram (stats::spec.pgram(), without taper) when the inputs are regularly spaced, and the Lomb-Scargle periodogram otherwise, computed on a grid four times finer than the frequency resolution \(1 / T\), for an input range \(T\), up to half the reciprocal median spacing. A mixture of \(Q\) Gaussians and a uniform background is fitted to the spectrum, as a density over frequency weighted by the power, by EM, starting from the \(Q\) largest local maxima. The background absorbs the flat spectrum of observation noise, which would otherwise widen the components.

Component \(q\) takes the mean and variance of Gaussian \(q\) as its frequency and spectral variance, and the share of the variance of \(y\) that the Gaussian explains as its weight; the background share is left to the noise. Spectral variances are at least \((1 / (4 T))^2\), the order of a spectral peak's width.

The likelihood of a spectral mixture is strongly multimodal, and this initialisation is only a starting point: optimize with several starts, for example n_starts = 5, which displace each frequency by one resolution bin in turn.

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

spectral_mixture_kernel().

Examples

set.seed(1)
x <- 0:199
y <- 2 * sin(2 * pi * 0.1 * x) + sin(2 * pi * 0.27 * x) +
  rnorm(200, sd = 0.5)

kernel <- initialize_spectral_mixture(x, y, n_components = 2)
kernel
#> SumKernel(
#>   SpectralComponent(weight=1.92925, frequency=0.0999983, spectral_variance=1.5625e-06)
#>   SpectralComponent(weight=0.444209, frequency=0.270004, spectral_variance=1.5625e-06)
#> )