##' Power spectrum
##'
##' Power spectrum computation and spectrum-matching for partially-observed
##' Markov processes.
##'
##' \code{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
##' (\code{spect.match}).
##'
##' A call to \code{spect} results in the estimation of the power spectrum for
##' the (transformed, detrended) data and \code{nsim} model simulations. The
##' results of these computations are stored in an object of class
##' \sQuote{spectd_pomp}.
##'
##' @name spect
##' @docType methods
##' @rdname spect
##' @aliases spect spect,missing-method spect,ANY-method
##' @family summary statistics methods
##' @family elementary POMP methods
##'
##' @include simulate.R pomp.R
##' @importFrom stats spec.pgram kernel .lm.fit
##'
##' @param vars optional; names of observed variables for which the power spectrum will be computed.
##' By default, the spectrum will be computed for all observables.
##' @param kernel.width width parameter for the smoothing kernel used for
##' calculating the estimate of the spectrum.
##' @param nsim number of model simulations to be computed.
##' @param seed optional; if non-\code{NULL}, the random number generator will
##' be initialized with this seed for simulations.
##' See \code{\link[=simulate-pomp]{simulate}}.
##' @param transform.data function; this transformation will be applied to the
##' observables prior to estimation of the spectrum, and prior to any
##' detrending.
##' @param 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.
##' @inheritParams pomp
##'
##' @return
##' An object of class \sQuote{spectd_pomp}, which contains the model, the data, and the results of the \code{spect} computation.
##' The following methods are available:
##' \describe{
##' \item{plot}{produces some diagnostic plots}
##' \item{summary}{displays a summary}
##' \item{logLik}{gives a measure of the agreement of the power spectra}
##' }
##'
##' @author Daniel C. Reuman, Cai GoGwilt, Aaron A. King
##'
##' @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. \emph{Proceedings of the National Academy
##' of Sciences} \bold{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. \emph{Ecology Letters}, \bold{11}, 820-830.
NULL
setClass(
"spectd_pomp",
contains="pomp",
slots=c(
nsim="integer",
seed="integer",
kernel.width="numeric",
transform.data="function",
vars="character",
freq="numeric",
datspec="array",
simspec="array",
pvals="numeric",
detrend="character"
)
)
setGeneric(
"spect",
function (data, ...)
standardGeneric("spect")
)
setMethod(
"spect",
signature=signature(data="missing"),
definition=function (...) {
reqd_arg("spect","data")
}
)
setMethod(
"spect",
signature=signature(data="ANY"),
definition=function (data, ...) {
undef_method("spect",data)
}
)
##' @name spect-data.frame
##' @aliases spect spect,data.frame-method
##' @rdname spect
##' @export
setMethod(
"spect",
signature(data="data.frame"),
function (data,
vars, kernel.width, nsim, seed = NULL, transform.data = identity,
detrend = c("none","mean","linear","quadratic"),
params, rinit, rprocess, rmeasure,
..., verbose = getOption("verbose", FALSE)) {
tryCatch(
spect.internal(
data,
vars=vars,
kernel.width=kernel.width,
nsim=nsim,
seed=seed,
transform.data=match.fun(transform.data),
detrend=match.arg(detrend),
params=params,
rinit=rinit,
rprocess=rprocess,
rmeasure=rmeasure,
...,
verbose=verbose
),
error = function (e) pStop("spect",conditionMessage(e))
)
}
)
##' @name spect-pomp
##' @aliases spect spect,pomp-method
##' @rdname spect
##' @export
setMethod(
"spect",
signature(data="pomp"),
function (data,
vars, kernel.width, nsim, seed = NULL, transform.data = identity,
detrend = c("none","mean","linear","quadratic"),
..., verbose = getOption("verbose", FALSE)) {
tryCatch(
spect.internal(
data,
vars=vars,
kernel.width=kernel.width,
nsim=nsim,
seed=seed,
transform.data=match.fun(transform.data),
detrend=match.arg(detrend),
...,
verbose=verbose
),
error = function (e) pStop("spect",conditionMessage(e))
)
}
)
##' @name spect-spectd_pomp
##' @aliases spect,spectd_pomp-method
##' @rdname spect
##' @export
setMethod(
"spect",
signature=signature(data="spectd_pomp"),
definition=function (data,
vars, kernel.width, nsim, seed = NULL, transform.data,
detrend,
..., verbose = getOption("verbose", FALSE)) {
if (missing(vars)) vars <- colnames(data@datspec)
if (missing(kernel.width)) kernel.width <- data@kernel.width
if (missing(nsim)) nsim <- data@nsim
if (missing(transform.data)) transform.data <- data@transform.data
if (missing(detrend)) detrend <- data@detrend
spect(
as(data,"pomp"),
vars=vars,
kernel.width=kernel.width,
nsim=nsim,
seed=seed,
transform.data=transform.data,
detrend=detrend,
...,
verbose=verbose
)
}
)
spect.internal <- function (object, vars, kernel.width, nsim, seed = NULL,
transform.data, detrend, ..., verbose) {
verbose <- as.logical(verbose)
object <- pomp(object,...,verbose=verbose)
if (undefined(object@rprocess) || undefined(object@rmeasure))
pStop_(paste(sQuote(c("rprocess","rmeasure")),collapse=", ")," are needed basic components.")
if (missing(vars)) vars <- rownames(object@data)
if (missing(kernel.width) || length(kernel.width) > 1 ||
!is.numeric(kernel.width) ||
!is.finite(kernel.width) || kernel.width < 0)
pStop_(sQuote("kernel.width")," must be a positive integer.")
if (missing(nsim) || length(nsim) > 1 || !is.numeric(nsim)||
!is.finite(nsim) || (nsim<1))
pStop_(sQuote("nsim")," must be a positive integer.")
nsim <- as.integer(nsim)
seed <- as.integer(seed)
ker <- reuman.kernel(kernel.width)
params <- coef(object)
pompLoad(object)
on.exit(pompUnload(object))
ds <- compute.spect.data(
object,
vars=vars,
transform.data=transform.data,
detrend=detrend,
ker=ker
)
freq <- ds$freq
datspec <- ds$spec
simspec <- compute.spect.sim(
object,
params=params,
vars=vars,
nsim=nsim,
seed=seed,
transform.data=transform.data,
detrend=detrend,
ker=ker
)
pvals <- numeric(length(vars)+1)
names(pvals) <- c(vars,"all")
mean.simspec <- colMeans(simspec) # mean spectrum of simulations
totdatdist <- 0
totsimdist <- 0
for (j in seq_along(vars)) {
## L-2 distance between data and mean simulated spectrum
datdist <- sum((datspec[,j]-mean.simspec[,j])^2)
## L-2 distance betw. each sim. and mean simulated spectrum
simdist <- vapply(
seq_len(nsim),
function(k)sum((simspec[k,,j]-mean.simspec[,j])^2),
numeric(1)
)
pvals[j] <- (nsim+1-sum(simdist<datdist))/(nsim+1)
totdatdist <- totdatdist+datdist
totsimdist <- totsimdist+simdist
}
pvals[length(vars)+1] <- (nsim+1-sum(totsimdist<totdatdist))/(nsim+1)
coef(object) <- params
new(
"spectd_pomp",
object,
nsim=nsim,
seed=seed,
kernel.width=kernel.width,
transform.data=transform.data,
vars=vars,
detrend=detrend,
freq=freq,
datspec=datspec,
simspec=simspec,
pvals=pvals
)
}
compute.spect.data <- function (object, vars, transform.data, detrend, ker) {
dat <- obs(object,vars)
if (any(!is.finite(dat)))
pStop_("missing or infinite values in the data.")
dt <- diff(time(object,t0=FALSE))
base.freq <- 1/mean(dt)
dt.tol <- 0.025
if (max(dt)-min(dt)>dt.tol*mean(dt))
pStop_(sQuote("spect")," assumes evenly spaced times.")
for (j in seq_along(vars)) {
sp <- spec.pgram(
pomp.detrend(transform.data(dat[j,]),type=detrend),spans=ker,taper=0,
pad=0,fast=FALSE,detrend=FALSE,plot=FALSE
)
if (j==1) {
freq <- base.freq*sp$freq
datspec <- array(
dim=c(length(freq),nrow(dat)),
dimnames=list(NULL,vars)
)
}
datspec[,j] <- log10(sp$spec)
}
list(freq=freq,spec=datspec)
}
compute.spect.sim <- function (object, params, vars, nsim, seed,
transform.data, detrend, ker) {
sims <- tryCatch(
{
s <- freeze(.Call(P_do_simulate,object,params,nsim,0L,gnsi=TRUE),seed=seed)
s$obs[vars,,,drop=FALSE]
},
error = function (e) pStop_("in simulation: ",
conditionMessage(e))
)
if (any(!is.finite(sims)))
pStop_("missing or infinite values in simulated data.")
nobs <- length(vars)
for (j in seq_len(nobs)) {
for (k in seq_len(nsim)) {
sp <- tryCatch(
spec.pgram(pomp.detrend(transform.data(sims[j,k,]),type=detrend),
spans=ker,taper=0,pad=0,fast=FALSE,detrend=FALSE,plot=FALSE),
error = function (e) pStop("spec.pgram",conditionMessage(e))
)
if ((j==1)&&(k==1)) {
simspec <- array(
dim=c(nsim,length(sp$freq),nobs),
dimnames=list(NULL,NULL,vars)
)
}
simspec[k,,j] <- log10(sp$spec)
}
}
simspec
}
## detrends in one of several ways, according to type.
## tseries is a numeric vector,
pomp.detrend <- function (tseries, type) {
switch(
type,
mean=tseries-mean(tseries),
linear={
m <- cbind(1,seq_along(tseries))
.lm.fit(m,tseries)$residuals
},
quadratic={
x <- seq_along(tseries)
m <- cbind(1,x,x*x)
.lm.fit(m,tseries)$residuals
},
tseries
)
}
## The default smoothing kernel for the R spec.pgram function is weird.
## This function creates a better one.
reuman.kernel <- function (kernel.width) {
ker <- kernel("modified.daniell",m=kernel.width,r=NA)
x <- seq.int(from=0,to=kernel.width,by=1)/kernel.width
ker[[1L]] <- (15/(16*2*pi))*((x-1)^2)*((x+1)^2)
ker[[1L]] <- ker[[1L]]/(2*sum(ker[[1L]][-1])+ker[[1L]][1L])
attr(ker,"name") <- NULL
ker
}
##' @name summary-spectd_pomp
##' @aliases summary,spectd_pomp-method
##' @rdname summary
##' @export
setMethod(
"summary",
signature=signature(object="spectd_pomp"),
definition=function (object) {
list(
coef=coef(object),
nsim=nrow(object@simspec),
pvals=object@pvals
)
}
)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.