R/mif1.R

## mif1 algorithm functions

## define the mif1 class
setClass(
  'mif1',
  contains='pfilterd.pomp',
  slots=c(
    transform = "logical",
    ivps = 'character',
    pars = 'character',
    Nmif = 'integer',
    var.factor='numeric',
    ic.lag='integer',
    cooling.type='character',
    cooling.fraction.50='numeric',
    method='character',
    random.walk.sd = 'numeric',
    conv.rec = 'matrix'
  )
)

mif1.particles <- function (Np, center, sd, ...) {
  matrix(
    data=rnorm(
      n=Np*length(center),
      mean=center,
      sd=sd
    ),
    nrow=length(center),
    ncol=Np,
    dimnames=list(
      variable=names(center),
      rep=NULL
    )
  )
}

mif1.cooling <- function (type, fraction, ntimes) {
  switch(
    type,
    geometric={
      factor <- fraction^(1/50)
      function (nt, m) {
        alpha <- factor^(m-1)
        list(alpha=alpha,gamma=alpha^2)
      }
    },
    hyperbolic={
      if (fraction < 1) {
        scal <- (50*fraction-1)/(1-fraction)
        function (nt, m) {
          alpha <- (1+scal)/(scal+m-1)
          list(alpha=alpha,gamma=alpha^2)
        }
      } else {
        function (nt, m) {
          list(alpha=1,gamma=1)
        }
      }
    }
  )
}

