[Pomp-commits] r1171 - pkg
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jun 4 14:22:29 CEST 2015
Author: kingaa
Date: 2015-06-04 14:22:28 +0200 (Thu, 04 Jun 2015)
New Revision: 1171
Added:
pkg/mif2.R
Log:
- work on integrating mif2 into pomp: stage 1
Added: pkg/mif2.R
===================================================================
--- pkg/mif2.R (rev 0)
+++ pkg/mif2.R 2015-06-04 12:22:28 UTC (rev 1171)
@@ -0,0 +1,520 @@
+## MIF2 algorithm functions
+
+## define a class of perturbation magnitudes
+setClass(
+ "mif2.perturb.sd",
+ slots=c(
+ sds="list",
+ cooling.fraction.50="numeric",
+ cooling.type="character",
+ ivps="character",
+ rwnames="character",
+ transform="logical"
+ ),
+ prototype=prototype(
+ sds=list(),
+ cooling.fraction.50=1.0,
+ cooling.type="hyperbolic",
+ ivps=character(0),
+ rwnames=character(0)
+ )
+ )
+
+mif2.sd <- function (..., .cooling.fraction.50 = 1,
+ .cooling.type = c("hyperbolic","geometric"),
+ .ivps = character(0), transform = FALSE) {
+ sds <- list(...)
+ if (length(sds)==0)
+ stop(sQuote("mif2.sd")," error: at least one parameter should be perturbed",
+ call.=FALSE)
+ rwnames <- names(sds)
+ if (is.null(rwnames) || any(names(rwnames)==""))
+ stop(sQuote("mif2.sd")," accepts only named arguments",call.=FALSE)
+ .cooling.fraction.50 <- as.numeric(.cooling.fraction.50)
+ if (.cooling.fraction.50 <= 0 || .cooling.fraction.50 > 1)
+ stop(sQuote(".cooling.fraction.50")," must be in (0,1]",call.=FALSE)
+ .cooling.type <- match.arg(.cooling.type)
+ .ivps <- as.character(.ivps)
+ transform <- as.logical(transform)
+ new("mif2.perturb.sd",sds=sds,cooling.fraction.50=.cooling.fraction.50,
+ cooling.type=.cooling.type,ivps=.ivps,rwnames=rwnames,transform=transform)
+}
+
+setMethod("cooling_fn",signature=signature(object="mif2.perturb.sd"),
+ definition=function (object, paramnames) {
+ if (!all(object at rwnames %in% paramnames)) {
+ unrec <- object at rwnames[!object at rwnames %in% paramnames]
+ stop(sQuote("mif2")," error: the following parameter(s), ",
+ "which are supposed to be estimated, are not present: ",
+ paste(sapply(sQuote,unrec),collapse=","),
+ call.=FALSE)
+ }
+ for (p in object at rwnames) {
+ sds[[p]] <- match.fun(sds[[p]])
+ }
+ })
+
+## define the mif2d.pomp class
+setClass(
+ 'mif2d.pomp',
+ contains='pfilterd.pomp',
+ slots=c(
+ transform = "logical",
+ Nmif = 'integer',
+ perturbn = 'function',
+ pertsd = 'mif2.perturb.sd',
+ conv.rec = 'matrix'
+ )
+ )
+
+mif2.pfilter <- function (object, params, Np,
+ tol, max.fail,
+ pred.mean, pred.var, filter.mean,
+ mifiter, perturb.fn,
+ transform, verbose,
+ .getnativesymbolinfo = TRUE) {
+
+ ptsi.inv <- ptsi.for <- gnsi.rproc <- gnsi.dmeas <- as.logical(.getnativesymbolinfo)
+ transform <- as.logical(transform)
+ mifiter <- as.integer(mifiter)
+ Np <- as.integer(Np)
+ Np <- c(Np,Np[1])
+
+ times <- time(object,t0=TRUE)
+ ntimes <- length(times)-1
+
+ paramnames <- rownames(params)
+ npars <- nrow(params)
+
+ ## perturb parameters
+ params <- apply(params,2L,perturb.fn,mifiter=mifiter,timeno=1L)
+
+ ## transform parameters if necessary
+ if (transform) {
+ tparams <- partrans(object,params,dir="forward",
+ .getnativesymbolinfo=ptsi.for)
+ ptsi.for <- FALSE
+ }
+
+ ## get initial states
+ x <- init.state(
+ object,
+ params=if (transform) tparams else params,
+ )
+ statenames <- rownames(x)
+ 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(c(statenames,paramnames),NULL)
+ )
+ 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(c(statenames,paramnames),NULL)
+ )
+ else
+ pred.v <- array(data=numeric(0),dim=c(0,0))
+
+ if (filter.mean)
+ filt.m <- matrix(
+ data=0,
+ nrow=nvars+length(paramnames),
+ ncol=ntimes,
+ dimnames=list(c(statenames,paramnames),NULL)
+ )
+ else
+ filt.m <- array(data=numeric(0),dim=c(0,0))
+
+ for (nt in seq_len(ntimes)) {
+
+ if (nt > 1) {
+ ## perturb parameters
+ params <- apply(params,2L,perturb.fn,mifiter=mifiter,timeno=nt)
+ ## transform parameters if necessary
+ if (transform) {
+ tparams <- partrans(object,params,dir="forward",
+ .getnativesymbolinfo=ptsi.for)
+ }
+ }
+
+ ## advance the state variables according to the process model
+ X <- try(
+ rprocess(
+ object,
+ xstart=x,
+ times=times[c(nt,nt+1)],
+ params=if (transform) tparams else params,
+ offset=1,
+ .getnativesymbolinfo=gnsi.rproc
+ ),
+ silent=FALSE
+ )
+ if (inherits(X,'try-error'))
+ stop(sQuote("mif2.pfilter")," error: process simulation error")
+ gnsi.rproc <- FALSE
+
+ ## determine the weights
+ weights <- try(
+ dmeasure(
+ object,
+ y=object at data[,nt,drop=FALSE],
+ x=X,
+ times=times[nt+1],
+ params=if (transform) tparams else params,
+ log=FALSE,
+ .getnativesymbolinfo=gnsi.dmeas
+ ),
+ silent=FALSE
+ )
+ if (inherits(weights,'try-error'))
+ stop(sQuote("mif2.pfilter")," error: error in calculation of weights")
+ if (any(!is.finite(weights))) {
+ stop(sQuote("mif2.pfilter")," error: ",sQuote("dmeasure"),
+ " returns non-finite value")
+ }
+ gnsi.dmeas <- FALSE
+
+ ## compute prediction mean, prediction variance, filtering mean,
+ ## effective sample size, log-likelihood
+ ## also do resampling if filtering has not failed
+ xx <- try(
+ .Call(
+ pomp:::pfilter_computations,
+ X,params,Np[nt+1],
+ FALSE,numeric(0),
+ pred.mean,pred.var,filter.mean,
+ FALSE,weights,tol
+ ),
+ silent=FALSE
+ )
+ if (inherits(xx,'try-error')) {
+ stop(sQuote("mif2.pfilter")," error",call.=FALSE)
+ }
+ all.fail <- xx$fail
+ loglik[nt] <- xx$loglik
+ eff.sample.size[nt] <- xx$ess
+
+ x <- xx$states
+ params <- xx$params
+
+ if (pred.mean)
+ pred.m[,nt] <- xx$pm
+ if (pred.var)
+ pred.v[,nt] <- xx$pv
+ if (filter.mean)
+ filt.m[,nt] <- xx$fm
+
+ 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(sQuote("mif2.pfilter")," error: too many filtering failures",call.=FALSE)
+ }
+
+ if (verbose && (nt%%5==0))
+ cat("mif2 pfilter timestep",nt,"of",ntimes,"finished\n")
+
+ }
+
+ if (nfail>0)
+ warning(sprintf(ngettext(nfail,msg1="%d filtering failure occurred in ",
+ msg2="%d filtering failures occurred in "),nfail),
+ sQuote("mif2.pfilter"),call.=FALSE)
+
+ new(
+ "pfilterd.pomp",
+ object,
+ pred.mean=pred.m,
+ pred.var=pred.v,
+ filter.mean=filt.m,
+ paramMatrix=params,
+ eff.sample.size=eff.sample.size,
+ cond.loglik=loglik,
+ saved.states=list(),
+ saved.params=list(),
+ seed=as.integer(NULL),
+ Np=as.integer(head(Np,-1)),
+ tol=tol,
+ nfail=as.integer(nfail),
+ loglik=sum(loglik)
+ )
+}
+
+mif2.internal <- function (object, Nmif,
+ start, Np, perturb.fn,
+ tol, max.fail,
+ verbose, transform, .ndone = 0L,
+ .getnativesymbolinfo = TRUE,
+ ...) {
+
+ gnsi <- as.logical(.getnativesymbolinfo)
+ transform <- as.logical(transform)
+
+ Nmif <- as.integer(Nmif)
+ if (Nmif<0)
+ stop("mif error: ",sQuote("Nmif")," must be a positive integer",call.=FALSE)
+
+ if (transform)
+ start <- partrans(object,start,dir="inverse")
+
+ ntimes <- length(time(object))
+
+ start.names <- names(start)
+ if (is.null(start.names))
+ stop("mif2 error: ",sQuote("start")," must be a named vector")
+
+ if (!is.function(perturb.fn))
+ stop("mif2 error: ",sQuote("perturb.fn")," must be a function")
+
+ if (is.function(Np)) {
+ Np <- try(
+ vapply(seq.int(from=1,to=ntimes,by=1),Np,numeric(1)),
+ silent=FALSE
+ )
+ if (inherits(Np,"try-error"))
+ stop("if ",sQuote("Np")," is a function, it must return a single positive integer")
+ }
+ if (length(Np)==1)
+ Np <- rep(Np,times=ntimes)
+ else if (length(Np)!=ntimes)
+ stop(sQuote("Np")," must have length 1 or length ",ntimes)
+ if (any(Np<=0))
+ stop("number of particles, ",sQuote("Np"),", must always be positive")
+ if (!is.numeric(Np))
+ stop(sQuote("Np")," must be a number, a vector of numbers, or a function")
+ Np <- as.integer(Np)
+
+ conv.rec <- matrix(
+ data=NA,
+ nrow=Nmif+1,
+ ncol=length(start)+2,
+ dimnames=list(
+ seq(.ndone,.ndone+Nmif),
+ c('loglik','nfail',names(start))
+ )
+ )
+ conv.rec[1L,] <- c(NA,NA,start)
+
+ if (.ndone > 0) { # call is from 'continue'
+ paramMatrix <- object at paramMatrix
+ } else if (Nmif > 0) { # initial call
+ paramMatrix <- array(data=start,dim=c(length(start),Np[1L]),
+ dimnames=list(names(start),NULL))
+ } else { # no work to do
+ paramMatrix <- array(data=numeric(0),dim=c(0,0))
+ }
+
+ object <- as(object,"pomp")
+
+ for (n in seq_len(Nmif)) { ## iterate the filtering
+
+ pfp <- try(
+ mif2.pfilter(
+ object=object,
+ params=paramMatrix,
+ Np=Np,
+ tol=tol,
+ max.fail=max.fail,
+ pred.mean=(n==Nmif),
+ pred.var=(n==Nmif),
+ filter.mean=(n==Nmif),
+ mifiter=.ndone+n,
+ perturb.fn=perturb.fn,
+ transform=transform,
+ verbose=verbose,
+ .getnativesymbolinfo=gnsi
+ ),
+ silent=FALSE
+ )
+ if (inherits(pfp,"try-error"))
+ stop("mif2 particle-filter error")
+
+ gnsi <- FALSE
+
+ theta <- rowMeans(pfp at paramMatrix[,,drop=FALSE])
+
+ conv.rec[n+1,-c(1,2)] <- theta
+ conv.rec[n,c(1,2)] <- c(pfp at loglik,pfp at nfail)
+
+ if (verbose) cat("MIF2 iteration ",n," of ",Nmif," completed\n")
+
+ } ### end of main loop
+
+ ## back transform the parameter estimate if necessary
+ if (transform) theta <- partrans(pfp,theta,dir="forward")
+
+ new(
+ "mif2d.pomp",
+ pfp,
+ params=theta,
+ paramMatrix=pfp at paramMatrix,
+ transform=transform,
+ Nmif=Nmif,
+ perturb.fn=perturb.fn,
+ conv.rec=conv.rec,
+ tol=tol
+ )
+}
+
+setMethod(
+ "mif2",
+ signature=signature(object="pomp"),
+ function (object, Nmif = 1,
+ start, Np, perturb.fn,
+ tol = 1e-17, max.fail = Inf,
+ verbose = getOption("verbose"),
+ transform = FALSE,
+ ...) {
+
+ if (missing(start)) start <- coef(object)
+ if (length(start)==0)
+ stop(
+ "mif2 error: ",sQuote("start")," must be specified if ",
+ sQuote("coef(object)")," is NULL",
+ call.=FALSE
+ )
+
+ if (missing(Np))
+ stop("mif2 error: ",sQuote("Np")," must be specified",call.=FALSE)
+
+ if (missing(perturb.fn))
+ stop("mif2 error: ",sQuote("perturb.fn")," must be specified",call.=FALSE)
+ perturb.fn <- match.fun(perturb.fn)
+ if (!all(c('params','mifiter','timeno','...')%in%names(formals(perturb.fn)))) {
+ stop(
+ "mif2 error: ",
+ sQuote("perturb.fn"),
+ " must be a function of prototype ",
+ sQuote("perturb.fn(params,mifiter,timeno,...)"),
+ call.=FALSE
+ )
+ }
+
+ mif2.internal(
+ object=object,
+ Nmif=Nmif,
+ start=start,
+ Np=Np,
+ perturb.fn=perturb.fn,
+ tol=tol,
+ max.fail=max.fail,
+ transform=transform,
+ verbose=verbose,
+ ...
+ )
+
+ }
+ )
+
+
+setMethod(
+ "mif2",
+ signature=signature(object="pfilterd.pomp"),
+ function (object, Nmif = 1, Np, tol, ...) {
+
+ if (missing(Np)) Np <- object at Np
+ if (missing(tol)) tol <- object at tol
+
+ mif2(
+ object=as(object,"pomp"),
+ Nmif=Nmif,
+ Np=Np,
+ tol=tol,
+ ...
+ )
+ }
+ )
+
+setMethod(
+ "mif2",
+ signature=signature(object="mif2d.pomp"),
+ function (object, Nmif, start,
+ Np, perturb.fn,
+ tol, transform,
+ ...) {
+
+ if (missing(Nmif)) Nmif <- object at Nmif
+ if (missing(start)) start <- coef(object)
+ if (missing(perturb.fn)) perturb.fn <- object at perturb.fn
+ if (missing(transform)) transform <- object at transform
+
+ if (missing(Np)) Np <- object at Np
+ if (missing(tol)) tol <- object at tol
+
+ mif2(
+ object=as(object,"pomp"),
+ Nmif=Nmif,
+ start=start,
+ Np=Np,
+ perturb.fn=perturb.fn,
+ tol=tol,
+ transform=transform,
+ ...
+ )
+ }
+ )
+
+setMethod(
+ 'continue',
+ signature=signature(object='mif2d.pomp'),
+ function (object, Nmif = 1,
+ ...) {
+
+ ndone <- object at Nmif
+
+ obj <- mif2(
+ object=object,
+ Nmif=Nmif,
+ .ndone=ndone,
+ ...
+ )
+
+ object at conv.rec[ndone+1,c('loglik','nfail')] <- obj at conv.rec[1L,c('loglik','nfail')]
+ obj at conv.rec <- rbind(
+ object at conv.rec,
+ obj at conv.rec[-1L,colnames(object at conv.rec)]
+ )
+ obj at Nmif <- as.integer(ndone+Nmif)
+
+ obj
+ }
+ )
+
+default.cooling <- function (object, fraction, par.sd, ic.sd) {
+
+ nT <- length(time(object))
+ theta <- (1-fraction)/fraction/(50*nT-1)
+
+ function (params, mifiter, timeno, ...) {
+ pert <- params
+ sigma <- 1/(1+theta*((mifiter-1)*nT+timeno-1))
+ if (timeno==1) {
+ pert[names(ic.sd)] <- rnorm(
+ n=length(ic.sd),
+ mean=pert[names(ic.sd)],
+ sd=ic.sd*sigma
+ )
+ }
+ pert[names(par.sd)] <- rnorm(
+ n=length(par.sd),
+ mean=pert[names(par.sd)],
+ sd=par.sd*sigma
+ )
+ pert
+ }
+}
More information about the pomp-commits
mailing list