[Pomp-commits] r179 - in pkg: R tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Dec 9 13:16:05 CET 2009
Author: kingaa
Date: 2009-12-09 13:16:05 +0100 (Wed, 09 Dec 2009)
New Revision: 179
Modified:
pkg/R/mif.R
pkg/tests/ou2-mif.R
Log:
- bug fix in 'mif' relating to non-standard particles functions
Modified: pkg/R/mif.R
===================================================================
--- pkg/R/mif.R 2009-11-16 11:18:55 UTC (rev 178)
+++ pkg/R/mif.R 2009-12-09 12:16:05 UTC (rev 179)
@@ -23,6 +23,7 @@
Np = NULL, cooling.factor = NULL, var.factor = NULL, ic.lag = NULL,
weighted = TRUE, tol = 1e-17, warn = TRUE, max.fail = 0,
verbose = FALSE, .ndone = 0) {
+ is.mif <- is(object,"mif")
if (is.null(start)) {
start <- coef(object)
if (length(start)==0)
@@ -35,8 +36,12 @@
if (is.null(start.names))
stop("mif error: ",sQuote("start")," must be a named vector",call.=FALSE)
- if (is.null(rw.sd))
- rw.sd <- object at random.walk.sd
+ if (is.null(rw.sd)) {
+ if (is.mif)
+ rw.sd <- object at random.walk.sd
+ else
+ stop("mif error: ",sQuote("rw.sd")," must be specified",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)
@@ -47,12 +52,19 @@
stop("mif error: ",sQuote("rw.sd")," must have one positive entry for each parameter to be estimated",call.=FALSE)
if (is.null(pars)) {
- pars <- object at pars
+ if (is.mif)
+ pars <- object at pars
+ else
+ stop("mif error: ",sQuote("pars")," must be specified",call.=FALSE)
}
if (length(pars)==0)
stop("mif error: at least one ordinary (non-IVP) parameter must be estimated",call.=FALSE)
- if (is.null(ivps))
- ivps <- object at ivps
+ if (is.null(ivps)) {
+ if (is.mif)
+ ivps <- object at ivps
+ else
+ stop("mif error: ",sQuote("ivps")," must be specified",call.=FALSE)
+ }
if (
!is.character(pars) ||
!is.character(ivps) ||
@@ -87,25 +99,45 @@
rw.sd <- rw.sd[c(pars,ivps)]
rw.names <- names(rw.sd)
- if (is.null(particles))
- particles <- object at particles
+ if (is.null(particles)) {
+ if (is.mif)
+ particles <- object at particles
+ else
+ stop("mif error: ",sQuote("particles")," must be specified",call.=FALSE)
+ }
- if (is.null(Np))
- Np <- object at alg.pars$Np
+ if (is.null(Np)) {
+ if (is.mif)
+ Np <- object at alg.pars$Np
+ else
+ stop("mif error: ",sQuote("Np")," must be specified",call.=FALSE)
+ }
Np <- as.integer(Np)
if ((length(Np)!=1)||(Np < 1))
stop("mif error: ",sQuote("Np")," must be a positive integer",call.=FALSE)
- if (is.null(ic.lag))
- ic.lag <- object at alg.pars$ic.lag
+ if (is.null(ic.lag)) {
+ if (is.mif)
+ ic.lag <- object at alg.pars$ic.lag
+ else
+ stop("mif error: ",sQuote("ic.lag")," must be specified",call.=FALSE)
+ }
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)
- if (is.null(cooling.factor))
- cooling.factor <- object at alg.pars$cooling.factor
+ if (is.null(cooling.factor)) {
+ if (is.mif)
+ cooling.factor <- object at alg.pars$cooling.factor
+ else
+ stop("mif error: ",sQuote("cooling.factor")," must be specified",call.=FALSE)
+ }
if ((length(cooling.factor)!=1)||(cooling.factor < 0)||(cooling.factor>1))
stop("mif error: ",sQuote("cooling.factor")," must be a number between 0 and 1",call.=FALSE)
- if (is.null(var.factor))
- var.factor <- object at alg.pars$var.factor
+ if (is.null(var.factor)) {
+ if (is.mif)
+ var.factor <- object at alg.pars$var.factor
+ else
+ stop("mif error: ",sQuote("var.factor")," must be specified",call.=FALSE)
+ }
if ((length(var.factor)!=1)||(var.factor < 0))
stop("mif error: ",sQuote("var.factor")," must be a positive number",call.=FALSE)
@@ -160,6 +192,23 @@
)
}
+ obj <- new(
+ "mif",
+ as(object,"pomp"),
+ ivps=ivps,
+ pars=pars,
+ Nmif=0L,
+ particles=particles,
+ alg.pars=list(
+ Np=Np,
+ cooling.factor=cooling.factor,
+ var.factor=var.factor,
+ ic.lag=ic.lag
+ ),
+ random.walk.sd=sigma[rw.names],
+ conv.rec=conv.rec
+ )
+
for (n in seq(length=Nmif)) { # main loop
## compute the cooled sigma
@@ -173,7 +222,7 @@
## initialize the particles' parameter portion...
P <- try(
- particles(object,Np=Np,center=theta,sd=sigma.n*var.factor),
+ particles(obj,Np=Np,center=theta,sd=sigma.n*var.factor),
silent=FALSE
)
if (inherits(P,'try-error'))
@@ -182,7 +231,7 @@
## run the particle filter
x <- try(
pfilter.internal(
- object=object,
+ object=obj,
params=P,
tol=tol,
warn=warn,
@@ -219,48 +268,20 @@
}
- coef(object) <- theta
+ coef(obj) <- theta
- if (Nmif>0)
- new(
- "mif",
- as(object,"pomp"),
- ivps=ivps,
- pars=pars,
- Nmif=as.integer(Nmif),
- particles=particles,
- alg.pars=list(
- Np=Np,
- cooling.factor=cooling.factor,
- var.factor=var.factor,
- ic.lag=ic.lag
- ),
- random.walk.sd=sigma[rw.names],
- pred.mean=x$pred.mean,
- pred.var=x$pred.var,
- filter.mean=x$filter.mean,
- conv.rec=conv.rec,
- eff.sample.size=x$eff.sample.size,
- cond.loglik=x$cond.loglik,
- loglik=x$loglik
- )
- else
- new(
- "mif",
- as(object,"pomp"),
- ivps=ivps,
- pars=pars,
- Nmif=0L,
- particles=particles,
- alg.pars=list(
- Np=Np,
- cooling.factor=cooling.factor,
- var.factor=var.factor,
- ic.lag=ic.lag
- ),
- random.walk.sd=sigma[rw.names],
- conv.rec=conv.rec
- )
+ if (Nmif>0) {
+ obj at Nmif <- as.integer(Nmif)
+ obj at conv.rec <- conv.rec
+ obj at pred.mean <- x$pred.mean
+ obj at pred.var <- x$pred.var
+ obj at filter.mean <- x$filter.mean
+ obj at eff.sample.size <- x$eff.sample.size
+ obj at cond.loglik <- x$cond.loglik
+ obj at loglik <- x$loglik
+ }
+
+ obj
}
setMethod(
Modified: pkg/tests/ou2-mif.R
===================================================================
--- pkg/tests/ou2-mif.R 2009-11-16 11:18:55 UTC (rev 178)
+++ pkg/tests/ou2-mif.R 2009-12-09 12:16:05 UTC (rev 179)
@@ -160,3 +160,23 @@
fit <- mif(ou2,Nmif=3,rw.sd=c(alpha.1=0.1,alpha.4=0.1),Np=1000,cooling.factor=0.98,var.factor=1,ic.lag=2)
fit <- continue(fit,Nmif=2,Np=2000)
fit <- continue(fit,ivps=c("x1.0"),rw.sd=c(alpha.1=0.1,alpha.4=0.1,x1.0=5,x2.0=5),Nmif=2)
+
+pp <- particles(fit,Np=10,center=coef(fit),sd=abs(0.1*coef(fit)))
+fit <- mif(
+ fit,
+ Nmif=10,
+ Np=1000,
+ particles=function(Np,center,sd,...){
+ matrix(
+ data=runif(
+ n=Np*length(center),
+ min=center-sd,
+ max=center+sd
+ ),
+ nrow=length(center),
+ ncol=Np,
+ dimnames=list(names(center),NULL)
+ )
+ }
+ )
+pp <- particles(fit,Np=10,center=coef(fit),sd=abs(0.1*coef(fit)))
More information about the pomp-commits
mailing list