| Title: | Superfast Likelihood 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. |
| Authors: | Yun Ling [aut], Martin Lysy [aut, cre] |
| Maintainer: | Martin Lysy <[email protected]> |
| License: | GPL-3 |
| Version: | 2.0.3 |
| Built: | 2024-10-31 05:36:30 UTC |
| Source: | https://github.com/mlysy/supergauss |
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.
While likelihood calculations with stationary Gaussian time series generally scale as O(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's include directory.
Maintainer: Martin Lysy [email protected]
Authors:
Yun Ling
Useful links:
# Superfast inference for the timescale parameter # of the exponential autocorrelation function exp_acf <- function(lambda) exp(-(1:N-1)/lambda) # simulate data lambda0 <- 1 N <- 1000 X <- rnormtz(n = 1, acf = exp_acf(lambda0)) # loglikelihood function # allocate memory for a NormalToeplitz distribution object NTz <- NormalToeplitz$new(N) loglik <- function(lambda) { NTz$logdens(z = X, acf = exp_acf(lambda)) ## dSnorm(X = X, acf = Toep, log = TRUE) } # maximum likelihood estimation optimize(f = loglik, interval = c(.2, 5), maximum = TRUE)# Superfast inference for the timescale parameter # of the exponential autocorrelation function exp_acf <- function(lambda) exp(-(1:N-1)/lambda) # simulate data lambda0 <- 1 N <- 1000 X <- rnormtz(n = 1, acf = exp_acf(lambda0)) # loglikelihood function # allocate memory for a NormalToeplitz distribution object NTz <- NormalToeplitz$new(N) loglik <- function(lambda) { NTz$logdens(z = X, acf = exp_acf(lambda)) ## dSnorm(X = X, acf = Toep, log = TRUE) } # maximum likelihood estimation optimize(f = loglik, interval = c(.2, 5), maximum = TRUE)
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)).
acf2incr(acf)acf2incr(acf)
acf |
Length- |
Length N-1 vector of increment autocorrelations.
acf2incr(acf = exp(-(0:10)))acf2incr(acf = exp(-(0:10)))
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].
acf2msd(acf)acf2msd(acf)
acf |
Length- |
Length-N MSD vector of the corresponding positions.
acf2msd(acf = exp(-(0:10)))acf2msd(acf = exp(-(0:10)))
Multiplies the Cholesky decomposition of the Toeplitz matrix with another matrix, or solves a system of equations with the Cholesky factor.
cholZX(Z, acf) cholXZ(X, acf)cholZX(Z, acf) cholXZ(X, acf)
Z |
Length- |
acf |
Length- |
X |
Length- |
If C == t(chol(toeplitz(acf))), then cholZX() computes C %*% Z and cholZX() computes solve(C, X). Both functions use the Durbin-Levinson algorithm.
Size N x p residual or observation matrix.
N <- 10 p <- 2 acf <- 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)N <- 10 p <- 2 acf <- 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.
new()
Class constructor.
Circulant$new(N, uacf, upsd)
NSize of Circulant matrix.
uacfOptional vector of Nu = floor(N/2)+1 unique elements of the autocorrelation.
upsdOptional vector of Nu = floor(N/2)+1 unique elements of the PSD.
A Circulant object.
size()
Get the size of the Circulant matrix.
Circulant$size()
Size of the Circulant matrix.
set_acf()
Set the autocorrelation of the Circulant matrix.
Circulant$set_acf(uacf)
uacfVector of Nu = floor(N/2)+1 unique elements of the autocorrelation.
get_acf()
Get the autocorrelation of the Circulant matrix.
Circulant$get_acf()
The complete autocorrelation vector of length N.
set_psd()
Set the PSD of the Circulant matrix.
The power spectral density (PSD) of a Circulant matrix Ct = Circulant(acf) is defined as psd = iFFT(acf).
Circulant$set_psd(upsd)
upsdVector of Nu = floor(N/2)+1 unique elements of the psd.
get_psd()
Get the PSD of the Circulant matrix.
Circulant$get_psd()
The complete PSD vector of length N.
has_acf()
Check whether the autocorrelation of the Circulant matrix has been set.
Circulant$has_acf()
Logical; TRUE if Circulant$set_acf() has been called.
prod()
Circulant matrix-matrix product.
Circulant$prod(x)
xVector or matrix with N rows.
The matrix product Ct %*% x.
solve()
Solve a Circulant system of equations.
Circulant$solve(x)
xOptional vector or matrix with N rows.
The solution in z to the system of equations Ct %*% z = x. If x is missing, returns the inverse of Ct.
log_det()
Calculate the log-determinant of the Circulant matrix.
Circulant$log_det()
The log-determinant log(det(Ct)).
clone()
The objects of this class are cloneable with this method.
Circulant$clone(deep = FALSE)
deepWhether to make a deep clone.
Density of a multivariate normal with Toeplitz variance matrix.
dnormtz(X, mu, acf, log = FALSE, method = c("gschur", "ltz"))dnormtz(X, mu, acf, log = FALSE, method = c("gschur", "ltz"))
X |
Vector of length |
mu |
Vector or matrix of mean values of compatible dimensions with |
acf |
Vector of length |
log |
Logical; whether to return the multivariate normal density on the log scale. |
method |
Which calculation method to use. Choices are: |
Vector of n (log-)densities, one for each column of X.
# simulate data N <- 10 # length of each time series n <- 3 # number of time series theta <- 0.1 lambda <- 2 mu <- theta^2 * rep(1, N) acf <- exp(-lambda * (1:N - 1)) X <- rnormtz(n, acf = acf) + mu # evaluate log-density dnormtz(X, mu, acf, log = TRUE)# simulate data N <- 10 # length of each time series n <- 3 # number of time series theta <- 0.1 lambda <- 2 mu <- theta^2 * rep(1, N) acf <- exp(-lambda * (1:N - 1)) X <- rnormtz(n, acf = acf) + mu # evaluate log-density dnormtz(X, mu, acf, log = TRUE)
Mean square displacement of fractional Brownian motion.
fbm_msd(tseq, H)fbm_msd(tseq, H)
tseq |
Length- |
H |
Hurst parameter (between 0 and 1). |
The mean squared displacement (MSD) of a stochastic process X_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).
Length-N vector of mean square displacements.
fbm_msd(tseq = 1:10, H = 0.4)fbm_msd(tseq = 1:10, H = 0.4)
Matern autocorrelation function.
matern_acf(tseq, lambda, nu)matern_acf(tseq, lambda, nu)
tseq |
Vector of |
lambda |
Timescale parameter. |
nu |
Smoothness parameter. |
The Matern autocorrelation is given by
where is the modified Bessel function of second kind.
An autocorrelation vector of length N.
matern_acf(tseq = 1:10, lambda = 1, nu = 3/2)matern_acf(tseq = 1:10, lambda = 1, nu = 3/2)
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)).
msd2acf(msd)msd2acf(msd)
msd |
Length- |
Length-N autocorrelation vector.
# autocorrelation of fBM increments msd2acf(msd = fbm_msd(tseq = 0:10, H = .3))# autocorrelation of fBM increments msd2acf(msd = fbm_msd(tseq = 0:10, H = .3))
Provides methods for the Normal-Circulant (NCt) distribution, which for a random vector z of length N is defined as
z ~ NCt(uacf) <=> z ~ Normal(0, toeplitz(acf)),
where uacf are the Nu = floor(N/2)+1 unique elements of the autocorrelation vector acf, i.e.,
acf = (uacf, rev(uacf[2:(Nu-1)]), N even,
= (uacf, rev(uacf[2:Nu])), N odd.
new()
Class constructor.
NormalCirculant$new(N)
NSize of the NCt random vector.
A NormalCirculant object.
size()
Get the size of the NCt random vector.
NormalCirculant$size()
Size of the NCt random vector.
logdens()
Log-density function.
NormalCirculant$logdens(z, uacf)
zDensity argument. A vector of length N or an N x n_obs matrix where each column is an N-dimensional observation.
uacfA vector of length Nu = floor(N/2) containing the first half of the autocorrelation (i.e., first row/column) of the Circulant variance matrix.
A scalar or vector of length n_obs containing the log-density of the NCt evaluated at its arguments.
grad_full()
Full gradient of log-density function.
NormalCirculant$grad_full(z, uacf, calc_dldz = TRUE, calc_dldu = TRUE)
zDensity argument. A vector of length N.
uacfA vector of length Nu = floor(N/2) containing the first half of the autocorrelation (i.e., first row/column) of the Circulant variance matrix.
calc_dldzWhether or not to calculate the gradient with respect to z.
calc_dlduWhether or not to calculate the gradient with respect to uacf.
A list with elements:
ldensThe log-density evaluated at z and uacf.
dldzThe length-N gradient vector with respect to z, if calc_dldz = TRUE.
dlduThe length-Nu = floor(N/2)+1 gradient vector with respect to uacf, if calc_dldu = TRUE.
clone()
The objects of this class are cloneable with this method.
NormalCirculant$clone(deep = FALSE)
deepWhether to make a deep clone.
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 variance Tz = toeplitz(acf).
new()
Class constructor.
NormalToeplitz$new(N)
NSize of the NTz random vector.
A NormalToeplitz object.
size()
Get the size of the NTz random vector.
NormalToeplitz$size()
Size of the NTz random vector.
logdens()
Log-density function.
NormalToeplitz$logdens(z, acf)
zDensity argument. A vector of length N or an N x n_obs matrix where each column is an N-dimensional observation.
acfA vector of length N containing the autocorrelation (i.e., first row/column) of the Toeplitz variance matrix.
A scalar or vector of length n_obs containing the log-density of the NTz evaluated at its arguments.
grad()
Gradient of the log-density with respect to parameters.
NormalToeplitz$grad(z, dz, acf, dacf, full_out = FALSE)
zDensity argument. A vector of length N.
dzAn N x n_theta matrix containing the gradient dz/dtheta.
acfA vector of length N containing the autocorrelation of the Toeplitz variance matrix.
dacfAn N x n_theta matrix containing the gradient dacf/dtheta.
full_outIf TRUE, returns the log-density as well (see 'Value').
A vector of length n_theta containing the gradient of the NTz log-density with respect to theta, or a list with elements ldens and grad consisting of the log-density and the gradient vector.
hess()
Hessian of log-density with respect to parameters.
NormalToeplitz$hess(z, dz, d2z, acf, dacf, d2acf, full_out = FALSE)
zDensity argument. A vector of length N.
dzAn N x n_theta matrix containing the gradient dz/dtheta.
d2zAn N x n_theta x n_theta array containing the Hessian d^2z/dtheta^2.
acfA vector of length N containing the autocorrelation of the Toeplitz variance matrix.
dacfAn N x n_theta matrix containing the gradient dacf/dtheta.
d2acfAn N x n_theta x n_theta array containing the Hessian dacf^2/dtheta^2.
full_outIf TRUE, returns the log-density and its gradient as well (see 'Value').
An n_theta x n_theta matrix containing the Hessian of the NTz log-density with respect to theta, or a list with elements ldens, grad, and hess consisting of the log-density, its gradient (a vector of size n_theta), and the Hessian matrix, respectively.
grad_full()
Full gradient of log-density function.
NormalToeplitz$grad_full(z, acf, calc_dldz = TRUE, calc_dlda = TRUE)
zDensity argument. A vector of length N.
acfA vector of length N containing the autocorrelation of the Toeplitz variance matrix.
calc_dldzWhether or not to calculate the gradient with respect to z.
calc_dldaWhether or not to calculate the gradient with respect to acf.
A list with elements:
ldensThe log-density evaluated at z and acf.
dldzThe length-N gradient vector with respect to z, if calc_dldz = TRUE.
dldaThe length-N gradient vector with respect to acf, if calc_dlda = TRUE.
clone()
The objects of this class are cloneable with this method.
NormalToeplitz$clone(deep = FALSE)
deepWhether to make a deep clone.
Power-exponential autocorrelation function.
pex_acf(tseq, lambda, rho)pex_acf(tseq, lambda, rho)
tseq |
Vector of |
lambda |
Timescale parameter. |
rho |
Power parameter. |
The power-exponential autocorrelation function is given by:
An autocorrelation vector of length N.
pex_acf(tseq = 1:10, lambda = 1, rho = 2)pex_acf(tseq = 1:10, lambda = 1, rho = 2)
Simulate a stationary Gaussian time series.
rnormtz(n = 1, acf, Z, fft = TRUE, nkeep, tol = 1e-06)rnormtz(n = 1, acf, Z, fft = TRUE, nkeep, tol = 1e-06)
n |
Number of time series to generate. |
acf |
Length- |
Z |
Optional size |
fft |
Logical; whether or not to use the |
nkeep |
Length of time series. Defaults to |
tol |
Relative tolerance on negative eigenvalues. See Details. |
The FFT method fails when the embedding circulant matrix is not positive definite. This is typically due to one of two things:
Roundoff error can make tiny eigenvalues appear negative. For this purpose, argument tol can be used to replace all negative eigenvalues by tol * ev_max, where ev_max is the largest eigenvalue.
The autocorrelation is decaying too slowly on the given timescale. To mitigate this, argument nkeep can be used to supply a longer acf than is required, and keep only the first nkeep time series observations. For consistency, nkeep also applies to Durbin-Levinson method.
Length-nkeep vector or size nkeep x n matrix with time series as columns.
N <- 10 acf <- exp(-(1:N - 1)/N) rnormtz(n = 3, acf = acf)N <- 10 acf <- exp(-(1:N - 1)/N) rnormtz(n = 3, acf = acf)
Defunct functions in SuperGauss.
rSnorm()Please use rnormtz() instead.
dSnorm()Please use dnormtz() instead.
Snorm.grad()Please use the grad() method in the NormalToeplitz class.
Snorm.hess()Please use the hess() method in the NormalToeplitz class.
Efficient matrix multiplication with Toeplitz matrix and arbitrary matrix or vector.
toep.mult(acf, X)toep.mult(acf, X)
acf |
Length- |
X |
Vector or matrix of compatible dimensions with |
An N-row matrix corresponding to toeplitz(acf) %*% X.
N <- 20 d <- 3 acf <- exp(-(1:N)) X <- matrix(rnorm(N*d), N, d) toep.mult(acf, X)N <- 20 d <- 3 acf <- exp(-(1:N)) X <- matrix(rnorm(N*d), N, d) toep.mult(acf, X)
The Toeplitz class contains efficient methods for linear algebra with symmetric positive definite (i.e., variance) Toeplitz matrices.
is.Toeplitz(x) as.Toeplitz(x) ## S3 method for class 'Toeplitz' dim(x)is.Toeplitz(x) as.Toeplitz(x) ## S3 method for class 'Toeplitz' dim(x)
x |
An R object. |
An N x N Toeplitz matrix Tz is defined by its length-N "autocorrelation" vector acf, i.e., first row/column Tz. Thus, for the function stats::toeplitz(), we have Tz = toeplitz(acf).
It is assumed that acf 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 contain NaNs).
as.Toeplitz(x) attempts to convert its argument to a Toeplitz object by calling Toeplitz$new(acf = x). is.Toeplitz(x) checks whether its argument is a Toeplitz object.
new()
Class constructor.
Toeplitz$new(N, acf)
NSize of Toeplitz matrix.
acfAutocorrelation vector of length N.
A Toeplitz object.
print()
Print method.
Toeplitz$print()
size()
Get the size of the Toeplitz matrix.
Toeplitz$size()
Size of the Toeplitz matrix. ncol(), nrow(), and dim() methods for Toeplitz objects also work as expected.
set_acf()
Set the autocorrelation of the Toeplitz matrix.
Toeplitz$set_acf(acf)
acfAutocorrelation vector of length N.
get_acf()
Get the autocorrelation of the Toeplitz matrix.
Toeplitz$get_acf()
The autocorrelation vector of length N.
has_acf()
Check whether the autocorrelation of the Toeplitz matrix has been set.
Toeplitz$has_acf()
Logical; TRUE if Toeplitz$set_acf() has been called.
prod()
Toeplitz matrix-matrix product.
Toeplitz$prod(x)
xVector or matrix with N rows.
The matrix product Tz %*% x. Tz %*% x and x %*% Tz also work as expected.
solve()
Solve a Toeplitz system of equations.
Toeplitz$solve(x, method = c("gschur", "pcg"), tol = 1e-10)xOptional vector or matrix with N rows.
methodSolve method to use. Choices are: gschur for a modified version of the Generalized Schur algorithm of Ammar & Gragg (1988), or pcg for the preconditioned conjugate gradient method of Chen et al (2006). The former is faster and obtains the log-determinant as a direct biproduct. The latter is more numerically stable for long-memory autocorrelations.
tolTolerance level for the pcg method.
The solution in z to the system of equations Tz %*% z = x. If x is missing, returns the inverse of Tz. solve(Tz, x) and solve(Tz, x, method, tol) also work as expected.
log_det()
Calculate the log-determinant of the Toeplitz matrix.
Toeplitz$log_det()
The log-determinant log(det(Tz)). determinant(Tz) also works as expected.
trace_grad()
Computes the trace-gradient with respect to Toeplitz matrices.
Toeplitz$trace_grad(acf2)
acf2Length-N autocorrelation vector of the second Toeplitz matrix. This matrix must be symmetric but not necessarily positive definite.
Computes the trace of
solve(Tz, toeplitz(acf2)).
This is used in the computation of the gradient of log(det(Tz(theta))) with respect to theta.
trace_hess()
Computes the trace-Hessian with respect to Toeplitz matrices.
Toeplitz$trace_hess(acf2, acf3)
acf2Length-N autocorrelation vector of the second Toeplitz matrix. This matrix must be symmetric but not necessarily positive definite.
acf3Length-N autocorrelation vector of the third Toeplitz matrix. This matrix must be symmetric but not necessarily positive definite.
Computes the trace of
solve(Tz, toeplitz(acf2)) %*% solve(Tz, toeplitz(acf3)).
This is used in the computation of the Hessian of log(det(Tz(theta))) with respect to theta.
clone()
The objects of this class are cloneable with this method.
Toeplitz$clone(deep = FALSE)
deepWhether to make a deep clone.
# construct a Toeplitz matrix acf <- exp(-(1:5)) Tz <- Toeplitz$new(acf = acf) # alternatively, can allocate space first Tz <- Toeplitz$new(N = length(acf)) Tz$set_acf(acf = acf) # basic methods Tz$get_acf() # extract the acf dim(Tz) # == c(nrow(Tz), ncol(Tz)) Tz # print method # linear algebra methods X <- matrix(rnorm(10), 5, 2) Tz %*% X t(X) %*% Tz solve(Tz, X) determinant(Tz) # log-determinant# construct a Toeplitz matrix acf <- exp(-(1:5)) Tz <- Toeplitz$new(acf = acf) # alternatively, can allocate space first Tz <- Toeplitz$new(N = length(acf)) Tz$set_acf(acf = acf) # basic methods Tz$get_acf() # extract the acf dim(Tz) # == c(nrow(Tz), ncol(Tz)) Tz # print method # linear algebra methods X <- matrix(rnorm(10), 5, 2) Tz %*% X t(X) %*% Tz solve(Tz, X) determinant(Tz) # log-determinant