mif1.pfilter <- function (object, params, Np,
                         tol, max.fail,
                         pred.mean = FALSE,
                         pred.var = FALSE,
                         cooling, cooling.m,
                         rw.sd,
                         verbose = FALSE,
                         .transform = FALSE,
                         .getnativesymbolinfo = TRUE) {

  object <- as(object,"pomp")
  ep <- paste0("in ",sQuote("mif1.pfilter"),": ")

  pompLoad(object,verbose=verbose)

  ptsi.for <- gnsi.rproc <- gnsi.dmeas <- as.logical(.getnativesymbolinfo)
  pred.mean <- as.logical(pred.mean)
  pred.var <- as.logical(pred.var)
  verbose <- as.logical(verbose)
  transform <- as.logical(.transform)

  times <- time(object,t0=TRUE)
  ntimes <- length(times)-1
  Np <- as.integer(Np)

  x <- init.state(
    object,
    params=if (transform) {
      partrans(object,params,dir="fromEstimationScale",
               .getnativesymbolinfo=ptsi.for)
    } else {
      params
    }
  )
  ptsi.for <- FALSE

  statenames <- rownames(x)
  paramnames <- rownames(params)

  rw.names <- names(rw.sd)
  sigma <- rw.sd
  npars <- length(rw.names)
  nvars <- nrow(x)

  loglik <- rep(NA,ntimes)
  eff.sample.size <- numeric(ntimes)
  nfail <- 0

  ## set up storage for prediction means, variances, etc.
  if (pred.mean) {
    pred.m <- matrix(
      data=0,
      nrow=nvars+npars,
      ncol=ntimes,
      dimnames=list(
        variable=c(statenames,rw.names),
        time=time(object))
    )
  } else {
    pred.m <- array(data=numeric(0),dim=c(0,0))
  }

  if (pred.var) {
    pred.v <- matrix(
      data=0,
      nrow=nvars+npars,
      ncol=ntimes,
      dimnames=list(
        variable=c(statenames,rw.names),
        time=time(object))
    )
  } else {
    pred.v <- array(data=numeric(0),dim=c(0,0))
  }

  filt.m <- matrix(
    data=0,
    nrow=nvars+length(paramnames),
    ncol=ntimes,
    dimnames=list(
      variable=c(statenames,paramnames),
      time=time(object)
    )
  )

  for (nt in seq_len(ntimes)) { ## main loop

    sigma1 <- sigma

    ## transform the parameters if necessary
    if (transform) tparams <- partrans(object,params,dir="fromEstimationScale",
                                       .getnativesymbolinfo=ptsi.for)
    ptsi.for <- FALSE

    ## advance the state variables according to the process model
    X <- tryCatch(
      rprocess(
        object,
        xstart=x,
        times=times[c(nt,nt+1)],
        params=if (transform) tparams else params,
        offset=1,
        .getnativesymbolinfo=gnsi.rproc
      ),
      error = function (e) {
        stop(ep,"process simulation error: ",conditionMessage(e),call.=FALSE)
      }
    )
    gnsi.rproc <- FALSE

    if (pred.var) { ## check for nonfinite state variables and parameters
      problem.indices <- unique(which(!is.finite(X),arr.ind=TRUE)[,1L])
      if (length(problem.indices)>0) {  # state variables
        stop(
          ep,"non-finite state variable(s): ",
          paste(rownames(X)[problem.indices],collapse=', '),
          call.=FALSE
        )
      }
      problem.indices <- unique(which(!is.finite(params[rw.names,,drop=FALSE]),arr.ind=TRUE)[,1L])
      if (length(problem.indices)>0) {
        stop(
          ep,"non-finite parameter(s): ",
          paste(rw.names[problem.indices],collapse=', '),
          call.=FALSE
        )
      }
    }

    ## determine the weights
    weights <- tryCatch(
      dmeasure(
        object,
        y=object@data[,nt,drop=FALSE],
        x=X,
        times=times[nt+1],
        params=if (transform) tparams else params,
        log=FALSE,
        .getnativesymbolinfo=gnsi.dmeas
      ),
      error = function (e) {
        stop(ep,"error in calculation of weights.",
             conditionMessage(e),call.=FALSE)
      }
    )
    if (!all(is.finite(weights))) {
      first <- which(!is.finite(weights))[1L]
      datvals <- object@data[,nt]
      weight <- weights[first]
      states <- X[,first,1L]
      params <- if (transform) tparams[,first] else params[,first]
      msg <- nonfinite_dmeasure_error(time=times[nt+1],lik=weight,datvals,states,params)
      stop(ep,msg,call.=FALSE)
    }
    gnsi.dmeas <- FALSE

    ## compute prediction mean, prediction variance, filtering mean,
    ## effective sample size, log-likelihood
    ## also do resampling if filtering has not failed
    xx <- tryCatch(
      .Call(
        pfilter3_computations,
        x=X,
        params=params,
        Np=Np[nt+1],
        rw_sd=sigma1,
        predmean=pred.mean,
        predvar=pred.var,
        filtmean=TRUE,
        trackancestry=FALSE,
        onepar=FALSE,
        weights=weights,
        tol=tol
      ),
      error = function (e) {
        stop(ep,"particle-filter error: ",conditionMessage(e),call.=FALSE) # nocov
      }
    )
    all.fail <- xx$fail
    loglik[nt] <- xx$loglik
    eff.sample.size[nt] <- xx$ess

    if (pred.mean)
      pred.m[,nt] <- xx$pm
    if (pred.var)
      pred.v[,nt] <- xx$pv
    filt.m[,nt] <- xx$fm

    x <- xx$states
    params <- .Call(randwalk_perturbation,xx$params,sigma1)

    if (all.fail) { ## all particles are lost
      nfail <- nfail+1
      if (verbose)
        message("filtering failure at time t = ",times[nt+1])
      if (nfail>max.fail)
        stop(ep,"too many filtering failures.",call.=FALSE)
    } else {
      if (pred.var)
        pred.v[rw.names,nt] <- pred.v[rw.names,nt]+sigma1^2
    }

    if (verbose && (nt%%5==0))
      cat("pfilter timestep",nt,"of",ntimes,"finished\n")

  } ## end of main loop

  if (nfail>0) {
    warning(
      ep,nfail,
      ngettext(
        nfail,
        msg1=" filtering failure occurred.",
        msg2=" filtering failures occurred."
      ),
      call.=FALSE
    )
  }

  pompUnload(object,verbose=verbose)

  new(
    "pfilterd.pomp",
    object,
    pred.mean=pred.m,
    pred.var=pred.v,
    filter.mean=filt.m,
    paramMatrix=array(data=numeric(0),dim=c(0,0)),
    eff.sample.size=eff.sample.size,
    cond.loglik=loglik,
    Np=as.integer(Np),
    tol=tol,
    nfail=as.integer(nfail),
    loglik=sum(loglik)
  )
}

