R/spect_match.R

In kingaa/pomp: Statistical Inference for Partially Observed Markov Processes

##' Spectrum matching
##'
##' Estimation of parameters by matching power spectra
##'
##' In spectrum matching, one attempts to minimize the discrepancy between a \acronym{POMP} model's predictions and data, as measured in the frequency domain by the power spectrum.
##'
##' \code{spect_objfun} constructs an objective function that measures the discrepancy.
##' It can be passed to any one of a variety of numerical optimization routines, which will adjust model parameters to minimize the discrepancies between the power spectrum of model simulations and that of the data.
##' 
##' @docType methods
##' @name spect_match
##' @rdname spect_match
##' @aliases spect_objfun spect_objfun,missing-method spect_objfun,ANY-method
##' @concept power-spectrum matching
##' @family estimation methods
##' @family summary statistic-based methods
##' @family methods based on maximization
##' @include spect.R probe_match.R loglik.R plot.R
##' @param weights optional numeric or function.
##' The mismatch between model and data is measured by a weighted average of mismatch at each frequency.
##' By default, all frequencies are weighted equally.
##' \code{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, \code{weights(freq)} must return a nonnegative weight for each frequency.
##' @inheritParams probe_match
##' @inheritParams spect
##' @inheritParams pomp
##' @return
##' \code{spect_objfun} constructs a stateful objective function for spectrum matching.
##' Specifically, \code{spect_objfun} returns an object of class \sQuote{spect_match_objfun}, which is a function suitable for use in an \code{\link[stats]{optim}}-like optimizer.
##' This function takes a single numeric-vector argument that is assumed to contain the parameters named in \code{est}, in that order.
##' When called, it will return the (optionally weighted) \eqn{L^2}{L2} distance between the data spectrum and simulated spectra.
##' It is a stateful function:
##' Each time it is called, it will remember the values of the parameters and the discrepancy measure.
##' @references
##'
##' \Reuman2006
##'
##' \Reuman2008
##' @inheritSection pomp Note for Windows users
##' @inheritSection objfun Important Note
##' @inheritSection objfun Warning! Objective functions based on C snippets
##' @seealso \code{\link{spect}} \code{\link[stats]{optim}}
##' \code{\link[subplex]{subplex}} \code{\link[nloptr]{nloptr}}
##' 
##' @example examples/spect_match.R
##'
NULL

setClass(
  "spect_match_objfun",
  contains="function",
  slots=c(
    env="environment",
    est="character"
  )
)

setGeneric(
  "spect_objfun",
  function (data, ...)
    standardGeneric("spect_objfun")
)

setMethod(
  "spect_objfun",
  signature=signature(data="missing"),
  definition=function (...) {
    reqd_arg("spect_objfun","data")
  }
)

setMethod(
  "spect_objfun",
  signature=signature(data="ANY"),
  definition=function (data, ...) {
    undef_method("spect_objfun",data)
  }
)

##' @rdname spect_match
##' @export
setMethod(
  "spect_objfun",
  signature=signature(data="data.frame"),
  definition=function(data, est = character(0),
    weights = 1, fail.value = NA,
    vars, kernel.width, nsim, seed = NULL, transform.data = identity,
    detrend = c("none","mean","linear","quadratic"),
    params, rinit, rprocess, rmeasure, partrans,
    ..., verbose = getOption("verbose", FALSE)) {

    tryCatch(
      smof_internal(
        data,
        est=est,
        weights=weights,
        fail.value=fail.value,
        vars=vars,
        kernel.width=kernel.width,
        nsim=nsim,
        seed=seed,
        transform.data=transform.data,
        detrend=detrend,
        params=params,
        rinit=rinit,
        rprocess=rprocess,
        rmeasure=rmeasure,
        partrans=partrans,
        ...,
        verbose=verbose
      ),
      error = function (e) pStop(who="spect_objfun",conditionMessage(e))
    )

  }
)

##' @rdname spect_match
##' @export
setMethod(
  "spect_objfun",
  signature=signature(data="pomp"),
  definition=function(data, est = character(0),
    weights = 1, fail.value = NA,
    vars, kernel.width, nsim, seed = NULL, transform.data = identity,
    detrend = c("none","mean","linear","quadratic"),
    ..., verbose = getOption("verbose", FALSE)) {

    tryCatch(
      smof_internal(
        data,
        est=est,
        weights=weights,
        fail.value=fail.value,
        vars=vars,
        kernel.width=kernel.width,
        nsim=nsim,
        seed=seed,
        transform.data=transform.data,
        detrend=detrend,
        ...,
        verbose=verbose
      ),
      error = function (e) pStop(who="spect_objfun",conditionMessage(e))
    )

  }
)

##' @rdname spect_match
##' @export
setMethod(
  "spect_objfun",
  signature=signature(data="spectd_pomp"),
  definition=function(data, est = character(0),
    weights = 1, fail.value = NA,
    vars, kernel.width, nsim, seed = NULL, transform.data = identity,
    detrend,
    ..., verbose = getOption("verbose", FALSE)) {

    if (missing(vars)) vars <- data@vars
    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_objfun(
      as(data,"pomp"),
      est=est,
      weights=weights,
      fail.value=fail.value,
      vars=vars,
      kernel.width=kernel.width,
      nsim=nsim,
      seed=seed,
      transform.data=transform.data,
      detrend=detrend,
      ...,
      verbose=verbose
    )

  }
)

