[Pomp-commits] r836 - branches/mif2/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Mar 11 15:22:36 CET 2013
Author: nxdao2000
Date: 2013-03-11 15:22:35 +0100 (Mon, 11 Mar 2013)
New Revision: 836
Modified:
branches/mif2/R/mif.R
Log:
rw.sd is either a vector, a list or a matrix will be converted to a matrix rwsdMat in mif
Modified: branches/mif2/R/mif.R
===================================================================
--- branches/mif2/R/mif.R 2013-03-11 14:20:26 UTC (rev 835)
+++ branches/mif2/R/mif.R 2013-03-11 14:22:35 UTC (rev 836)
@@ -2,55 +2,55 @@
default.pomp.particles.fun <- function (Np, center, sd, ...) {
matrix(
- data=rnorm(
- n=Np*length(center),
- mean=center,
- sd=sd
- ),
- nrow=length(center),
- ncol=Np,
- dimnames=list(
- names(center),
- NULL
- )
- )
+ data=rnorm(
+ n=Np*length(center),
+ mean=center,
+ sd=sd
+ ),
+ nrow=length(center),
+ ncol=Np,
+ dimnames=list(
+ names(center),
+ NULL
+ )
+ )
}
cooling.function <- function (type, perobs, fraction, ntimes) {
switch(
- type,
- geometric={
- factor <- fraction^(1/50)
- if (perobs) {
- function (nt, m) {
- alpha <- factor^(nt/ntimes+m-1)
- list(alpha=alpha,gamma=alpha^2)
- }
- } else {
- function (nt, m) {
- alpha <- factor^(m-1)
- list(alpha=alpha,gamma=alpha^2)
- }
- }
- },
- hyperbolic={
- if (perobs) {
- scal <- (50*ntimes*fraction-1)/(1-fraction)
- function (nt, m) {
- alpha <- (1+scal)/(scal+nt+ntimes*(m-1))
- list(alpha=alpha,gamma=alpha^2)
- }
- } else {
- scal <- (50*fraction-1)/(1-fraction)
- function (nt, m) {
- alpha <- (1+scal)/(scal+m-1)
- list(alpha=alpha,gamma=alpha^2)
- }
-
- }
- },
- stop("unrecognized cooling schedule type ",sQuote(type))
- )
+ type,
+ geometric={
+ factor <- fraction^(1/50)
+ if (perobs) {
+ function (nt, m) {
+ alpha <- factor^(nt/ntimes+m-1)
+ list(alpha=alpha,gamma=alpha^2)
+ }
+ } else {
+ function (nt, m) {
+ alpha <- factor^(m-1)
+ list(alpha=alpha,gamma=alpha^2)
+ }
+ }
+ },
+ hyperbolic={
+ if (perobs) {
+ scal <- (50*ntimes*fraction-1)/(1-fraction)
+ function (nt, m) {
+ alpha <- (1+scal)/(scal+nt+ntimes*(m-1))
+ list(alpha=alpha,gamma=alpha^2)
+ }
+ } else {
+ scal <- (50*fraction-1)/(1-fraction)
+ function (nt, m) {
+ alpha <- (1+scal)/(scal+m-1)
+ list(alpha=alpha,gamma=alpha^2)
+ }
+
+ }
+ },
+ stop("unrecognized cooling schedule type ",sQuote(type))
+ )
}
mif.cooling <- function (factor, n) { # default geometric cooling schedule
@@ -90,15 +90,15 @@
.getnativesymbolinfo = TRUE) {
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
- )
+ "mif error: ",sQuote("start")," must be specified if ",
+ sQuote("coef(object)")," is NULL",
+ call.=FALSE
+ )
if (transform)
start <- partrans(object,start,dir="inverse")
@@ -107,7 +107,7 @@
if (is.null(start.names))
stop("mif error: ",sQuote("start")," must be a named vector",call.=FALSE)
- rw.names <- names(rw.sd)
+ rw.names <- colnames(rw.sd)
if (is.null(rw.names))
stop("mif error: ",sQuote("rw.sd")," must be a named non-negative numerical vector",call.=FALSE)
if (!all(rw.names%in%start.names))
@@ -120,7 +120,7 @@
if (missing(ivps)) stop("mif error: ",sQuote("ivps")," must be specified",call.=FALSE)
if (
- !is.character(pars) ||
+ !is.character(pars) ||
!is.character(ivps) ||
!all(pars%in%start.names) ||
!all(ivps%in%start.names) ||
@@ -128,34 +128,34 @@
any(ivps%in%pars) ||
!all(pars%in%rw.names) ||
!all(ivps%in%rw.names)
- )
+ )
stop(
- "mif error: ",
- sQuote("pars")," and ",sQuote("ivps"),
- " must be mutually disjoint subsets of ",
- sQuote("names(start)"),
- " and must have a positive random-walk SDs specified in ",
- sQuote("rw.sd"),
- call.=FALSE
- )
+ "mif error: ",
+ sQuote("pars")," and ",sQuote("ivps"),
+ " must be mutually disjoint subsets of ",
+ sQuote("names(start)"),
+ " and must have a positive random-walk SDs specified in ",
+ sQuote("rw.sd"),
+ call.=FALSE
+ )
if (!all(rw.names%in%c(pars,ivps))) {
extra.rws <- rw.names[!(rw.names%in%c(pars,ivps))]
warning(
- ngettext(length(extra.rws),"mif warning: the variable ",
- "mif warning: the variables "),
- paste(sQuote(extra.rws),collapse=", "),
- ngettext(length(extra.rws)," has positive random-walk SD specified, but is included in neither ",
- " have positive random-walk SDs specified, but are included in neither "),
- sQuote("pars")," nor ",sQuote("ivps"),
- ngettext(length(extra.rws),". This random walk SD will be ignored.",
- ". These random walk SDs will be ignored."),
- call.=FALSE
- )
+ ngettext(length(extra.rws),"mif warning: the variable ",
+ "mif warning: the variables "),
+ paste(sQuote(extra.rws),collapse=", "),
+ ngettext(length(extra.rws)," has positive random-walk SD specified, but is included in neither ",
+ " have positive random-walk SDs specified, but are included in neither "),
+ sQuote("pars")," nor ",sQuote("ivps"),
+ ngettext(length(extra.rws),". This random walk SD will be ignored.",
+ ". These random walk SDs will be ignored."),
+ call.=FALSE
+ )
}
- rw.sd <- rw.sd[c(pars,ivps)]
- rw.names <- names(rw.sd)
-
+ #rw.sd <- rw.sd[c(pars,ivps)]
+ rw.names <- colnames(rw.sd)
+ rwsdMat <-rw.sd
if (missing(particles))
stop("mif error: ",sQuote("particles")," must be specified",call.=FALSE)
@@ -163,9 +163,9 @@
if (missing(Np)) stop("mif error: ",sQuote("Np")," must be specified",call.=FALSE)
if (is.function(Np)) {
Np <- try(
- vapply(seq.int(from=0,to=ntimes,by=1),Np,numeric(1)),
- silent=FALSE
- )
+ vapply(seq.int(from=0,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")
}
@@ -185,20 +185,20 @@
stop("mif error: ",sQuote("ic.lag")," must be a positive integer",call.=FALSE)
if (ic.lag>ntimes) {
warning(
- "mif warning: ",sQuote("ic.lag")," = ",ic.lag," > ",ntimes,
- " = length(time(",sQuote("object"),"))",
- " is nonsensical. Setting ",sQuote("ic.lag")," = ",ntimes,".",
- call.=FALSE
- )
+ "mif warning: ",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(
- "mif warning: 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
- )
+ "mif warning: 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
+ )
}
## the following deals with the deprecated option 'cooling.factor'
@@ -217,7 +217,7 @@
call.=FALSE)
}
}
-
+
if (missing(cooling.fraction))
stop("mif error: ",sQuote("cooling.fraction")," must be specified",call.=FALSE)
cooling.fraction <- as.numeric(cooling.fraction)
@@ -225,12 +225,12 @@
stop("mif error: ",sQuote("cooling.fraction")," must be a number between 0 and 1",call.=FALSE)
cooling <- cooling.function(
- type=cooling.type,
- perobs=(method=="mif2"),
- fraction=cooling.fraction,
- ntimes=ntimes
- )
-
+ type=cooling.type,
+ perobs=(method=="mif2"),
+ fraction=cooling.fraction,
+ 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"))
@@ -240,88 +240,31 @@
Nmif <- as.integer(Nmif)
if (Nmif<0)
stop("mif error: ",sQuote("Nmif")," must be a positive integer",call.=FALSE)
-
+
theta <- start
- dtheta <- length(start)
- rwsdMat <- data.frame(matrix(rep(0, dtheta*(ntimes+1)), ncol=dtheta, nrow = (ntimes+1)))
- names(rwsdMat) <- names(start)
-
- sigma <- rep(0,length(start))
- names(sigma) <- start.names
-
- rw.names <- c(pars,ivps)
-
- #sigma[rw.names] <- rw.sd
-
- for (i in 1:length(rw.names))
- { if (rw.names[i] %in% ivps)
- {
- if (length((rw.sd[[rw.names[i]]])==1))
- {
- rwsdMat[1,rw.names[i]] <- rw.sd[[rw.names[i]]][1]
- sigma[rw.names[i]] <- rw.sd[[rw.names[i]]][1]
- }
- else if (length(rw.sd[[rw.names[i]]])==(ntimes+1))
- {
- sigma[rw.names[i]] <- rw.sd[[rw.names[i]]][1]
-
- for (j in 1:(ntimes+1) )
- rwsdMat[j,rw.names[i]] <- rw.sd[[rw.names[i]]][j]
- }
- else
- {
- stop(sQuote("rw.sd")," must have length 1 or length ",ntimes+1)
-
- }
- }
- else if (rw.names[i] %in% pars)
- {
- if (length(rw.sd[[rw.names[i]]])==1)
- {
- sigma[rw.names[i]] <- rw.sd[[rw.names[i]]]
- #rwsdMat[1,rw.names[i]] <- rw.sd[rw.names[i]]
- for (j in 1:(ntimes+1) )
- rwsdMat[j,rw.names[i]] <- rw.sd[[rw.names[i]]]
- }
- else if (length(rw.sd[[rw.names[i]]])==(ntimes+1))
- {
- sigma[rw.names[i]] <- rw.sd[[rw.names[i]]][1]
-
- for (j in 1:(ntimes+1) )
- rwsdMat[j,rw.names[i]] <- rw.sd[[rw.names[i]]][j]
- }
- else
- {
- stop(sQuote("rw.sd")," must have length 1 or length ",ntimes+1)
-
- }
- }
-
- }
-
conv.rec <- matrix(
- data=NA,
- nrow=Nmif+1,
- ncol=length(theta)+2,
- dimnames=list(
- seq(.ndone,.ndone+Nmif),
- c('loglik','nfail',names(theta))
- )
- )
+ data=NA,
+ nrow=Nmif+1,
+ ncol=length(theta)+2,
+ dimnames=list(
+ seq(.ndone,.ndone+Nmif),
+ c('loglik','nfail',names(theta))
+ )
+ )
conv.rec[1L,] <- c(NA,NA,theta)
if (!all(is.finite(theta[c(pars,ivps)]))) {
stop(
- 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
- )
+ 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
+ )
}
obj <- as(object,"pomp")
@@ -333,282 +276,353 @@
}
have.parmat <- !(missing(paramMatrix) || length(paramMatrix)==0)
-
+
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 <- as.numeric(rwsdMat[1,])*cool.sched$alpha
- names(sigma.n)<-names(start)
+ sigma.n <- rwsdMat[1,]*cool.sched$alpha
+
## initialize the parameter portions of the particles
P <- try(
- particles(
- tmp.mif,
- Np=Np[1L],
- center=theta,
- sd=sigma.n*var.factor
- ),
- silent = FALSE
- )
+ particles(
+ tmp.mif,
+ Np=Np[1L],
+ center=theta,
+ sd=sigma.n*var.factor
+ ),
+ silent = FALSE
+ )
if (inherits(P,"try-error"))
stop("mif error: error in ",sQuote("particles"),call.=FALSE)
-
+
if ((method=="mif2") && ((n>1) || have.parmat)) {
## use pre-existing particle matrix
P[pars,] <- paramMatrix[pars,]
}
- names(rwsdMat) <- names(start)
+ colnames(rwsdMat) <- names(start)
pfp <- try(
- pfilter.internal(
- object=obj,
- params=P,
- Np=Np,
- tol=tol,
- max.fail=max.fail,
- pred.mean=(n==Nmif),
- pred.var=((method=="mif")||(n==Nmif)),
- filter.mean=TRUE,
- cooling=cooling,
- cooling.m=.ndone+n,
- .mif2=(method=="mif2"),
- .rw.sd=rwsdMat,
- .transform=transform,
- save.states=FALSE,
- save.params=FALSE,
- verbose=verbose,
- .getnativesymbolinfo=gnsi
- ),
- silent=FALSE
- )
+ pfilter.internal(
+ object=obj,
+ params=P,
+ Np=Np,
+ tol=tol,
+ max.fail=max.fail,
+ pred.mean=(n==Nmif),
+ pred.var=((method=="mif")||(n==Nmif)),
+ filter.mean=TRUE,
+ cooling=cooling,
+ cooling.m=.ndone+n,
+ .mif2=(method=="mif2"),
+ .rw.sd=rwsdMat*cool.sched$alpha,
+ .transform=transform,
+ save.states=FALSE,
+ save.params=FALSE,
+ verbose=verbose,
+ .getnativesymbolinfo=gnsi
+ ),
+ silent=FALSE
+ )
if (inherits(pfp,"try-error"))
stop("mif error: error in ",sQuote("pfilter"),call.=FALSE)
-
+
gnsi <- FALSE
## update parameters
switch(
- method,
- mif={ # original Ionides et al. (2006) average
- theta <- .Call(mif_update,pfp,theta,cool.sched$gamma,var.factor,sigma,pars)
- },
- unweighted={ # unweighted average
- theta[pars] <- rowMeans(pfp at filter.mean[pars,,drop=FALSE])
- },
- fp={ # fixed-point iteration
- theta[pars] <- pfp at filter.mean[pars,ntimes,drop=FALSE]
- },
- mif2={ # "efficient" iterated filtering
- paramMatrix <- pfp at paramMatrix
- theta[pars] <- rowMeans(paramMatrix[pars,,drop=FALSE])
- },
- stop("unrecognized method ",sQuote(method))
- )
- theta[ivps] <- pfp at filter.mean[ivps,ic.lag]
+ method,
+ mif={ # original Ionides et al. (2006) average
+ theta <- .Call(mif_update,pfp,theta,cool.sched$gamma,var.factor,rwsdMat[1,],pars)
+ theta[ivps] <- pfp at filter.mean[ivps,ic.lag]
+ },
+ unweighted={ # unweighted average
+ theta[pars] <- rowMeans(pfp at filter.mean[pars,,drop=FALSE])
+ theta[ivps] <- pfp at filter.mean[ivps,ic.lag]
+ },
+ fp={ # fixed-point iteration
+ theta[pars] <- pfp at filter.mean[pars,ntimes,drop=FALSE]
+ theta[ivps] <- pfp at filter.mean[ivps,ic.lag]
+ },
+ mif2={ # "efficient" iterated filtering
+ paramMatrix <- pfp at paramMatrix
+ theta[pars] <- rowMeans(paramMatrix[pars,,drop=FALSE])
+ theta[ivps] <- pfp at filter.mean[ivps,ntimes]
+ },
+ stop("unrecognized method ",sQuote(method))
+ )
+
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="forward")
new(
- "mif",
- pfp,
- transform=transform,
- params=theta,
- ivps=ivps,
- pars=pars,
- Nmif=Nmif,
- particles=particles,
- 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=cooling.fraction,
- paramMatrix=if (method=="mif2") paramMatrix else array(data=numeric(0),dim=c(0,0))
- )
+ "mif",
+ pfp,
+ transform=transform,
+ params=theta,
+ ivps=ivps,
+ pars=pars,
+ Nmif=Nmif,
+ particles=particles,
+ var.factor=var.factor,
+ ic.lag=ic.lag,
+ random.walk.sd=rwsdMat,
+ tol=tol,
+ conv.rec=conv.rec,
+ method=method,
+ cooling.type=cooling.type,
+ cooling.fraction=cooling.fraction,
+ paramMatrix=if (method=="mif2") paramMatrix else array(data=numeric(0),dim=c(0,0))
+ )
}
setGeneric('mif',function(object,...)standardGeneric("mif"))
setMethod(
- "mif",
- signature=signature(object="pomp"),
- function (object, Nmif = 1,
- start,
- pars, ivps = character(0),
- particles, rw.sd,
- Np, ic.lag, var.factor,
- cooling.type = c("geometric","hyperbolic"),
- cooling.fraction, cooling.factor,
- method = c("mif","unweighted","fp","mif2"),
- tol = 1e-17, max.fail = Inf,
- verbose = getOption("verbose"),
- transform = FALSE,
- ...) {
+ "mif",
+ signature=signature(object="pomp"),
+ function (object, Nmif = 1,
+ start,
+ pars, ivps = character(0),
+ particles, rw.sd,
+ Np, ic.lag, var.factor,
+ cooling.type = c("geometric","hyperbolic"),
+ cooling.fraction, cooling.factor,
+ method = c("mif","unweighted","fp","mif2"),
+ tol = 1e-17, max.fail = Inf,
+ verbose = getOption("verbose"),
+ transform = FALSE,
+ ...) {
+
+ transform <- as.logical(transform)
+ method <- match.arg(method)
+
+
+ ntimes <- length(time(object))
+ if (missing(start)) start <- coef(object)
+ if (missing(rw.sd))
+ stop("mif error: ",sQuote("rw.sd")," must be specified",call.=FALSE)
+ if (missing(pars)) {
+ stop("mif error: ",sQuote("par")," must be specified",call.=FALSE)
+ }
+ rw.names <- c(pars,ivps)
+ dtheta <- length(start)
+ rwsdMat <- matrix(rep(0, dtheta*(ntimes+1)), ncol=dtheta, nrow = (ntimes+1))
+ colnames(rwsdMat) <- names(start)
+
+ if (is.matrix(rw.sd)) rwsdMat<-rw.sd
+ else if (is.list(rw.sd))
+ {
+ for (i in 1:length(rw.names))
+ { if (rw.names[i] %in% ivps)
+ {
+ if (length((rw.sd[[rw.names[i]]])==1))
+ {
+ rwsdMat[1,rw.names[i]] <- rw.sd[[rw.names[i]]][1]
+ }
+ else
+ {
+ stop(sQuote("rw.sd")," must have length 1 or length ",ntimes+1)
+
+ }
+ }
+ else if (rw.names[i] %in% pars)
+ {
+ if (length(rw.sd[[rw.names[i]]])==1)
+ {
+ rwsdMat[,rw.names[i]] <- rep(rw.sd[[rw.names[i]]],(ntimes+1))
+ }
+ else if (length(rw.sd[[rw.names[i]]])==(ntimes+1))
+ {
+ rwsdMat[,rw.names[i]] <- rw.sd[[rw.names[i]]]
+ }
+ else
+ {
+ stop(sQuote("rw.sd")," must have length 1 or length ",ntimes+1)
- transform <- as.logical(transform)
- 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)
- if (missing(pars)) {
- rw.names <- names(rw.sd)[rw.sd>0]
- pars <- rw.names[!(rw.names%in%ivps)]
- }
- if (missing(Np))
- stop("mif error: ",sQuote("Np")," must be specified",call.=FALSE)
- if (missing(ic.lag) && length(ivps)>0)
- stop("mif error: ",sQuote("ic.lag")," must be specified if ",sQuote("ivps")," are",call.=FALSE)
- if (missing(var.factor))
- stop("mif error: ",sQuote("var.factor")," must be specified",call.=FALSE)
-
- cooling.type <- match.arg(cooling.type)
-
- if (missing(particles)) { # use default: normal distribution
- particles <- default.pomp.particles.fun
- } else {
- particles <- match.fun(particles)
- if (!all(c('Np','center','sd','...')%in%names(formals(particles))))
- stop(
- "mif error: ",
- sQuote("particles"),
- " must be a function of prototype ",
- sQuote("particles(Np,center,sd,...)"),
- call.=FALSE
- )
- }
-
- mif.internal(
- object=object,
- Nmif=Nmif,
- start=start,
- pars=pars,
- ivps=ivps,
- particles=particles,
- rw.sd=rw.sd,
- Np=Np,
- cooling.type=cooling.type,
- cooling.factor=cooling.factor,
- cooling.fraction=cooling.fraction,
- var.factor=var.factor,
- ic.lag=ic.lag,
- method=method,
- tol=tol,
- max.fail=max.fail,
- verbose=verbose,
- transform=transform
- )
-
}
- )
+ }
+
+ }
+
+ }
+ else if (is.vector(rw.sd))
+ {
+ if (missing(pars)) {
+ rw.names <- names(rw.sd)[rw.sd>0]
+ pars <- rw.names[!(rw.names%in%ivps)]
+ }
+ for (i in 1:length(rw.names))
+ { if (rw.names[i] %in% ivps)
+ {
+ rwsdMat[1,rw.names[i]] <- rw.sd[rw.names[i]]
+
+ }
+ else if (rw.names[i] %in% pars)
+ {
+ rwsdMat[,rw.names[i]] <- rep(rw.sd[rw.names[i]],(ntimes+1))
+ }
+
+ }
+
+
+
+ }
+
+
+
+ if (missing(Np))
+ stop("mif error: ",sQuote("Np")," must be specified",call.=FALSE)
+ if (missing(ic.lag) && length(ivps)>0)
+ stop("mif error: ",sQuote("ic.lag")," must be specified if ",sQuote("ivps")," are",call.=FALSE)
+ if (missing(var.factor))
+ stop("mif error: ",sQuote("var.factor")," must be specified",call.=FALSE)
+
+ cooling.type <- match.arg(cooling.type)
+
+ if (missing(particles)) { # use default: normal distribution
+ particles <- default.pomp.particles.fun
+ } else {
+ particles <- match.fun(particles)
+ if (!all(c('Np','center','sd','...')%in%names(formals(particles))))
+ stop(
+ "mif error: ",
+ sQuote("particles"),
+ " must be a function of prototype ",
+ sQuote("particles(Np,center,sd,...)"),
+ call.=FALSE
+ )
+ }
+
+ mif.internal(
+ object=object,
+ Nmif=Nmif,
+ start=start,
+ pars=pars,
+ ivps=ivps,
+ particles=particles,
+ rw.sd=rwsdMat,
+ Np=Np,
+ cooling.type=cooling.type,
+ cooling.factor=cooling.factor,
+ cooling.fraction=cooling.fraction,
+ var.factor=var.factor,
+ ic.lag=ic.lag,
+ method=method,
+ tol=tol,
+ max.fail=max.fail,
+ verbose=verbose,
+ transform=transform
+ )
+
+ }
+)
setMethod(
- "mif",
- 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,
- Np=Np,
- tol=tol,
- ...
- )
- }
- )
+ "mif",
+ 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,
+ Np=Np,
+ tol=tol,
+ ...
+ )
+ }
+)
setMethod(
- "mif",
- signature=signature(object="mif"),
- function (object, Nmif,
- start,
- pars, ivps,
- particles, rw.sd,
- Np, ic.lag, var.factor,
- cooling.type, cooling.fraction,
- method,
- tol,
- transform,
- ...) {
-
- if (missing(Nmif)) Nmif <- object at Nmif
- if (missing(start)) start <- coef(object)
- if (missing(pars)) pars <- object at pars
- if (missing(ivps)) ivps <- object at ivps
- if (missing(particles)) particles <- object at particles
- if (missing(rw.sd)) rw.sd <- object at random.walk.sd
- if (missing(ic.lag)) ic.lag <- object at ic.lag
- if (missing(var.factor)) var.factor <- object at var.factor
- if (missing(cooling.type)) cooling.type <- object at cooling.type
- if (missing(cooling.fraction)) cooling.fraction <- object at cooling.fraction
- if (missing(method)) method <- object at method
- if (missing(transform)) transform <- object at transform
- transform <- as.logical(transform)
+ "mif",
+ signature=signature(object="mif"),
+ function (object, Nmif,
+ start,
+ pars, ivps,
+ particles, rw.sd,
+ Np, ic.lag, var.factor,
+ cooling.type, cooling.fraction,
+ method,
+ tol,
+ transform,
+ ...) {
+
+ if (missing(Nmif)) Nmif <- object at Nmif
+ if (missing(start)) start <- coef(object)
+ if (missing(pars)) pars <- object at pars
+ if (missing(ivps)) ivps <- object at ivps
+ if (missing(particles)) particles <- object at particles
+ if (missing(rw.sd)) rw.sd <- object at random.walk.sd
+ if (missing(ic.lag)) ic.lag <- object at ic.lag
+ if (missing(var.factor)) var.factor <- object at var.factor
+ if (missing(cooling.type)) cooling.type <- object at cooling.type
+ if (missing(cooling.fraction)) cooling.fraction <- object at cooling.fraction
+ if (missing(method)) method <- object at method
+ if (missing(transform)) transform <- object at transform
+ transform <- as.logical(transform)
+
+ if (missing(Np)) Np <- object at Np
+ if (missing(tol)) tol <- object at tol
+
+ mif(
+ object=as(object,"pomp"),
+ Nmif=Nmif,
+ start=start,
+ pars=pars,
+ ivps=ivps,
+ particles=particles,
+ rw.sd=rw.sd,
+ Np=Np,
+ cooling.type=cooling.type,
+ cooling.fraction=cooling.fraction,
+ var.factor=var.factor,
+ ic.lag=ic.lag,
+ method=method,
+ tol=tol,
+ transform=transform,
+ ...
+ )
+ }
+)
- if (missing(Np)) Np <- object at Np
- if (missing(tol)) tol <- object at tol
-
- mif(
- object=as(object,"pomp"),
- Nmif=Nmif,
- start=start,
- pars=pars,
- ivps=ivps,
- particles=particles,
- rw.sd=rw.sd,
- Np=Np,
- cooling.type=cooling.type,
- cooling.fraction=cooling.fraction,
- var.factor=var.factor,
- ic.lag=ic.lag,
- method=method,
- tol=tol,
- transform=transform,
- ...
- )
- }
- )
-
setMethod(
- 'continue',
- signature=signature(object='mif'),
- function (object, Nmif = 1,
- ...) {
-
- ndone <- object at Nmif
-
- obj <- mif(
- object=object,
- Nmif=Nmif,
- .ndone=ndone,
- 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,
- obj at conv.rec[-1L,colnames(object at conv.rec)]
- )
- obj at Nmif <- as.integer(ndone+Nmif)
-
- obj
- }
- )
+ 'continue',
+ signature=signature(object='mif'),
+ function (object, Nmif = 1,
+ ...) {
+
+ ndone <- object at Nmif
+
+ obj <- mif(
+ object=object,
+ Nmif=Nmif,
+ .ndone=ndone,
+ 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,
+ obj at conv.rec[-1L,colnames(object at conv.rec)]
+ )
+ obj at Nmif <- as.integer(ndone+Nmif)
+
+ obj
+ }
+)
mif.profileDesign <- function (object, profile, lower, upper, nprof, ivps,
rw.sd, Np, ic.lag, var.factor, cooling.factor,option, cooling.fraction, paramMatrix, ...)
@@ -633,24 +647,24 @@
ans <- vector(mode="list",length=nrow(pd))
for (k in seq_len(nrow(pd))) {
ans[[k]] <- list(
- mf=mif(
- object,
- Nmif=0,
- start=unlist(pd[k,]),
- pars=pars,
- ivps=ivps,
- rw.sd=rw.sd,
- Np=Np,
- ic.lag=ic.lag,
- var.factor=var.factor,
- cooling.factor=cooling.factor,
- option=option,
- cooling.fraction=cooling.fraction,
- paramMatrix=paramMatrix,
- ...
- )
- )
+ mf=mif(
+ object,
+ Nmif=0,
+ start=unlist(pd[k,]),
+ pars=pars,
+ ivps=ivps,
+ rw.sd=rw.sd,
+ Np=Np,
+ ic.lag=ic.lag,
+ var.factor=var.factor,
+ cooling.factor=cooling.factor,
+ option=option,
+ cooling.fraction=cooling.fraction,
+ paramMatrix=paramMatrix,
+ ...
+ )
+ )
}
ans
-}
+}
\ No newline at end of file
More information about the pomp-commits
mailing list