Movatterモバイル変換

Superfast inference for stationary Gaussian time series.

Description

Likelihood evaluations for stationary Gaussian time series are typically obtained via the Durbin-Levinson algorithm, which scales as O(n^2) in the number of time series observations. This package provides a "superfast" O(n log^2 n) algorithm written in C++, crossing over with Durbin-Levinson around n = 300. Efficient implementations of the score and Hessian functions are also provided, leading to superfast versions of inference algorithms such as Newton-Raphson and Hamiltonian Monte Carlo. The C++ code provides a Toeplitz matrix class packaged as a header-only library, to simplify low-level usage in other packages and outside of R.

Details

While likelihood calculations with stationary Gaussian time series generally scale asO(N^2) in the number of observations, this package implements an algorithm which scales as⁠O(N log^2 N)⁠. "Superfast" algorithms for loglikelihood gradients and Hessians are also provided. The underlying C++ code is distributed through a header-only library found in the installed package'sinclude directory.

Author(s)

Maintainer: Martin Lysymlysy@uwaterloo.ca

Authors:

• Yun Ling

See Also

Useful links:

• https://github.com/mlysy/SuperGauss

• Report bugs athttps://github.com/mlysy/SuperGauss/issues

Examples

# Superfast inference for the timescale parameter# of the exponential autocorrelation functionexp_acf <- function(lambda) exp(-(1:N-1)/lambda)# simulate datalambda0 <- 1N <- 1000X <- rnormtz(n = 1, acf = exp_acf(lambda0))# loglikelihood function# allocate memory for a NormalToeplitz distribution objectNTz <- NormalToeplitz$new(N)loglik <- function(lambda) {  NTz$logdens(z = X, acf = exp_acf(lambda))  ## dSnorm(X = X, acf = Toep, log = TRUE)}# maximum likelihood estimationoptimize(f = loglik, interval = c(.2, 5), maximum = TRUE)

Cholesky multiplication with Toeplitz variance matrices.

Description

Multiplies the Cholesky decomposition of the Toeplitz matrix with another matrix, or solves a system of equations with the Cholesky factor.

Usage

cholZX(Z, acf)cholXZ(X, acf)

Arguments

Details

IfC == t(chol(toeplitz(acf))), thencholZX() computesC %*% Z andcholZX() computessolve(C, X). Both functions use the Durbin-Levinson algorithm.

Value

Size⁠N x p⁠ residual or observation matrix.

Examples

N <- 10p <- 2acf <- exp(-(1:N - 1))Z <- matrix(rnorm(N * p), N, p)cholZX(Z = Z, acf = acf) - (t(chol(toeplitz(acf))) %*%  Z)X <- matrix(rnorm(N * p), N, p)cholXZ(X = X, acf = acf) - solve(t(chol(toeplitz(acf))), X)

Constructor and methods for Circulant matrix objects.

Description

Constructor and methods for Circulant matrix objects.

Methods

Public methods

• Circulant$new()

• Circulant$size()

• Circulant$set_acf()

• Circulant$get_acf()

• Circulant$set_psd()

• Circulant$get_psd()

• Circulant$has_acf()

• Circulant$prod()

• Circulant$solve()

• Circulant$log_det()

• Circulant$clone()

Methodnew()

Class constructor.

Usage

Circulant$new(N, uacf, upsd)

Arguments

Returns

ACirculant object.

Methodsize()

Get the size of the Circulant matrix.

Usage

Circulant$size()

Returns

Size of the Circulant matrix.

Methodset_acf()

Set the autocorrelation of the Circulant matrix.

Usage

Circulant$set_acf(uacf)

Arguments

Methodget_acf()

Get the autocorrelation of the Circulant matrix.

Usage

Circulant$get_acf()

Returns

The complete autocorrelation vector of lengthN.

Methodset_psd()

Set the PSD of the Circulant matrix.

The power spectral density (PSD) of a Circulant matrixCt = Circulant(acf) is defined aspsd = iFFT(acf).

Usage

Circulant$set_psd(upsd)

Arguments

Methodget_psd()

Get the PSD of the Circulant matrix.

Usage

Circulant$get_psd()

Returns

The complete PSD vector of lengthN.

Methodhas_acf()

Check whether the autocorrelation of the Circulant matrix has been set.

Usage

Circulant$has_acf()

Returns

Logical;TRUE ifCirculant$set_acf() has been called.

Methodprod()

Circulant matrix-matrix product.

Usage

Circulant$prod(x)

Arguments

Returns

The matrix productCt %*% x.

Methodsolve()

Solve a Circulant system of equations.

Usage

Circulant$solve(x)

Arguments

Returns

The solution inz to the system of equationsCt %*% z = x. Ifx is missing, returns the inverse ofCt.

Methodlog_det()

Calculate the log-determinant of the Circulant matrix.