##' @rdname spect_match
##' @export
setMethod(
  "spect_objfun",
  signature=signature(data="spect_match_objfun"),
  definition=function(data, est, weights, fail.value,
    seed = NULL, ..., verbose = getOption("verbose", FALSE)) {

    if (missing(est)) est <- data@est
    if (missing(weights)) weights <-data@env$weights
    if (missing(fail.value)) fail.value <- data@env$fail.value

    spect_objfun(
      data@env$object,
      est=est,
      weights=weights,
      fail.value=fail.value,
      seed=seed,
      ...,
      verbose=verbose
    )

  }
)

smof_internal <- function (object,
  est, weights, fail.value,
  vars, kernel.width, nsim, seed, transform.data, detrend,
  ..., verbose) {

  verbose <- as.logical(verbose)

  object <- spect(object,vars=vars,kernel.width=kernel.width,
    nsim=nsim,seed=seed, transform.data=transform.data,detrend=detrend,
    ...,verbose=verbose)

  fail.value <- as.numeric(fail.value)

  if (is.numeric(weights)) {
    if (length(weights)==1) {
      weights <- rep(weights,length(object@freq))
    } else if ((length(weights) != length(object@freq)))
      pStop_("if ",sQuote("weights"),
        " is provided as a vector, it must have length ",
        length(object@freq))
  } else if (is.function(weights)) {
    weights <- tryCatch(
      vapply(object@freq,weights,numeric(1)),
      error = function (e)
        pStop_(sQuote("weights")," function: ",conditionMessage(e))
    )
  } else {
    pStop_(sQuote("weights"),
      " must be specified as a vector or as a function")
  }

  if (any(!is.finite(weights) | weights<0))
    pStop_(sQuote("weights")," should be nonnegative and finite")
  weights <- weights/mean(weights)

  params <- coef(object,transform=TRUE)

  est <- as.character(est)
  est <- est[nzchar(est)]

  idx <- match(est,names(params))
  if (any(is.na(idx))) {
    missing <- est[is.na(idx)]
    pStop_(ngettext(length(missing),"parameter","parameters")," ",
      paste(sQuote(missing),collapse=","),
      " not found in ",sQuote("params"))
  }

  pompLoad(object)

  ker <- reuman_kernel(kernel.width)
  discrep <- spect_discrep(object,ker=ker,weights=weights)
  .gnsi <- TRUE

  ofun <- function (par = numeric(0)) {
    params[idx] <- par
    coef(object,transform=TRUE,.gnsi=.gnsi) <<- params
    object@simspec <- compute_spect_sim(
      object,
      vars=object@vars,
      params=object@params,
      nsim=object@nsim,
      seed=object@seed,
      transform.data=object@transform.data,
      detrend=object@detrend,
      ker=ker,
      .gnsi=.gnsi
    )
    discrep <<- spect_discrep(object,ker=ker,weights=weights)
    .gnsi <<- FALSE
    if (is.finite(discrep) || is.na(fail.value)) discrep else fail.value #nocov
  }

  environment(ofun) <- list2env(
    list(object=object,fail.value=fail.value,
      params=params,idx=idx,discrep=discrep,seed=seed,ker=ker,
      .gnsi=.gnsi,weights=weights),
    parent=parent.frame(2)
  )

  new("spect_match_objfun",ofun,env=environment(ofun),est=est)

}

## compute a measure of the discrepancies between simulations and data
spect_discrep <- function (object, ker, weights) {
  discrep <- array(dim=c(length(object@freq),length(object@vars)))
  sim.means <- colMeans(object@simspec)
  for (j in seq_along(object@freq)) {
    for (k in seq_along(object@vars)) {
      discrep[j,k] <- ((object@datspec[j,k]-sim.means[j,k])^2)/
        mean((object@simspec[,j,k]-sim.means[j,k])^2)
    }
    discrep[j,] <- weights[j]*discrep[j,]
  }
  sum(discrep)
}

##' @rdname spect
##' @details
##' When \code{spect} operates on a spectrum-matching objective function (a \sQuote{spect_match_objfun} object), by default, the
##' random-number generator seed is fixed at the value given when the objective function was constructed.
##' Specifying \code{NULL} or an integer for \code{seed} overrides this behavior.
##' @export
setMethod(
  "spect",
  signature=signature(data="spect_match_objfun"),
  definition=function (data, seed,
    ..., verbose=getOption("verbose", FALSE)) {

    if (missing(seed)) seed <- data@env$seed

    spect(
      data@env$object,
      seed=seed,
      ...,
      verbose=verbose
    )

  }
)

setAs(
  from="spect_match_objfun",
  to="spectd_pomp",
  def = function (from) {
    from@env$object
  }
)

##' @rdname plot
##' @export
setMethod(
  "plot",
  signature=signature(x="spect_match_objfun"),
  definition=function (x, ...) {
    plot(as(x,"spectd_pomp"),...)
  }
)