mif1.internal <- function (object, Nmif,
                          start, ivps,
                          rw.sd,
                          Np, var.factor, ic.lag,
                          cooling.type, cooling.fraction.50,
                          method,
                          tol, max.fail,
                          verbose, transform, .ndone = 0L,
                          paramMatrix = NULL,
                          .getnativesymbolinfo = TRUE,
                          ...) {

  ep <- paste0("in ",sQuote("mif1"),": ")

  if (method=="mif2") {
    stop(ep,"method=",sQuote("mif2")," has been removed.\n",
         "Use ",sQuote("mif2")," instead.",call.=FALSE)
  }

  pompLoad(object,verbose=verbose)

  gnsi <- as.logical(.getnativesymbolinfo)

  transform <- as.logical(transform)

  if (length(start)==0)
    stop(
      ep,sQuote("start"),
      " must be specified if ",
      sQuote("coef(object)")," is NULL",
      call.=FALSE
    )

  if (transform)
    start <- partrans(object,start,dir="toEstimationScale")

  start.names <- names(start)
  if (is.null(start.names))
    stop(ep,sQuote("start"),
         " must be a named vector",call.=FALSE)

  rw.names <- names(rw.sd)
  if (is.null(rw.names) || any(rw.sd<0))
    stop(ep,sQuote("rw.sd"),
         " must be a named non-negative numerical vector",call.=FALSE)
  if (!all(rw.names%in%start.names))
    stop(ep,"all the names of ",sQuote("rw.sd"),
         " must be names of ",sQuote("start"),call.=FALSE)
  rw.names <- rw.names[rw.sd>0]
  rw.sd <- rw.sd[rw.sd>0]
  if (length(rw.names) == 0)
    stop(ep,sQuote("rw.sd"),
         " must have one positive entry for each parameter to be estimated",call.=FALSE)

  pars <- rw.names[!(rw.names%in%ivps)]

  if (!is.character(ivps) || !all(ivps%in%start.names))
    stop(ep,sQuote("ivps"),
         " must name model parameters",call.=FALSE)

  ntimes <- length(time(object))
  if (is.null(Np)) stop(ep,sQuote("Np"),
                        " must be specified",call.=FALSE)
  if (is.function(Np)) {
    Np <- tryCatch(
      vapply(seq.int(from=0,to=ntimes,by=1),Np,numeric(1)),
      error = function (e) {
        stop(ep,"if ",sQuote("Np")," is a function, ",
             "it must return a single positive integer",call.=FALSE)
      }
    )
  }
  if (length(Np)==1)
    Np <- rep(Np,times=ntimes+1)
  else if (length(Np)!=(ntimes+1))
    stop(ep,sQuote("Np")," must have length 1 or length ",ntimes+1,call.=FALSE)
  if (any(Np<=0))
    stop(ep,"number of particles, ",sQuote("Np"),", must always be positive",call.=FALSE)
  if (!is.numeric(Np))
    stop(ep,sQuote("Np"),
         " must be a number, a vector of numbers, or a function",
         call.=FALSE)
  Np <- as.integer(Np)

  ic.lag <- as.integer(ic.lag)
  if ((length(ic.lag)!=1)||(ic.lag<1))
    stop(ep,sQuote("ic.lag"),
         " must be a positive integer",call.=FALSE)
  if (ic.lag>ntimes) {
    warning(
      ep,sQuote("ic.lag")," = ",ic.lag," > ",ntimes,
      " = length(time(",sQuote("object"),"))",
      " is nonsensical.  Setting ",sQuote("ic.lag")," = ",ntimes,".",
      call.=FALSE
    )
    ic.lag <- length(time(object))
  }
  if ((length(pars)==0)&&(ic.lag<length(time(object)))) {
    warning(
      ep,"only IVPs are to be estimated, yet ",
      sQuote("ic.lag")," = ",ic.lag,
      " < ",ntimes," = length(time(",sQuote("object"),")),",
      " so unnecessary work is to be done.",
      call.=FALSE
    )
  }

  if (missing(cooling.fraction.50))
    stop(ep,sQuote("cooling.fraction.50"),
         " must be specified",call.=FALSE)
  cooling.fraction.50 <- as.numeric(cooling.fraction.50)
  if ((length(cooling.fraction.50)!=1)||(cooling.fraction.50<0)||(cooling.fraction.50>1))
    stop(ep,sQuote("cooling.fraction.50"),
         " must be a number between 0 and 1",call.=FALSE)

  cooling <- mif1.cooling(
    type=cooling.type,
    fraction=cooling.fraction.50,
    ntimes=ntimes
  )

  if ((length(var.factor)!=1)||(var.factor < 0))
    stop(ep,sQuote("var.factor")," must be a positive number",call.=FALSE)

  Nmif <- as.integer(Nmif)
  if (Nmif<0)
    stop(ep,sQuote("Nmif")," must be a positive integer",call.=FALSE)

  theta <- start

  sigma <- rep(0,length(start))
  names(sigma) <- start.names

  rw.sd <- rw.sd[c(pars,ivps)]
  rw.names <- names(rw.sd)

  sigma[rw.names] <- rw.sd

  conv.rec <- matrix(
    data=NA,
    nrow=Nmif+1,
    ncol=length(theta)+2,
    dimnames=list(
      iteration=seq(.ndone,.ndone+Nmif),
      variable=c('loglik','nfail',names(theta))
    )
  )
  conv.rec[1L,] <- c(NA,NA,theta)

  if (!all(is.finite(theta[c(pars,ivps)]))) {
    stop(
      ep,sQuote("mif")," cannot estimate non-finite parameters.\n",
      "The following ",if (transform) "transformed ", "parameters are non-finite: ",
      paste(
        sQuote(c(pars,ivps)[!is.finite(theta[c(pars,ivps)])]),
        collapse=","
      ),
      call.=FALSE
    )
  }

  pfp <- as(object,"pomp")

  for (n in seq_len(Nmif)) { ## iterate the filtering

    ## get the intensity of artificial noise from the cooling schedule
    cool.sched <- cooling(nt=1,m=.ndone+n)
    sigma.n <- sigma*cool.sched$alpha

    ## initialize the parameter portions of the particles
    P <- mif1.particles(
      Np=Np[1L],
      center=theta,
      sd=sigma.n*var.factor
    )

    pfp <- tryCatch(
      mif1.pfilter(
        object=pfp,
        params=P,
        Np=Np,
        tol=tol,
        max.fail=max.fail,
        pred.mean=(n==Nmif),
        pred.var=((method=="mif")||(n==Nmif)),
        cooling=cooling,
        cooling.m=.ndone+n,
        rw.sd=sigma.n[pars],
        .transform=transform,
        verbose=verbose,
        .getnativesymbolinfo=gnsi
      ),
      error = function (e) {
        stop(ep,conditionMessage(e),call.=FALSE)
      }
    )

    gnsi <- FALSE

    ## update parameters
    switch(
      method,
      mif={                                 # mif update rule
        v <- pfp$pred.var[pars,,drop=FALSE] # the prediction variance
        v1 <- cool.sched$gamma*(1+var.factor^2)*sigma[pars]^2
        theta.hat <- cbind(theta[pars],pfp$filter.mean[pars,,drop=FALSE])
        theta[pars] <- theta[pars]+colSums(apply(theta.hat,1,diff)/t(v))*v1
      },
      unweighted={                 # unweighted average
        theta[pars] <- rowMeans(pfp@filter.mean[pars,,drop=FALSE])
      },
      fp={                         # fixed-point iteration
        theta[pars] <- pfp@filter.mean[pars,ntimes,drop=FALSE]
      },
      stop(ep,"unrecognized method ",sQuote(method),call.=FALSE)
    )
    theta[ivps] <- pfp@filter.mean[ivps,ic.lag]
    conv.rec[n+1,-c(1,2)] <- theta
    conv.rec[n,c(1,2)] <- c(pfp@loglik,pfp@nfail)

    if (verbose) cat("mif iteration",n,"of",Nmif,"completed\n")

  } ### end of main loop

  ## back transform the parameter estimate if necessary
  if (transform) theta <- partrans(pfp,theta,dir="fromEstimationScale")

  pompUnload(object,verbose=verbose)

  new(
    "mif1",
    pfp,
    transform=transform,
    params=theta,
    ivps=ivps,
    pars=pars,
    Nmif=Nmif,
    var.factor=var.factor,
    ic.lag=ic.lag,
    random.walk.sd=sigma[rw.names],
    tol=tol,
    conv.rec=conv.rec,
    method=method,
    cooling.type=cooling.type,
    cooling.fraction.50=cooling.fraction.50,
    paramMatrix=array(data=numeric(0),dim=c(0,0))
  )
}