Usage

Circulant$log_det()

Returns

The log-determinantlog(det(Ct)).

Methodclone()

The objects of this class are cloneable with this method.

Usage

Circulant$clone(deep = FALSE)

Arguments

Multivariate normal with Circulant variance matrix.

Description

Provides methods for the Normal-Circulant (NCt) distribution, which for a random vectorz of lengthN is defined as

z ~ NCt(uacf)   <=>   z ~ Normal(0, toeplitz(acf)),

whereuacf are theNu = floor(N/2)+1 unique elements of the autocorrelation vectoracf, i.e.,

acf = (uacf, rev(uacf[2:(Nu-1)]),   N even,    = (uacf, rev(uacf[2:Nu])),      N odd.

Methods

Public methods

• NormalCirculant$new()

• NormalCirculant$size()

• NormalCirculant$logdens()

• NormalCirculant$grad_full()

• NormalCirculant$clone()

Methodnew()

Class constructor.

Usage

NormalCirculant$new(N)

Arguments

Returns

ANormalCirculant object.

Methodsize()

Get the size of the NCt random vector.

Usage

NormalCirculant$size()

Returns

Size of the NCt random vector.

Methodlogdens()

Log-density function.

Usage

NormalCirculant$logdens(z, uacf)

Arguments

Returns

A scalar or vector of lengthn_obs containing the log-density of the NCt evaluated at its arguments.

Methodgrad_full()

Full gradient of log-density function.

Usage

NormalCirculant$grad_full(z, uacf, calc_dldz = TRUE, calc_dldu = TRUE)

Arguments

Returns

A list with elements:

ldens

The log-density evaluated atz anduacf.

dldz

The length-N gradient vector with respect toz, ifcalc_dldz = TRUE.

dldu

The length-Nu = floor(N/2)+1 gradient vector with respect touacf, ifcalc_dldu = TRUE.

Methodclone()

The objects of this class are cloneable with this method.

Usage

NormalCirculant$clone(deep = FALSE)

Arguments

Multivariate normal with Toeplitz variance matrix.

Description

Provides methods for the Normal-Toeplitz (NTz) distribution defined as

z ~ NTz(acf)   <=>   z ~ Normal(0, toeplitz(acf)),

i.e., for a multivariate normal with mean zero and varianceTz = toeplitz(acf).

Methods

Public methods

• NormalToeplitz$new()

• NormalToeplitz$size()

• NormalToeplitz$logdens()

• NormalToeplitz$grad()

• NormalToeplitz$hess()

• NormalToeplitz$grad_full()

• NormalToeplitz$clone()

Methodnew()

Class constructor.

Usage

NormalToeplitz$new(N)

Arguments

Returns

ANormalToeplitz object.

Methodsize()

Get the size of the NTz random vector.

Usage

NormalToeplitz$size()

Returns

Size of the NTz random vector.

Methodlogdens()

Log-density function.

Usage

NormalToeplitz$logdens(z, acf)

Arguments

Returns

A scalar or vector of lengthn_obs containing the log-density of the NTz evaluated at its arguments.

Methodgrad()

Gradient of the log-density with respect to parameters.

Usage

NormalToeplitz$grad(z, dz, acf, dacf, full_out = FALSE)

Arguments

Returns

A vector of lengthn_theta containing the gradient of the NTz log-density with respect totheta, or a list with elementsldens andgrad consisting of the log-density and the gradient vector.

Methodhess()

Hessian of log-density with respect to parameters.

Usage

NormalToeplitz$hess(z, dz, d2z, acf, dacf, d2acf, full_out = FALSE)

Arguments

Returns

An⁠n_theta x n_theta⁠ matrix containing the Hessian of the NTz log-density with respect totheta, or a list with elementsldens,grad, andhess consisting of the log-density, its gradient (a vector of sizen_theta), and the Hessian matrix, respectively.

Methodgrad_full()

Full gradient of log-density function.

Usage

NormalToeplitz$grad_full(z, acf, calc_dldz = TRUE, calc_dlda = TRUE)

Arguments

Returns

A list with elements:

ldens

The log-density evaluated atz andacf.

dldz

The length-N gradient vector with respect toz, ifcalc_dldz = TRUE.

dlda

The length-N gradient vector with respect toacf, ifcalc_dlda = TRUE.

Methodclone()

The objects of this class are cloneable with this method.

Usage

NormalToeplitz$clone(deep = FALSE)

Arguments

Defunct functions inSuperGauss.

Description

Defunct functions inSuperGauss.

The following functions have been removed from theSuperGauss package

rSnorm()

Please usernormtz() instead.

dSnorm()

Please usednormtz() instead.

Snorm.grad()

Please use thegrad() method in theNormalToeplitz class.

Snorm.hess()

Please use thehess() method in theNormalToeplitz class.

Constructor and methods for Toeplitz matrix objects.

Description

TheToeplitz class contains efficient methods for linear algebra with symmetric positive definite (i.e., variance) Toeplitz matrices.

Usage

is.Toeplitz(x)as.Toeplitz(x)## S3 method for class 'Toeplitz'dim(x)

Arguments

Details

An⁠N x N⁠ Toeplitz matrixTz is defined by its length-N "autocorrelation" vectoracf, i.e., first row/columnTz. Thus, for the functionstats::toeplitz(), we haveTz = toeplitz(acf).

It is assumed thatacf defines a valid (i.e., positive definite) variance matrix. The matrix multiplication methods still work when this is not the case but the other methods do not (return values typically containNaNs).

as.Toeplitz(x) attempts to convert its argument to aToeplitz object by callingToeplitz$new(acf = x).is.Toeplitz(x) checks whether its argument is aToeplitz object.

Methods

Public methods

• Toeplitz$new()

• Toeplitz$print()

• Toeplitz$size()

• Toeplitz$set_acf()

• Toeplitz$get_acf()

• Toeplitz$has_acf()

• Toeplitz$prod()

• Toeplitz$solve()

• Toeplitz$log_det()

• Toeplitz$trace_grad()

• Toeplitz$trace_hess()

• Toeplitz$clone()

Methodnew()

Class constructor.

Usage

Toeplitz$new(N, acf)

Arguments

Returns

AToeplitz object.

Methodprint()

Print method.

Usage

Toeplitz$print()

Methodsize()

Get the size of the Toeplitz matrix.

Usage

Toeplitz$size()

Returns

Size of the Toeplitz matrix.ncol(),nrow(), anddim() methods forToeplitz objects also work as expected.

Methodset_acf()

Set the autocorrelation of the Toeplitz matrix.

Usage

Toeplitz$set_acf(acf)

Arguments

Methodget_acf()

Get the autocorrelation of the Toeplitz matrix.

Usage

Toeplitz$get_acf()

Returns

The autocorrelation vector of lengthN.

Methodhas_acf()

Check whether the autocorrelation of the Toeplitz matrix has been set.

Usage

Toeplitz$has_acf()

Returns

Logical;TRUE ifToeplitz$set_acf() has been called.

Methodprod()

Toeplitz matrix-matrix product.

Usage

Toeplitz$prod(x)

Arguments

Returns

The matrix productTz %*% x.Tz %*% x andx %*% Tz also work as expected.

Methodsolve()

Solve a Toeplitz system of equations.

Usage

Toeplitz$solve(x, method = c("gschur", "pcg"), tol = 1e-10)

Arguments

Returns

The solution inz to the system of equationsTz %*% z = x. Ifx is missing, returns the inverse ofTz.solve(Tz, x) andsolve(Tz, x, method, tol) also work as expected.

Methodlog_det()

Calculate the log-determinant of the Toeplitz matrix.

Usage

Toeplitz$log_det()

Returns

The log-determinantlog(det(Tz)).determinant(Tz) also works as expected.

Methodtrace_grad()

Computes the trace-gradient with respect to Toeplitz matrices.

Usage

Toeplitz$trace_grad(acf2)

Arguments

Returns

Computes the trace of

solve(Tz, toeplitz(acf2)).

This is used in the computation of the gradient oflog(det(Tz(theta))) with respect totheta.

Methodtrace_hess()

Computes the trace-Hessian with respect to Toeplitz matrices.

Usage

Toeplitz$trace_hess(acf2, acf3)

Arguments

Returns

Computes the trace of

solve(Tz, toeplitz(acf2)) %*% solve(Tz, toeplitz(acf3)).

This is used in the computation of the Hessian oflog(det(Tz(theta))) with respect totheta.

Methodclone()

The objects of this class are cloneable with this method.

Usage

Toeplitz$clone(deep = FALSE)

Arguments

Examples

# construct a Toeplitz matrixacf <- exp(-(1:5))Tz <- Toeplitz$new(acf = acf)# alternatively, can allocate space firstTz <- Toeplitz$new(N = length(acf))Tz$set_acf(acf = acf)# basic methodsTz$get_acf() # extract the acfdim(Tz) # == c(nrow(Tz), ncol(Tz))Tz # print method# linear algebra methodsX <- matrix(rnorm(10), 5, 2)Tz %*% Xt(X) %*% Tzsolve(Tz, X)determinant(Tz) # log-determinant

Convert position autocorrelations to increment autocorrelations.

Description

Convert the autocorrelation of a stationary sequence⁠x = (x_1, ..., x_N)⁠ to that of its increments,⁠dx = (x_2 - x_1, ..., x_N - x_(N-1))⁠.

Usage

acf2incr(acf)

Arguments

Value

LengthN-1 vector of increment autocorrelations.

Examples

acf2incr(acf = exp(-(0:10)))

Convert autocorrelation of stationary increments to mean squared displacement of posititions.

Description

Converts the autocorrelation of a stationary increment sequence⁠dx = (x_1 - x_0, ..., x_N - x_(N-1))⁠ to the mean squared displacement (MSD) of the corresponding positions, i.e.,MSD_i = E[(x_i - x_0)^2].

Usage

acf2msd(acf)

Arguments

Value

Length-N MSD vector of the corresponding positions.

Examples

acf2msd(acf = exp(-(0:10)))

Density of a multivariate normal with Toeplitz variance matrix.

Description

Density of a multivariate normal with Toeplitz variance matrix.

Usage

dnormtz(X, mu, acf, log = FALSE, method = c("gschur", "ltz"))

Arguments

Value

Vector ofn (log-)densities, one for each column ofX.

Examples

# simulate dataN <- 10 # length of each time seriesn <- 3 # number of time seriestheta <- 0.1lambda <- 2mu <- theta^2 * rep(1, N)acf <- exp(-lambda * (1:N - 1))X <- rnormtz(n, acf = acf) + mu# evaluate log-densitydnormtz(X, mu, acf, log = TRUE)

Mean square displacement of fractional Brownian motion.

Description

Mean square displacement of fractional Brownian motion.

Usage

fbm_msd(tseq, H)

Arguments

Details

The mean squared displacement (MSD) of a stochastic processX_t is defined as

MSD(t) = E[(X_t - X_0)^2].

Fractional Brownian motion (fBM) is a continuous Gaussian process with stationary increments, such that its covariance function is entirely defined the MSD, which in this case is⁠MSD(t) = |t|^(2H)⁠.

Value

Length-N vector of mean square displacements.

Examples

fbm_msd(tseq = 1:10, H = 0.4)

Matern autocorrelation function.

Description

Matern autocorrelation function.

Usage

matern_acf(tseq, lambda, nu)

Arguments

Details

The Matern autocorrelation is given by

\mathrm{\scriptsize ACF}(t) = \frac{2^{1-\nu}}{\Gamma(\nu)} \left(\sqrt{2\nu}\frac{t}{\lambda}\right)^\nu K_\nu\left(\sqrt{2\nu} \frac{t}{\lambda}\right),

whereK_\nu(x) is the modified Bessel function of second kind.

Value

An autocorrelation vector of lengthN.

Examples

matern_acf(tseq = 1:10, lambda = 1, nu = 3/2)

Convert mean square displacement of positions to autocorrelation of increments.

Description

Converts the mean squared displacement (MSD) of a stationary increments sequence⁠x = (x_0, x_1, ..., x_N)⁠ positions to the autocorrelation of the corresponding increments⁠dx = (x_1 - x_0, ..., x_N - x_(N-1))⁠.

Usage

msd2acf(msd)

Arguments

Value

Length-N autocorrelation vector.

Examples

# autocorrelation of fBM incrementsmsd2acf(msd = fbm_msd(tseq = 0:10, H = .3))

Power-exponential autocorrelation function.

Description

Power-exponential autocorrelation function.

Usage

pex_acf(tseq, lambda, rho)

Arguments

Details

The power-exponential autocorrelation function is given by:

\mathrm{\scriptsize ACF}(t) = \exp \left\{-(t/\lambda)^\rho\right\}.

Value

An autocorrelation vector of lengthN.

Examples

pex_acf(tseq = 1:10, lambda = 1, rho = 2)

Simulate a stationary Gaussian time series.

Description

Simulate a stationary Gaussian time series.

Usage

rnormtz(n = 1, acf, Z, fft = TRUE, nkeep, tol = 1e-06)

Arguments

Details

The FFT method fails when the embedding circulant matrix is not positive definite. This is typically due to one of two things:

1. Roundoff error can make tiny eigenvalues appear negative. For this purpose, argumenttol can be used to replace all negative eigenvalues bytol * ev_max, whereev_max is the largest eigenvalue.

2. The autocorrelation is decaying too slowly on the given timescale. To mitigate this, argumentnkeep can be used to supply a longeracf than is required, and keep only the firstnkeep time series observations. For consistency,nkeep also applies to Durbin-Levinson method.

Value

Length-nkeep vector or size⁠nkeep x n⁠ matrix with time series as columns.

Examples

N <- 10acf <- exp(-(1:N - 1)/N)rnormtz(n = 3, acf = acf)

Toeplitz matrix multiplication.

Description

Efficient matrix multiplication with Toeplitz matrix and arbitrary matrix or vector.

Usage

toep.mult(acf, X)

Arguments

Value

AnN-row matrix corresponding totoeplitz(acf) %*% X.

Examples

N <- 20d <- 3acf <- exp(-(1:N))X <- matrix(rnorm(N*d), N, d)toep.mult(acf, X)