[Pomp-commits] r1178 - in pkg/pomp: . R tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jun 4 20:07:50 CEST 2015
Author: kingaa
Date: 2015-06-04 20:07:50 +0200 (Thu, 04 Jun 2015)
New Revision: 1178
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/NAMESPACE
pkg/pomp/R/mif.R
pkg/pomp/R/mif2.R
pkg/pomp/tests/ou2-mif2.R
pkg/pomp/tests/ou2-mif2.Rout.save
Log:
- refactor mif2 to separate perturbation rw.sd from cooling
- export mif2 methods
- rejigger the way that rw.sd are specified for mif2
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2015-06-04 13:06:28 UTC (rev 1177)
+++ pkg/pomp/DESCRIPTION 2015-06-04 18:07:50 UTC (rev 1178)
@@ -1,7 +1,7 @@
Package: pomp
Type: Package
Title: Statistical Inference for Partially Observed Markov Processes
-Version: 0.66-1
+Version: 0.66-2
Date: 2015-06-04
Authors at R: c(person(given=c("Aaron","A."),family="King",
role=c("aut","cre"),email="kingaa at umich.edu"),
Modified: pkg/pomp/NAMESPACE
===================================================================
--- pkg/pomp/NAMESPACE 2015-06-04 13:06:28 UTC (rev 1177)
+++ pkg/pomp/NAMESPACE 2015-06-04 18:07:50 UTC (rev 1178)
@@ -45,6 +45,7 @@
pomp,
pfilterd.pomp,
mif,mifList,
+ mif2d.pomp,mif2List,
pmcmc,pmcmcList,
traj.matched.pomp,
nlfd.pomp,
@@ -64,7 +65,8 @@
time,"time<-",timezero,"timezero<-",
simulate,pfilter,
eff.sample.size,cond.logLik,
- particles,mif,continue,states,trajectory,
+ states,trajectory,
+ particles,mif,mif2,continue,
pred.mean,pred.var,filter.mean,conv.rec,
values,
bsmc2,bsmc,pmcmc,abc,nlf,
@@ -89,6 +91,7 @@
onestep.dens,
gillespie.sim,
mvn.diag.rw,mvn.rw,
+ rw.sd,
sobol,
sobolDesign,
sliceDesign,
Modified: pkg/pomp/R/mif.R
===================================================================
--- pkg/pomp/R/mif.R 2015-06-04 13:06:28 UTC (rev 1177)
+++ pkg/pomp/R/mif.R 2015-06-04 18:07:50 UTC (rev 1178)
@@ -37,7 +37,7 @@
)
}
-mif1.cooling.function <- function (type, perobs, fraction, ntimes) {
+mif.cooling.function <- function (type, perobs, fraction, ntimes) {
switch(
type,
geometric={
@@ -55,13 +55,6 @@
}
},
hyperbolic={
- if (fraction>=1)
- stop(
- "mif error: ",sQuote("cooling.fraction.50"),
- " must be < 1 when cooling.type = ",
- sQuote("hyperbolic"),
- call.=FALSE
- )
if (perobs) {
scal <- (50*ntimes*fraction-1)/(1-fraction)
function (nt, m) {
@@ -93,27 +86,27 @@
paramMatrix = NULL,
.getnativesymbolinfo = TRUE,
...) {
-
+
pompLoad(object)
gnsi <- as.logical(.getnativesymbolinfo)
transform <- as.logical(transform)
-
+
if (length(start)==0)
stop(
"mif error: ",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("mif error: ",sQuote("start")," must be a named vector",call.=FALSE)
-
+
rw.names <- names(rw.sd)
if (is.null(rw.names) || any(rw.sd<0))
stop("mif error: ",sQuote("rw.sd")," must be a named non-negative numerical vector",call.=FALSE)
@@ -122,12 +115,12 @@
rw.names <- names(rw.sd[rw.sd>0])
if (length(rw.names) == 0)
stop("mif error: ",sQuote("rw.sd")," must have one positive entry for each parameter to be estimated",call.=FALSE)
-
+
if (is.null(pars))
pars <- rw.names[!(rw.names%in%ivps)]
else
warning("mif warning: argument ",sQuote("pars")," is redundant and deprecated. It will be removed in a future release.",call.=FALSE)
-
+
if (
!is.character(pars) ||
!is.character(ivps) ||
@@ -147,7 +140,7 @@
sQuote("rw.sd"),
call.=FALSE
)
-
+
if (!all(rw.names%in%c(pars,ivps))) {
extra.rws <- rw.names[!(rw.names%in%c(pars,ivps))]
warning(
@@ -164,7 +157,7 @@
}
rw.sd <- rw.sd[c(pars,ivps)]
rw.names <- names(rw.sd)
-
+
ntimes <- length(time(object))
if (is.null(Np)) stop("mif error: ",sQuote("Np")," must be specified",call.=FALSE)
if (is.function(Np)) {
@@ -184,7 +177,7 @@
if (!is.numeric(Np))
stop(sQuote("Np")," must be a number, a vector of numbers, or a function")
Np <- as.integer(Np)
-
+
ic.lag <- as.integer(ic.lag)
if ((length(ic.lag)!=1)||(ic.lag<1))
stop("mif error: ",sQuote("ic.lag")," must be a positive integer",call.=FALSE)
@@ -205,40 +198,40 @@
call.=FALSE
)
}
-
+
if (missing(cooling.fraction.50))
stop("mif error: ",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("mif error: ",sQuote("cooling.fraction.50")," must be a number between 0 and 1",call.=FALSE)
-
- cooling <- mif1.cooling.function(
- type=cooling.type,
- perobs=(method=="mif2"),
- fraction=cooling.fraction.50,
- ntimes=ntimes
- )
+ cooling <- mif.cooling.function(
+ type=cooling.type,
+ perobs=(method=="mif2"),
+ fraction=cooling.fraction.50,
+ ntimes=ntimes
+ )
+
if ((method=="mif2")&&(Np[1L]!=Np[ntimes+1]))
stop("the first and last values of ",sQuote("Np")," must agree when method = ",sQuote("mif2"))
-
+
if ((length(var.factor)!=1)||(var.factor < 0))
stop("mif error: ",sQuote("var.factor")," must be a positive number",call.=FALSE)
-
+
Nmif <- as.integer(Nmif)
if (Nmif<0)
stop("mif error: ",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,
@@ -249,7 +242,7 @@
)
)
conv.rec[1L,] <- c(NA,NA,theta)
-
+
if (!all(is.finite(theta[c(pars,ivps)]))) {
stop(
sQuote("mif")," cannot estimate non-finite parameters.\n",
@@ -261,15 +254,15 @@
call.=FALSE
)
}
-
+
obj <- as(object,"pomp")
-
+
if (Nmif>0) {
tmp.mif <- new("mif",object,particles=particles,Np=Np[1L])
} else {
pfp <- obj
}
-
+
have.parmat <- !(is.null(paramMatrix) || length(paramMatrix)==0)
for (n in seq_len(Nmif)) { ## iterate the filtering
@@ -277,7 +270,7 @@
## 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 <- try(
particles(
@@ -288,7 +281,7 @@
),
silent = FALSE
)
- if (inherits(P,"try-error"))
+ if (inherits(P,"try-error"))
stop("mif error: error in ",sQuote("particles"),call.=FALSE)
if ((method=="mif2") && ((n>1) || have.parmat)) {
@@ -299,7 +292,7 @@
pfp <- try(
pfilter.internal(
object=obj,
- params=P,
+ params=P,
Np=Np,
tol=tol,
max.fail=max.fail,
@@ -311,18 +304,18 @@
.mif2=(method=="mif2"),
.rw.sd=sigma.n[pars],
.transform=transform,
- save.states=FALSE,
+ save.states=FALSE,
save.params=FALSE,
verbose=verbose,
.getnativesymbolinfo=gnsi
),
silent=FALSE
)
- if (inherits(pfp,"try-error"))
+ if (inherits(pfp,"try-error"))
stop("mif error: error in ",sQuote("pfilter"),call.=FALSE)
gnsi <- FALSE
-
+
## update parameters
switch(
method,
@@ -344,14 +337,14 @@
theta[ivps] <- pfp at filter.mean[ivps,ic.lag]
conv.rec[n+1,-c(1,2)] <- theta
conv.rec[n,c(1,2)] <- c(pfp at loglik,pfp at 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)
new(
@@ -390,9 +383,9 @@
verbose = getOption("verbose"),
transform = FALSE,
...) {
-
+
method <- match.arg(method)
-
+
if (missing(start)) start <- coef(object)
if (missing(rw.sd))
stop("mif error: ",sQuote("rw.sd")," must be specified",call.=FALSE)
@@ -410,7 +403,7 @@
stop("mif error: ",sQuote("Np")," must be specified",call.=FALSE)
cooling.type <- match.arg(cooling.type)
-
+
if (missing(particles)) { # use default: normal distribution
particles <- default.mif.particles.fun
} else {
@@ -424,7 +417,7 @@
call.=FALSE
)
}
-
+
mif.internal(
object=object,
Nmif=Nmif,
@@ -444,7 +437,7 @@
transform=transform,
...
)
-
+
}
)
@@ -454,10 +447,10 @@
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
-
+
mif(
object=as(object,"pomp"),
Nmif=Nmif,
@@ -481,7 +474,7 @@
tol,
transform,
...) {
-
+
if (missing(Nmif)) Nmif <- object at Nmif
if (missing(start)) start <- coef(object)
if (missing(ivps)) ivps <- object at ivps
@@ -522,9 +515,9 @@
signature=signature(object='mif'),
function (object, Nmif = 1,
...) {
-
+
ndone <- object at Nmif
-
+
obj <- mif(
object=object,
Nmif=Nmif,
@@ -532,7 +525,7 @@
paramMatrix=object at paramMatrix,
...
)
-
+
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,
@@ -540,7 +533,7 @@
)
names(dimnames(obj at conv.rec)) <- c("iteration","variable")
obj at Nmif <- as.integer(ndone+Nmif)
-
+
obj
}
)
Modified: pkg/pomp/R/mif2.R
===================================================================
--- pkg/pomp/R/mif2.R 2015-06-04 13:06:28 UTC (rev 1177)
+++ pkg/pomp/R/mif2.R 2015-06-04 18:07:50 UTC (rev 1178)
@@ -1,143 +1,49 @@
## 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)
- )
- )
+rw.sd <- function (...) {
+ as.list(match.call())[-1L]
+}
+pkern.sd <- function (rw.sd, time, paramnames) {
+ if (!all(names(rw.sd) %in% paramnames)) {
+ unrec <- names(rw.sd)[!names(rw.sd) %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)
+ }
+ ivp <- function (sd, lag = 1L) {
+ sd*(seq_along(time)==lag)
+ }
+ sds <- lapply(rw.sd,eval,envir=list(time=time,ivp=ivp))
+ for (n in names(sds)) {
+ len <- length(sds[[n]])
+ if (len!=1 && len!=length(time))
+ stop(sQuote("mif2")," error: ",sQuote("rw.sd")," spec for parameter ",
+ sQuote(n)," does not evaluate to a vector of the correct length (",
+ length(time),")",call.=FALSE)
+ }
+ do.call(rbind,sds)
+}
+
## define the mif2d.pomp class
setClass(
'mif2d.pomp',
contains='pfilterd.pomp',
slots=c(
Nmif = 'integer',
+ rw.sd = 'list',
+ cooling.type = 'character',
+ cooling.fraction.50 = 'numeric',
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,
+ mifiter, rw.sd, cooling.fn,
tol = 1e-17, max.fail = Inf,
- transform = FALSE, verbose = FALSE,
- filter.mean = FALSE,
+ transform, verbose, filter.mean,
.getnativesymbolinfo = TRUE) {
gnsi <- as.logical(.getnativesymbolinfo)
@@ -152,6 +58,7 @@
paramnames <- rownames(params)
npars <- nrow(params)
+ rwnames <- rownames(rw.sd)
loglik <- rep(NA,ntimes)
eff.sample.size <- numeric(ntimes)
@@ -160,7 +67,9 @@
for (nt in seq_len(ntimes)) {
## perturb parameters
- params <- perturb.fn(params,sd=cooling.fn(mifiter,nt))
+ pmag <- cooling.fn(nt,mifiter)$alpha*rw.sd[,nt]
+ params[rwnames,] <- rnorm(n=length(rwnames)*ncol(params),
+ mean=params[rwnames,],sd=pmag)
if (transform)
tparams <- partrans(object,params,dir="fromEstimationScale",
@@ -280,29 +189,27 @@
)
}
-mif2.internal <- function (object, Nmif, start, Np, rw.sd, perturb.fn = NULL,
- tol = 1e17, max.fail = Inf, transform = FALSE,
+mif2.internal <- function (object, Nmif, start, Np, rw.sd, transform = FALSE,
+ cooling.type, cooling.fraction.50,
+ tol = 1e17, max.fail = Inf,
verbose = FALSE, .ndone = 0L,
.getnativesymbolinfo = TRUE, ...) {
pompLoad(object)
+ transform <- as.logical(transform)
+ verbose <- as.logical(verbose)
gnsi <- as.logical(.getnativesymbolinfo)
- Nmif <- as.integer(Nmif)
- if (Nmif<0) stop(sQuote("mif2")," error: ",sQuote("Nmif"),
- " must be a positive integer",call.=FALSE)
+ cooling.fn <- mif.cooling.function(
+ type=cooling.type,
+ perobs=TRUE,
+ fraction=cooling.fraction.50,
+ ntimes=length(time(object))
+ )
- ntimes <- length(time(object))
+ rw.sd.mat <- pkern.sd(rw.sd,time=time(object),paramnames=names(start))
- 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))))
@@ -333,7 +240,7 @@
Np=Np,
mifiter=.ndone+n,
cooling.fn=cooling.fn,
- perturb.fn=perturb.fn,
+ rw.sd=rw.sd.mat,
tol=tol,
max.fail=max.fail,
verbose=verbose,
@@ -367,10 +274,10 @@
pfp,
Nmif=Nmif,
rw.sd=rw.sd,
- perturb.fn=perturb.fn,
+ cooling.type=cooling.type,
+ cooling.fraction.50=cooling.fraction.50,
transform=transform,
- conv.rec=conv.rec,
- tol=tol
+ conv.rec=conv.rec
)
}
@@ -379,10 +286,17 @@
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,
+ definition = function (object, Nmif = 1, start, Np,
+ rw.sd, transform = FALSE,
+ cooling.type = c("hyperbolic", "geometric"),
+ cooling.fraction.50,
+ tol = 1e-17, max.fail = Inf,
verbose = getOption("verbose"),...) {
+ Nmif <- as.integer(Nmif)
+ if (Nmif<0) stop(sQuote("mif2")," error: ",sQuote("Nmif"),
+ " must be a positive integer",call.=FALSE)
+
if (missing(start)) start <- coef(object)
if (length(start)==0)
stop(
@@ -390,8 +304,13 @@
sQuote("coef(object)")," is NULL",
call.=FALSE
)
+ if (is.null(names(start)))
+ stop(sQuote("mif2")," error: ",sQuote("start")," must be a named vector",
+ 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)) {
@@ -418,34 +337,24 @@
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
- )
- }
- }
+ cooling.type <- match.arg(cooling.type)
+ cooling.fraction.50 <- as.numeric(cooling.fraction.50)
+ if (cooling.fraction.50 <= 0 || cooling.fraction.50 > 1)
+ stop(sQuote("mif2")," error: ",
+ sQuote("cooling.fraction.50")," must be in (0,1]",call.=FALSE)
+
mif2.internal(
object=object,
Nmif=Nmif,
start=start,
Np=Np,
rw.sd=rw.sd,
- perturb.fn=perturb.fn,
+ transform=transform,
+ cooling.type=cooling.type,
+ cooling.fraction.50=cooling.fraction.50,
tol=tol,
max.fail=max.fail,
- transform=transform,
verbose=verbose,
...
)
@@ -462,35 +371,33 @@
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,
- ...
- )
+ f <- selectMethod("mif2","pomp")
+ f(object=object,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, ...) {
+ definition = function (object, Nmif, start, Np,
+ rw.sd, transform, cooling.type, cooling.fraction.50,
+ tol, ...) {
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(cooling.type)) cooling.type <- object at cooling.type
+ if (missing(cooling.fraction.50)) cooling.fraction.50 <- object at cooling.fraction.50
+
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,...)
+ f(object,Nmif=Nmif,start=start,Np=Np,rw.sd=rw.sd,transform=transform,
+ cooling.type=cooling.type,cooling.fraction.50=cooling.fraction.50,
+ tol=tol,...)
}
)
@@ -611,7 +518,7 @@
xx <- z[[1]]
parnames <- names(coef(xx,transform=xx at transform))
estnames <- parnames
-
+
## plot filter means
filt.diag <- rbind("eff. sample size"=xx at eff.sample.size,filter.mean(xx))
filtnames <- rownames(filt.diag)
Modified: pkg/pomp/tests/ou2-mif2.R
===================================================================
--- pkg/pomp/tests/ou2-mif2.R 2015-06-04 13:06:28 UTC (rev 1177)
+++ pkg/pomp/tests/ou2-mif2.R 2015-06-04 18:07:50 UTC (rev 1178)
@@ -108,18 +108,16 @@
guess2[c('x1.0','x2.0','alpha.2','alpha.3')] <- 1.5*guess1[c('x1.0','x2.0','alpha.2','alpha.3')]
m1 <- pomp:::mif2(ou2,Nmif=100,start=guess1,Np=1000,
- rw.sd=pomp:::mif2.sd(
- x1.0=pomp:::ivphypcool(0.5,0.05),
- x2.0=pomp:::ivphypcool(0.5,0.05),
- alpha.2=pomp:::hypcool(0.1,0.05,ntimes=100),
- alpha.3=pomp:::hypcool(0.1,0.05,ntimes=100)))
+ cooling.type="hyperbolic",cooling.fraction.50=0.05,
+ rw.sd=pomp:::rw.sd(
+ x1.0=ivp(0.5),x2.0=ivp(0.5),
+ alpha.2=0.1,alpha.3=0.1))
m2 <- pomp:::mif2(ou2,Nmif=100,start=guess2,Np=1000,
- rw.sd=pomp:::mif2.sd(
- x1.0=pomp:::ivphypcool(0.5,0.05),
- x2.0=pomp:::ivphypcool(0.5,0.05),
- alpha.2=pomp:::hypcool(0.1,0.05,ntimes=100),
- alpha.3=pomp:::hypcool(0.1,0.05,ntimes=100)))
+ cooling.type="hyperbolic",cooling.fraction.50=0.05,
+ rw.sd=pomp:::rw.sd(
+ x1.0=ivp(0.5),x2.0=ivp(0.5),
+ alpha.2=0.1,alpha.3=0.1))
plot(c(m1,m2))
Modified: pkg/pomp/tests/ou2-mif2.Rout.save
===================================================================
--- pkg/pomp/tests/ou2-mif2.Rout.save 2015-06-04 13:06:28 UTC (rev 1177)
+++ pkg/pomp/tests/ou2-mif2.Rout.save 2015-06-04 18:07:50 UTC (rev 1178)
@@ -139,32 +139,30 @@
> guess2[c('x1.0','x2.0','alpha.2','alpha.3')] <- 1.5*guess1[c('x1.0','x2.0','alpha.2','alpha.3')]
>
> m1 <- pomp:::mif2(ou2,Nmif=100,start=guess1,Np=1000,
-+ rw.sd=pomp:::mif2.sd(
-+ x1.0=pomp:::ivphypcool(0.5,0.05),
-+ x2.0=pomp:::ivphypcool(0.5,0.05),
-+ alpha.2=pomp:::hypcool(0.1,0.05,ntimes=100),
-+ alpha.3=pomp:::hypcool(0.1,0.05,ntimes=100)))
++ cooling.type="hyperbolic",cooling.fraction.50=0.05,
++ rw.sd=pomp:::rw.sd(
++ x1.0=ivp(0.5),x2.0=ivp(0.5),
++ alpha.2=0.1,alpha.3=0.1))
>
> m2 <- pomp:::mif2(ou2,Nmif=100,start=guess2,Np=1000,
-+ rw.sd=pomp:::mif2.sd(
-+ x1.0=pomp:::ivphypcool(0.5,0.05),
-+ x2.0=pomp:::ivphypcool(0.5,0.05),
-+ alpha.2=pomp:::hypcool(0.1,0.05,ntimes=100),
-+ alpha.3=pomp:::hypcool(0.1,0.05,ntimes=100)))
++ cooling.type="hyperbolic",cooling.fraction.50=0.05,
++ rw.sd=pomp:::rw.sd(
++ x1.0=ivp(0.5),x2.0=ivp(0.5),
++ alpha.2=0.1,alpha.3=0.1))
>
> plot(c(m1,m2))
>
> rbind(mle1=c(coef(m1),loglik=logLik(pfilter(m1,Np=1000))),
+ mle2=c(coef(m2),loglik=logLik(pfilter(m1,Np=1000))),
+ truth=c(coef(ou2),loglik=logLik(pfilter(m1,Np=1000))))
- alpha.1 alpha.2 alpha.3 alpha.4 sigma.1 sigma.2 sigma.3 tau x1.0
-mle1 0.8 -0.4342775 0.341218 0.9 3 -0.5 2 1 -2.445886
-mle2 0.8 -0.5406991 0.326855 0.9 3 -0.5 2 1 -4.132440
-truth 0.8 -0.5000000 0.300000 0.9 3 -0.5 2 1 -3.000000
- x2.0 loglik
-mle1 1.897252 -481.8861
-mle2 2.828333 -485.8904
-truth 4.000000 -482.0197
+ alpha.1 alpha.2 alpha.3 alpha.4 sigma.1 sigma.2 sigma.3 tau
+mle1 0.8 -0.4489983 0.3184109 0.9 3 -0.5 2 1
+mle2 0.8 -0.5080060 0.2439865 0.9 3 -0.5 2 1
+truth 0.8 -0.5000000 0.3000000 0.9 3 -0.5 2 1
+ x1.0 x2.0 loglik
+mle1 -1.442418 1.585629 -481.6106
+mle2 -2.696841 3.157511 -482.3154
+truth -3.000000 4.000000 -481.9939
>
> dev.off()
null device
@@ -172,4 +170,4 @@
>
> proc.time()
user system elapsed
- 68.784 0.092 69.269
+ 70.828 0.064 71.333
More information about the pomp-commits
mailing list