[Pomp-commits] r1176 - in pkg: . pomp pomp/R pomp/tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jun 4 15:06:26 CEST 2015
Author: kingaa
Date: 2015-06-04 15:06:26 +0200 (Thu, 04 Jun 2015)
New Revision: 1176
Added:
pkg/pomp/R/mif2.R
Removed:
pkg/mif2.R
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/tests/ou2-mif2.R
Log:
- mif2 integrated into pomp, but not yet exported
Deleted: pkg/mif2.R
===================================================================
--- pkg/mif2.R 2015-06-04 12:22:37 UTC (rev 1175)
+++ pkg/mif2.R 2015-06-04 13:06:26 UTC (rev 1176)
@@ -1,841 +0,0 @@
-require(pomp)
-
-## MIF2 algorithm functions
-
-## define a class of perturbation magnitudes
-setClass(
- "mif2.perturb.sd",
- slots=c(
- sds="list",
- rwnames="character"
- ),
- prototype=prototype(
- sds=list(),
- rwnames=character(0)
- )
-)
-
-## define the mif2d.pomp class
-setClass(
- 'mif2d.pomp',
- contains='pfilterd.pomp',
- slots=c(
- Nmif = 'integer',
- transform = 'logical',
- perturb.fn = 'function',
- rw.sd = 'mif2.perturb.sd',
- conv.rec = 'matrix'
- )
-)
-
-mif2.sd <- function (...) {
- 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)
- for (n in rwnames) {
- sds[[n]] <- try(match.fun(sds[[n]]),silent=TRUE)
- if (inherits(sds[[n]],"try-error"))
- stop(sQuote("mif2.sd")," error: ",sQuote(n)," is not a function",call.=FALSE)
- }
- new("mif2.perturb.sd",sds=sds,rwnames=rwnames)
-}
-
-setGeneric("cooling_fn",function(object,...)standardGeneric("cooling_fn"))
-
-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)
- }
- function (mifiter, timept) {
- rw.sd <- setNames(numeric(length(object at rwnames)),object at rwnames)
- for (nm in object at rwnames) {
- rw.sd[nm] <- object at sds[[nm]](mifiter,timept)
- }
- rw.sd
- }
- })
-
-geomcool <- function (sd, cooling.fraction.50 = 1) {
- if (missing(sd))
- stop(sQuote("geomcool")," error: ",sQuote("sd")," must be supplied",call.=FALSE)
- sd <- as.numeric(sd)
- if (sd <= 0)
- stop(sQuote("geomcool")," error: ",sQuote("sd")," must be non-negative",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)
- factor <- cooling.fraction.50^(1/50)
- function (mifiter, timept) {
- sd*(factor^(mifiter-1))
- }
-}
-
-hypcool <- function (sd, cooling.fraction.50 = 1, ntimes = NULL) {
- if (missing(sd))
- stop(sQuote("hypcool")," error: ",sQuote("sd")," must be supplied",call.=FALSE)
- sd <- as.numeric(sd)
- if (sd <= 0)
- stop(sQuote("hypcool")," error: ",sQuote("sd")," must be non-negative",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)
- if (is.null(ntimes)) {
- scal <- (50*cooling.fraction.50-1)/(1-cooling.fraction.50)
- function (mifiter, timept) {
- (1+scal)/(mifiter+scal)
- }
- } else {
- scal <- (50*ntimes*cooling.fraction.50-1)/(1-cooling.fraction.50)
- function (mifiter, timept) {
- sd*(1+scal)/(timept+ntimes*(mifiter-1)+scal)
- }
- }
-}
-
-ivphypcool <- function (sd, cooling.fraction.50 = 1) {
- if (missing(sd))
- stop(sQuote("ivphypcool")," error: ",sQuote("sd")," must be supplied",call.=FALSE)
- sd <- as.numeric(sd)
- if (sd <= 0)
- stop(sQuote("ivphypcool")," error: ",sQuote("sd")," must be non-negative",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)
- scal <- (50*cooling.fraction.50-1)/(1-cooling.fraction.50)
- function (mifiter, timept) {
- if (timept==1L)
- sd*(1+scal)/(mifiter+scal)
- else 0.0
- }
-}
-
-ivpgeomcool <- function (sd, cooling.fraction.50 = 1) {
- if (missing(sd))
- stop(sQuote("ivpgeomcool")," error: ",sQuote("sd")," must be supplied",call.=FALSE)
- sd <- as.numeric(sd)
- if (sd <= 0)
- stop(sQuote("ivpgeomcool")," error: ",sQuote("sd")," must be non-negative",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)
- factor <- cooling.fraction.50^(1/50)
- function (mifiter, timept) {
- if (timept==1L)
- sd*factor^(mifiter-1)
- else 0.0
- }
-}
-
-mif2.pfilter <- function (object, params, Np,
- mifiter, cooling.fn, perturb.fn,
- tol = 1e-17, max.fail = Inf,
- transform = FALSE, verbose = FALSE,
- filter.mean = FALSE,
- .getnativesymbolinfo = TRUE) {
-
- gnsi <- as.logical(.getnativesymbolinfo)
- transform <- as.logical(transform)
- verbose <- as.logical(verbose)
- filter.mean <- as.logical(filter.mean)
- mifiter <- as.integer(mifiter)
- Np <- as.integer(Np)
-
- times <- time(object,t0=TRUE)
- ntimes <- length(times)-1
-
- paramnames <- rownames(params)
- npars <- nrow(params)
-
- loglik <- rep(NA,ntimes)
- eff.sample.size <- numeric(ntimes)
- nfail <- 0
-
- for (nt in seq_len(ntimes)) {
-
- ## perturb parameters
- params <- perturb.fn(params,sd=cooling.fn(mifiter,nt))
-
- if (transform)
- tparams <- partrans(object,params,dir="fromEstimationScale",
- .getnativesymbolinfo=gnsi)
-
- if (nt == 1L) {
- ## get initial states
- x <- init.state(object,params=if (transform) tparams else params)
-
- if (filter.mean)
- filt.m <- array(dim=c(nrow(x),ntimes),
- dimnames=list(rownames(x),NULL))
- else
- filt.m <- array(dim=c(0,0))
- }
-
- ## 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
- ),
- silent=FALSE
- )
- if (inherits(X,'try-error'))
- stop(sQuote("mif2.pfilter")," error: process simulation error")
-
- ## 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
- ),
- silent=FALSE
- )
- if (inherits(weights,'try-error'))
- stop(sQuote("mif2.pfilter")," error: error in calculation of weights",call.=FALSE)
- if (any(!is.finite(weights))) {
- stop(sQuote("mif2.pfilter")," error: ",sQuote("dmeasure"),
- " returns non-finite value",call.=FALSE)
- }
- gnsi <- FALSE
-
- ## compute weighted mean at last timestep
- if (nt == ntimes)
- coef(object,transform=transform) <- apply(params,1L,weighted.mean,w=weights)
-
- ## compute effective sample size, log-likelihood
- ## also do resampling if filtering has not failed
- xx <- try(
- .Call(
- pomp:::pfilter_computations,
- x=X,
- params=params,
- Np=Np[nt+1],
- rw=FALSE,
- rw_sd=numeric(0),
- predmean=FALSE,
- predvar=FALSE,
- filtmean=filter.mean,
- onepar=FALSE,
- weights=weights,
- tol=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 (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,
- paramMatrix=params,
- eff.sample.size=eff.sample.size,
- cond.loglik=loglik,
- filter.mean=filt.m,
- Np=Np,
- tol=tol,
- nfail=as.integer(nfail),
- loglik=sum(loglik)
- )
-}
-
-mif2.internal <- function (object, Nmif, start, Np, rw.sd, perturb.fn = NULL,
- tol = 1e17, max.fail = Inf, transform = FALSE,
- verbose = FALSE, .ndone = 0L,
- .getnativesymbolinfo = TRUE, ...) {
-
- pompLoad(object)
-
- gnsi <- as.logical(.getnativesymbolinfo)
-
- Nmif <- as.integer(Nmif)
- if (Nmif<0) stop(sQuote("mif2")," error: ",sQuote("Nmif"),
- " must be a positive integer",call.=FALSE)
-
- ntimes <- length(time(object))
-
- Np <- as.integer(Np)
-
- if (is.null(names(start)))
- stop(sQuote("mif2")," error: ",sQuote("start")," must be a named vector",
- call.=FALSE)
-
- cooling.fn <- cooling_fn(rw.sd,paramnames=names(start))
-
- conv.rec <- array(data=NA,dim=c(Nmif+1,length(start)+2),
- dimnames=list(seq.int(.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(dim=c(0,0))
- }
-
- object <- as(object,"pomp")
-
- if (transform)
- paramMatrix <- partrans(object,paramMatrix,dir="toEstimationScale",
- .getnativesymbolinfo=gnsi)
-
- ## iterate the filtering
- for (n in seq_len(Nmif)) {
-
- pfp <- try(
- mif2.pfilter(
- object=object,
- params=paramMatrix,
- Np=Np,
- mifiter=.ndone+n,
- cooling.fn=cooling.fn,
- perturb.fn=perturb.fn,
- tol=tol,
- max.fail=max.fail,
- verbose=verbose,
- filter.mean=(n==Nmif),
- transform=transform,
- .getnativesymbolinfo=gnsi
- ),
- silent=FALSE
- )
- if (inherits(pfp,"try-error"))
- stop("mif2 particle-filter error")
-
- gnsi <- FALSE
-
- paramMatrix <- pfp at paramMatrix
- conv.rec[n+1,-c(1,2)] <- coef(pfp)
- conv.rec[n,c(1,2)] <- c(pfp at loglik,pfp at nfail)
-
- if (verbose) cat("mif2 iteration ",n," of ",Nmif," completed\n")
-
- }
-
- if (transform)
- pfp at paramMatrix <- partrans(object,paramMatrix,dir="fromEstimationScale",
- .getnativesymbolinfo=gnsi)
-
- pompUnload(object)
-
- new(
- "mif2d.pomp",
- pfp,
- Nmif=Nmif,
- rw.sd=rw.sd,
- perturb.fn=perturb.fn,
- transform=transform,
- conv.rec=conv.rec,
- tol=tol
- )
-}
-
-setGeneric("mif2",function(object,...)standardGeneric("mif2"))
-
-setMethod(
- "mif2",
- signature=signature(object="pomp"),
- definition = function (object, Nmif = 1, start, Np, rw.sd, perturb.fn,
- tol = 1e-17, max.fail = Inf, transform = FALSE,
- verbose = getOption("verbose"),...) {
-
- if (missing(start)) start <- coef(object)
- if (length(start)==0)
- stop(
- sQuote("mif2")," error: ",sQuote("start")," must be specified if ",
- sQuote("coef(object)")," is NULL",
- call.=FALSE
- )
-
- ntimes <- length(time(object))
- if (missing(Np)) {
- stop(sQuote("mif2")," error: ",sQuote("Np")," must be specified",call.=FALSE) }
- else if (is.function(Np)) {
- Np <- try(
- vapply(seq.int(1,ntimes),Np,numeric(1)),
- silent=FALSE
- )
- if (inherits(Np,"try-error"))
- stop(sQuote("mif2")," error: if ",sQuote("Np"),
- " is a function, it must return a single positive integer")
- } else if (!is.numeric(Np))
- stop(sQuote("mif2")," error: ",sQuote("Np"),
- " must be a number, a vector of numbers, or a function")
- if (length(Np)==1) {
- Np <- rep(Np,times=ntimes+1)
- } else if (length(Np)==ntimes) {
- Np <- c(Np,Np[1L])
- } else if (length(Np)>ntimes) {
- if (Np[1L] != Np[ntimes+1])
- stop(sQuote("mif2")," error: Np[ntimes+1] != Np[1]")
- if (length(Np)>ntimes+1)
- warning("in ",sQuote("mif2"),": Np[k] ignored for k > ntimes")
- }
- if (any(Np<=0))
- stop("number of particles, ",sQuote("Np"),", must always be positive")
-
- if (missing(perturb.fn)) {
- perturb.fn <- function (theta, sd) {
- theta[names(sd),] <- rnorm(n=length(sd)*ncol(theta),mean=theta[names(sd),],sd=sd)
- theta
- }
- } else {
- perturb.fn <- match.fun(perturb.fn)
- if (!all(c('theta','sd')%in%names(formals(perturb.fn)))) {
- stop(
- sQuote("mif2")," error: ",
- sQuote("perturb.fn"),
- " must be a function of prototype ",
- sQuote("perturb.fn(theta,sd)"),
- call.=FALSE
- )
- }
- }
-
- mif2.internal(
- object=object,
- Nmif=Nmif,
- start=start,
- Np=Np,
- rw.sd=rw.sd,
- perturb.fn=perturb.fn,
- tol=tol,
- max.fail=max.fail,
- transform=transform,
- verbose=verbose,
- ...
- )
-
- }
-)
-
-
-setMethod(
- "mif2",
- signature=signature(object="pfilterd.pomp"),
- definition = 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"),
- definition = function (object, Nmif, start, Np, rw.sd, perturb.fn, tol,
- transform, ...) {
-
- if (missing(Nmif)) Nmif <- object at Nmif
- if (missing(start)) start <- coef(object)
- if (missing(rw.sd)) rw.sd <- object at rw.sd
- 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
-
- f <- selectMethod("mif2","pomp")
-
- f(object,Nmif=Nmif,start=start,Np=Np,rw.sd=rw.sd,
- perturb.fn=perturb.fn,tol=tol,transform=transform,...)
- }
-)
-
-setMethod(
- 'continue',
- signature=signature(object='mif2d.pomp'),
- definition = 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
- }
-)
-
-## extract the estimated log likelihood
-setMethod('logLik','mif2d.pomp',function(object,...)object at loglik)
-
-setMethod('conv.rec','mif2d.pomp',
- function (object, pars, transform = FALSE, ...) {
- pomp:::conv.rec.internal(object=object,pars=pars,transform=transform,...)
- }
-)
-
-## mif2List class
-setClass(
- 'mif2List',
- contains='list',
- validity=function (object) {
- if (!all(sapply(object,is,'mif2d.pomp'))) {
- retval <- paste0(
- "error in ",sQuote("c"),
- ": dissimilar objects cannot be combined"
- )
- return(retval)
- }
- d <- sapply(object,function(x)dim(x at conv.rec))
- if (!all(apply(d,1,diff)==0)) {
- retval <- paste0(
- "error in ",sQuote("c"),
- ": to be combined, ",sQuote("mif2d.pomp"),
- " objects must equal numbers of iterations"
- )
- return(retval)
- }
- TRUE
- }
-)
-
-setMethod(
- 'c',
- signature=signature(x='mif2d.pomp'),
- definition=function (x, ...) {
- y <- list(...)
- if (length(y)==0) {
- new("mif2List",list(x))
- } else {
- p <- sapply(y,is,'mif2d.pomp')
- pl <- sapply(y,is,'mif2List')
- if (any(!(p||pl)))
- stop("cannot mix ",sQuote("mif2d.pomp"),
- " and non-",sQuote("mif2d.pomp")," objects")
- y[p] <- lapply(y[p],list)
- y[pl] <- lapply(y[pl],as,"list")
- new("mif2List",c(list(x),y,recursive=TRUE))
- }
- }
-)
-
-setMethod(
- 'c',
- signature=signature(x='mif2List'),
- definition=function (x, ...) {
- y <- list(...)
- if (length(y)==0) {
- x
- } else {
- p <- sapply(y,is,'mif2d.pomp')
- pl <- sapply(y,is,'mif2List')
- if (any(!(p||pl)))
- stop("cannot mix ",sQuote("mif2d.pomp"),
- " and non-",sQuote("mif2d.pomp")," objects")
- y[p] <- lapply(y[p],list)
- y[pl] <- lapply(y[pl],as,"list")
- new("mif2List",c(as(x,"list"),y,recursive=TRUE))
- }
- }
-)
-
-setMethod(
- "[",
- signature=signature(x="mif2List"),
- definition=function(x, i, ...) {
- new('mif2List',as(x,"list")[i])
- }
-)
-
-setMethod(
- 'conv.rec',
- signature=signature(object='mif2List'),
- definition=function (object, ...) {
- lapply(object,conv.rec,...)
- }
-)
-
-mif2.diagnostics <- function (z) {
- ## assumes that z is a list of mif2d.pomps with identical structure
- mar.multi <- c(0,5.1,0,2.1)
- oma.multi <- c(6,0,5,0)
- xx <- z[[1]]
- ivpnames <- xx at ivps
- estnames <- c(xx at pars,ivpnames)
- parnames <- names(coef(xx,transform=xx at transform))
- unestnames <- parnames[-match(estnames,parnames)]
-
- ## plot filter means
- filt.diag <- rbind("eff. sample size"=xx at eff.sample.size,filter.mean(xx))
- filtnames <- rownames(filt.diag)
- plotnames <- if(length(unestnames)>0) filtnames[-match(unestnames,filtnames)] else filtnames
- lognames <- filtnames[1] # eff. sample size
- nplots <- length(plotnames)
- n.per.page <- min(nplots,10)
- if(n.per.page<=4) nc <- 1 else nc <- 2
- nr <- ceiling(n.per.page/nc)
- oldpar <- par(mar=mar.multi,oma=oma.multi,mfcol=c(nr,nc),ask=dev.interactive(orNone=TRUE))
- on.exit(par(oldpar))
- low <- 1
- hi <- 0
- time <- time(xx)
- while (hi<nplots) {
- hi <- min(low+n.per.page-1,nplots)
- for (i in seq(from=low,to=hi,by=1)) {
- n <- i-low+1
- logplot <- if (plotnames[i]%in%lognames) "y" else ""
- dat <- sapply(
- z,
- function(po, label) {
- if (label=="eff. sample size")
- po at eff.sample.size
- else
- filter.mean(po,label)
- },
- label=plotnames[i]
- )
- matplot(
- y=dat,
- x=time,
- axes = FALSE,
- xlab = "",
- log=logplot,
- ylab = "",
- type = "l"
- )
- box()
- y.side <- 2
- axis(y.side, xpd = NA)
- mtext(plotnames[i], y.side, line = 3)
- do.xax <- (n%%nr==0||n==n.per.page)
- if (do.xax) axis(1,xpd=NA)
- if (do.xax) mtext("time",side=1,line=3)
- }
- low <- hi+1
- mtext("Filter diagnostics (last iteration)",3,line=2,outer=TRUE)
- }
-
- ## plot mif convergence diagnostics
- other.diagnostics <- c("loglik", "nfail")
- plotnames <- c(other.diagnostics,estnames)
- nplots <- length(plotnames)
- n.per.page <- min(nplots,10)
- nc <- if (n.per.page<=4) 1 else 2
- nr <- ceiling(n.per.page/nc)
- par(mar=mar.multi,oma=oma.multi,mfcol=c(nr,nc))
- ## on.exit(par(oldpar))
- low <- 1
- hi <- 0
- iteration <- seq(0,xx at Nmif)
- while (hi<nplots) {
- hi <- min(low+n.per.page-1,nplots)
- for (i in seq(from=low,to=hi,by=1)) {
- n <- i-low+1
- dat <- sapply(z,function(po,label) conv.rec(po,label),label=plotnames[i])
- matplot(
- y=dat,
- x=iteration,
- axes = FALSE,
- xlab = "",
- ylab = "",
- type = "l"
- )
- box()
- y.side <- 2
- axis(y.side,xpd=NA)
- mtext(plotnames[i],y.side,line=3)
- do.xax <- (n%%nr==0||n==n.per.page)
- if (do.xax) axis(1,xpd=NA)
- if (do.xax) mtext("MIF iteration",side=1,line=3)
- }
- low <- hi+1
- mtext("MIF convergence diagnostics",3,line=2,outer=TRUE)
- }
- invisible(NULL)
-}
-
-
-setMethod(
- "plot",
- "mif2d.pomp",
- function (x, y, ...) {
- if (!missing(y)) {
- y <- substitute(y)
- warning(sQuote(y)," is ignored")
- }
- pomp:::mif.diagnostics(list(x))
- }
- )
-
-setMethod(
- "plot",
- signature=signature(x='mif2List'),
- definition=function (x, y, ...) {
- if (!missing(y)) {
- y <- substitute(y)
- warning(sQuote(y)," is ignored")
- }
- pomp:::mif.diagnostics(x)
- }
- )
-
-
-
-require(ggplot2)
-require(plyr)
-require(reshape2)
-require(magrittr)
-theme_set(theme_bw())
-require(pomp)
-stopifnot(packageVersion("pomp")>="0.65-1")
-options(
- keep.source=TRUE,
- stringsAsFactors=FALSE,
- encoding="UTF-8",
- scipen=5,
- cores=5
-)
-
-## ----gompertz-init,cache=FALSE-------------------------------------------
-require(pomp)
-pompExample(gompertz)
-theta <- coef(gompertz)
-theta.true <- theta
-
-## ----gompertz-multi-mif2-eval,results='hide'-----------------------------
-require(foreach)
-require(doMC)
-registerDoMC()
-
-save.seed <- .Random.seed
-set.seed(334388458L,kind="L'Ecuyer")
-
-estpars <- c("r","sigma","tau")
-mf <- foreach(i=1:5,
- .inorder=FALSE,
- .options.multicore=list(set.seed=TRUE)
-) %dopar%
-{
- theta.guess <- theta.true
- theta.guess[estpars] <- rlnorm(
- n=length(estpars),
- meanlog=log(theta.guess[estpars]),
- sdlog=1
- )
- m1 <- mif2(
- gompertz,
- Nmif=50,
- start=theta.guess,
- transform=TRUE,
- rw.sd=mif2.sd(
- r=hypcool(0.02,0.95,ntimes=101),
- sigma=hypcool(0.02,0.95,ntimes=101),
- tau=hypcool(0.05,0.95,ntimes=101)
- ),
- Np=2000
- )
- m1 <- continue(m1,Nmif=50,rw.sd=mif2.sd(
- r=hypcool(0.02,0.8,ntimes=101),
- sigma=hypcool(0.02,0.8,ntimes=101),
- tau=hypcool(0.05,0.8,ntimes=101)
- ))
- m1 <- continue(m1,Nmif=50,rw.sd=mif2.sd(
- r=hypcool(0.02,0.6,ntimes=101),
- sigma=hypcool(0.02,0.6,ntimes=101),
- tau=hypcool(0.05,0.6,ntimes=101)
- ))
- m1 <- continue(m1,Nmif=50,rw.sd=mif2.sd(
- r=hypcool(0.02,0.2,ntimes=101),
- sigma=hypcool(0.02,0.2,ntimes=101),
- tau=hypcool(0.05,0.2,ntimes=101)
- ))
- ll <- replicate(n=10,logLik(pfilter(m1,Np=10000)))
- list(mif=m1,ll=ll)
-}
-
-## ----gompertz-post-mif2--------------------------------------------------
-loglik.true <- replicate(n=10,logLik(pfilter(gompertz,Np=10000)))
-loglik.true <- logmeanexp(loglik.true,se=TRUE)
-theta.mif <- t(sapply(mf,function(x)coef(x$mif)))
-loglik.mif <- t(sapply(mf,function(x)logmeanexp(x$ll,se=TRUE)))
-best <- which.max(loglik.mif[,1])
-theta.mif <- theta.mif[best,]
-loglik.mif <- loglik.mif[best,]
-rbind(
- mle=c(signif(theta.mif[estpars],3),loglik=round(loglik.mif,2)),
- truth=c(signif(theta.true[estpars],3),loglik=round(loglik.true,2))
-) -> results.table
-
-.Random.seed <<- save.seed
-
-## ----mif2-plot,echo=FALSE,cache=FALSE,fig.height=6-----------------------
-op <- par(mfrow=c(4,1),mar=c(3,3,0,0),mgp=c(2,1,0),bty='l')
-loglik <- sapply(mf,function(x)conv.rec(x$mif,"loglik"))
-log.r <- sapply(mf,function(x)log(conv.rec(x$mif,"r")))
-log.sigma <- sapply(mf,function(x)log(conv.rec(x$mif,"sigma")))
-log.tau <- sapply(mf,function(x)log(conv.rec(x$mif,"tau")))
-matplot(loglik,type='l',lty=1,xlab="",ylab=expression(log~L),xaxt='n',ylim=max(loglik,na.rm=T)+c(-12,3))
-matplot(log.r,type='l',lty=1,xlab="",ylab=expression(log~r),xaxt='n')
-abline(h=log(theta.true["r"]),lty=2)
-matplot(log.sigma,type='l',lty=1,xlab="",ylab=expression(log~sigma),xaxt='n')
-abline(h=log(theta.true["sigma"]),lty=2)
-matplot(log.tau,type='l',lty=1,xlab="MIF iteration",ylab=expression(log~tau))
-abline(h=log(theta.true["tau"]),lty=2)
-par(op)
-
-## ----first-mif-results-table,echo=FALSE,cache=FALSE----------------------
-print(results.table)
-
-# plot(do.call(c,lapply(mf,getElement,"mif")))
-
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2015-06-04 12:22:37 UTC (rev 1175)
+++ pkg/pomp/DESCRIPTION 2015-06-04 13:06:26 UTC (rev 1176)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical Inference for Partially Observed Markov Processes
-Version: 0.65-3
-Date: 2015-06-02
+Version: 0.66-1
+Date: 2015-06-04
Authors at R: c(person(given=c("Aaron","A."),family="King",
role=c("aut","cre"),email="kingaa at umich.edu"),
person(given=c("Edward","L."),family="Ionides",role=c("aut")),
@@ -34,7 +34,7 @@
simulate-pomp.R trajectory-pomp.R plot-pomp.R
pfilter.R pfilter-methods.R minim.R traj-match.R
bsmc.R bsmc2.R
- mif.R mif-methods.R
+ mif.R mif-methods.R mif2.R
proposals.R pmcmc.R pmcmc-methods.R
nlf-funcs.R nlf-guts.R nlf-objfun.R nlf.R
probe.R probe-match.R basic-probes.R spect.R spect-match.R
Copied: pkg/pomp/R/mif2.R (from rev 1175, pkg/mif2.R)
===================================================================
--- pkg/pomp/R/mif2.R (rev 0)
+++ pkg/pomp/R/mif2.R 2015-06-04 13:06:26 UTC (rev 1176)
@@ -0,0 +1,727 @@
+## MIF2 algorithm functions
+
+## define a class of perturbation magnitudes
+setClass(
+ "mif2.perturb.sd",
+ slots=c(
+ sds="list",
+ rwnames="character"
+ ),
+ prototype=prototype(
+ sds=list(),
+ rwnames=character(0)
+ )
+ )
+
+## define the mif2d.pomp class
+setClass(
+ 'mif2d.pomp',
+ contains='pfilterd.pomp',
+ slots=c(
+ Nmif = 'integer',
+ transform = 'logical',
+ perturb.fn = 'function',
+ rw.sd = 'mif2.perturb.sd',
+ conv.rec = 'matrix'
+ )
+ )
+
+mif2.sd <- function (...) {
+ 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)
+ for (n in rwnames) {
+ sds[[n]] <- try(match.fun(sds[[n]]),silent=TRUE)
+ if (inherits(sds[[n]],"try-error"))
+ stop(sQuote("mif2.sd")," error: ",sQuote(n)," is not a function",call.=FALSE)
+ }
+ new("mif2.perturb.sd",sds=sds,rwnames=rwnames)
+}
+
+setGeneric("cooling_fn",function(object,...)standardGeneric("cooling_fn"))
+
+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)
+ }
+ function (mifiter, timept) {
+ rw.sd <- setNames(numeric(length(object at rwnames)),object at rwnames)
+ for (nm in object at rwnames) {
+ rw.sd[nm] <- object at sds[[nm]](mifiter,timept)
+ }
+ rw.sd
+ }
+ })
+
+geomcool <- function (sd, cooling.fraction.50 = 1) {
+ if (missing(sd))
+ stop(sQuote("geomcool")," error: ",sQuote("sd")," must be supplied",call.=FALSE)
+ sd <- as.numeric(sd)
+ if (sd <= 0)
+ stop(sQuote("geomcool")," error: ",sQuote("sd")," must be non-negative",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)
+ factor <- cooling.fraction.50^(1/50)
+ function (mifiter, timept) {
+ sd*(factor^(mifiter-1))
+ }
+}
+
+hypcool <- function (sd, cooling.fraction.50 = 1, ntimes = NULL) {
+ if (missing(sd))
+ stop(sQuote("hypcool")," error: ",sQuote("sd")," must be supplied",call.=FALSE)
+ sd <- as.numeric(sd)
+ if (sd <= 0)
+ stop(sQuote("hypcool")," error: ",sQuote("sd")," must be non-negative",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)
+ if (is.null(ntimes)) {
+ scal <- (50*cooling.fraction.50-1)/(1-cooling.fraction.50)
+ function (mifiter, timept) {
+ (1+scal)/(mifiter+scal)
+ }
+ } else {
+ scal <- (50*ntimes*cooling.fraction.50-1)/(1-cooling.fraction.50)
+ function (mifiter, timept) {
+ sd*(1+scal)/(timept+ntimes*(mifiter-1)+scal)
+ }
+ }
+}
+
+ivphypcool <- function (sd, cooling.fraction.50 = 1) {
+ if (missing(sd))
+ stop(sQuote("ivphypcool")," error: ",sQuote("sd")," must be supplied",call.=FALSE)
+ sd <- as.numeric(sd)
+ if (sd <= 0)
+ stop(sQuote("ivphypcool")," error: ",sQuote("sd")," must be non-negative",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)
+ scal <- (50*cooling.fraction.50-1)/(1-cooling.fraction.50)
+ function (mifiter, timept) {
+ if (timept==1L)
+ sd*(1+scal)/(mifiter+scal)
+ else 0.0
+ }
+}
+
+ivpgeomcool <- function (sd, cooling.fraction.50 = 1) {
+ if (missing(sd))
+ stop(sQuote("ivpgeomcool")," error: ",sQuote("sd")," must be supplied",call.=FALSE)
+ sd <- as.numeric(sd)
+ if (sd <= 0)
+ stop(sQuote("ivpgeomcool")," error: ",sQuote("sd")," must be non-negative",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)
+ factor <- cooling.fraction.50^(1/50)
+ function (mifiter, timept) {
+ if (timept==1L)
+ sd*factor^(mifiter-1)
+ else 0.0
+ }
+}
+
+mif2.pfilter <- function (object, params, Np,
+ mifiter, cooling.fn, perturb.fn,
+ tol = 1e-17, max.fail = Inf,
+ transform = FALSE, verbose = FALSE,
+ filter.mean = FALSE,
+ .getnativesymbolinfo = TRUE) {
+
+ gnsi <- as.logical(.getnativesymbolinfo)
+ transform <- as.logical(transform)
+ verbose <- as.logical(verbose)
+ filter.mean <- as.logical(filter.mean)
+ mifiter <- as.integer(mifiter)
+ Np <- as.integer(Np)
+
+ times <- time(object,t0=TRUE)
+ ntimes <- length(times)-1
+
+ paramnames <- rownames(params)
+ npars <- nrow(params)
+
+ loglik <- rep(NA,ntimes)
+ eff.sample.size <- numeric(ntimes)
+ nfail <- 0
+
+ for (nt in seq_len(ntimes)) {
+
+ ## perturb parameters
+ params <- perturb.fn(params,sd=cooling.fn(mifiter,nt))
+
+ if (transform)
+ tparams <- partrans(object,params,dir="fromEstimationScale",
+ .getnativesymbolinfo=gnsi)
+
+ if (nt == 1L) {
+ ## get initial states
+ x <- init.state(object,params=if (transform) tparams else params)
+
+ if (filter.mean)
+ filt.m <- array(dim=c(nrow(x),ntimes),
+ dimnames=list(rownames(x),NULL))
+ else
+ filt.m <- array(dim=c(0,0))
+ }
+
+ ## 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
+ ),
+ silent=FALSE
+ )
+ if (inherits(X,'try-error'))
+ stop(sQuote("mif2.pfilter")," error: process simulation error")
+
+ ## 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
+ ),
+ silent=FALSE
+ )
+ if (inherits(weights,'try-error'))
+ stop(sQuote("mif2.pfilter")," error: error in calculation of weights",call.=FALSE)
+ if (any(!is.finite(weights))) {
+ stop(sQuote("mif2.pfilter")," error: ",sQuote("dmeasure"),
+ " returns non-finite value",call.=FALSE)
+ }
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 1176
More information about the pomp-commits
mailing list