setMethod(
  "mif1",
  signature=signature(object="pomp"),
  function (object, Nmif = 1,
            start,
            ivps = character(0),
            rw.sd, Np, ic.lag, var.factor = 1,
            cooling.type = c("geometric","hyperbolic"),
            cooling.fraction.50,
            method = c("mif","unweighted","fp","mif2"),
            tol = 1e-17, max.fail = Inf,
            verbose = getOption("verbose"),
            transform = FALSE,
            ...) {

    ep <- paste0("in ",sQuote("mif"),": ")

    method <- match.arg(method)

    if (missing(start)) start <- coef(object)
    if (missing(rw.sd))
      stop(ep,sQuote("rw.sd")," must be specified.",call.=FALSE)
    if (!is.numeric(rw.sd))
      stop(ep,sQuote("rw.sd")," must be a named numeric vector.",call.=FALSE)
    if (missing(ic.lag)) {
      if (length(ivps)>0) {
        stop(ep,sQuote("ic.lag"),
             " must be specified if ",sQuote("ivps"),
             " are.",call.=FALSE)
      } else {
        ic.lag <- length(time(object))
      }
    }

    if (missing(Np))
      stop(ep,sQuote("Np")," must be specified",call.=FALSE)

    cooling.type <- match.arg(cooling.type)

    mif1.internal(
      object=object,
      Nmif=Nmif,
      start=start,
      ivps=ivps,
      rw.sd=rw.sd,
      Np=Np,
      cooling.type=cooling.type,
      cooling.fraction.50=cooling.fraction.50,
      var.factor=var.factor,
      ic.lag=ic.lag,
      method=method,
      tol=tol,
      max.fail=max.fail,
      verbose=verbose,
      transform=transform,
      ...
    )

  }
)


