Learn R Programming

pomp (version 1.19)

Power spectrum computation and matching: Power spectrum computation and spectrum-matching for partially-observed Markov processes

Description

spect estimates the power spectrum of time series data and model simulations and compares the results. It can be used to diagnose goodness of fit and/or as the basis for frequency-domain parameter estimation (spect.match).

spect.match tries to match the power spectrum of the model to that of the data. It calls an optimizer to adjust model parameters to minimize the discrepancy between simulated and actual data.

Usage

# S4 method for pomp
spect(object, params, vars, kernel.width, nsim, seed = NULL,
     transform.data,detrend = c("none","mean","linear","quadratic"),
     …)
# S4 method for spect.pomp
spect(object, params, vars, kernel.width, nsim, seed = NULL,
     transform.data, detrend, …)
# S4 method for pomp
spect.match(object, start, est = character(0), vars, nsim,
     seed = NULL, kernel.width, transform.data,
     detrend = c("none","mean","linear","quadratic"), weights = 1,
     method = c("subplex","Nelder-Mead","SANN"), verbose = getOption("verbose"),
     fail.value = NA, …)
# S4 method for spect.pomp
spect.match(object, start, est = character(0), vars, nsim,
     seed = NULL, kernel.width, transform.data, detrend, weights = 1,
     method = c("subplex","Nelder-Mead","SANN"), verbose = getOption("verbose"),
     fail.value = NA, …)
# S4 method for pomp
spect.match(object, start, est, vars, nsim,
     seed = NULL, kernel.width, transform.data, detrend,
     weights, method, verbose = getOption("verbose"),
     fail.value, …)

Arguments

object

An object of class pomp.

params

optional named numeric vector of model parameters. By default, params=coef(object).

vars

optional; names of observed variables for which the power spectrum will be computed. This must be a subset of rownames(obs(object)). By default, the spectrum will be computed for all observables.

kernel.width

width parameter for the smoothing kernel used for calculating the estimate of the spectrum.

nsim

number of model simulations to be computed.

seed

optional; if non-NULL, the random number generator will be initialized with this seed for simulations. See simulate.

transform.data

function; this transformation will be applied to the observables prior to estimation of the spectrum, and prior to any detrending.

detrend

de-trending operation to perform. Options include no detrending, and subtraction of constant, linear, and quadratic trends from the data. Detrending is applied to each data series and to each model simulation independently.

weights

optional. The mismatch between model and data is measured by a weighted average of mismatch at each frequency. By default, all frequencies are weighted equally. weights can be specified either as a vector (which must have length equal to the number of frequencies) or as a function of frequency. If the latter, weights(freq) must return a nonnegative weight for each frequency.

start

named numeric vector; the initial guess of parameters.

est

character vector; the names of parameters to be estimated.

method

Optimization method. Choices are subplex, sannbox, and any of the methods used by optim.

verbose

logical; print diagnostic messages?

fail.value

optional scalar; if non-NA, this value is substituted for non-finite values of the objective function.

Additional arguments. In the case of spect, these are currently ignored. In the case of spect.match, these are passed to optim or subplex in the control list.

Value

spect returns an object of class spect.pomp, which is derived from class pomp and therefore has all the slots of that class. In addition, spect.pomp objects have the following slots:

kernel.width

width parameter of the smoothing kernel used.

transform

transformation function used.

freq

numeric vector of the frequencies at which the power spectrum is estimated.

datspec, simspec

estimated power spectra for data and simulations, respectively.

pvals

one-sided p-values: fraction of the simulated spectra that differ more from the mean simulated spectrum than does the data. The metric used is \(L^2\) distance.

detrend

detrending option used.

spect.match returns an object of class spect.matched.pomp, which is derived from class spect.pomp and therefore has all the slots of that class. In addition, spect.matched.pomp objects have the following slots:

est, weights, fail.value

values of the corresponding arguments in the call to spect.match.

evals

number of function and gradient evaluations by the optimizer. See optim.

value

Value of the objective function.

convergence, msg

Convergence code and message from the optimizer. See optim.

Details

A call to spect results in the estimation of the power spectrum for the (transformed, detrended) data and nsim model simulations. The results of these computations are stored in an object of class spect.pomp.

A call to spect.match results in an attempt to optimize the agreement between model and data spectrum over the parameters named in est. The results, including coefficients of the fitted model and power spectra of fitted model and data, are stored in an object of class spect.matched.pomp.

References

D.C. Reuman, R.A. Desharnais, R.F. Costantino, O. Ahmad, J.E. Cohen (2006) Power spectra reveal the influence of stochasticity on nonlinear population dynamics. Proceedings of the National Academy of Sciences 103, 18860-18865.

D.C. Reuman, R.F. Costantino, R.A. Desharnais, J.E. Cohen (2008) Color of environmental noise affects the nonlinear dynamics of cycling, stage-structured populations. Ecology Letters, 11, 820-830.

See Also

pomp, probe

Examples

Run this code
# NOT RUN {
pompExample(ou2)
good <- spect(
              ou2,
              vars=c("y1","y2"),
              kernel.width=3,
              detrend="mean",
              nsim=500
              )
summary(good)
plot(good)

ou2.bad <- ou2
coef(ou2.bad,c("x1.0","x2.0","alpha.1","alpha.4")) <- c(0,0,0.1,0.2)
bad <- spect(
             ou2.bad,
             vars=c("y1","y2"),
             kernel.width=3,
             detrend="mean",
             nsim=500
             )
summary(bad)
plot(bad)
# }

Run the code above in your browser using DataLab