setMethod(
  "mif1",
  signature=signature(object="pfilterd.pomp"),
  function (object, Nmif = 1, Np, tol,
            ...) {

    if (missing(Np)) Np <- object@Np
    if (missing(tol)) tol <- object@tol

    f <- selectMethod("mif1","pomp")
    f(object=object,Nmif=Nmif,Np=Np,tol=tol,...)
  }
)

setMethod(
  "mif1",
  signature=signature(object="mif1"),
  function (object, Nmif,
            start, ivps, rw.sd,
            Np, ic.lag, var.factor,
            cooling.type, cooling.fraction.50,
            method,
            tol,
            transform,
            ...) {

    if (missing(Nmif)) Nmif <- object@Nmif
    if (missing(start)) start <- coef(object)
    if (missing(ivps)) ivps <- object@ivps
    if (missing(rw.sd)) rw.sd <- object@random.walk.sd
    if (missing(ic.lag)) ic.lag <- object@ic.lag
    if (missing(var.factor)) var.factor <- object@var.factor
    if (missing(cooling.type)) cooling.type <- object@cooling.type
    if (missing(cooling.fraction.50)) cooling.fraction.50 <- object@cooling.fraction.50
    if (missing(method)) method <- object@method
    if (missing(transform)) transform <- object@transform

    if (missing(Np)) Np <- object@Np
    if (missing(tol)) tol <- object@tol

    f <- selectMethod("mif1","pomp")
    f(object=object,
      Nmif=Nmif,
      start=start,
      ivps=ivps,
      rw.sd=rw.sd,
      Np=Np,
      cooling.type=cooling.type,
      cooling.fraction.50=cooling.fraction.50,
      var.factor=var.factor,
      ic.lag=ic.lag,
      method=method,
      tol=tol,
      transform=transform,
      ...
    )
  }
)

setMethod(
  'continue',
  signature=signature(object='mif1'),
  function (object, Nmif = 1,
            ...) {

    ndone <- object@Nmif

    f <- selectMethod("mif","mif")

    obj <- f(
      object=object,
      Nmif=Nmif,
      .ndone=ndone,
      paramMatrix=object@paramMatrix,
      ...
    )

    object@conv.rec[ndone+1,c('loglik','nfail')] <- obj@conv.rec[1L,c('loglik','nfail')]
    obj@conv.rec <- rbind(
      object@conv.rec,
      obj@conv.rec[-1L,colnames(object@conv.rec)]
    )
    names(dimnames(obj@conv.rec)) <- c("iteration","variable")
    obj@Nmif <- as.integer(ndone+Nmif)

    obj
  }
)
nxdao2000/is2 documentation built on May 24, 2019, 11:50 a.m.