[Pomp-commits] r733 - branches/mif2/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jun 18 20:07:51 CEST 2012
Author: nxdao2000
Date: 2012-06-18 20:07:50 +0200 (Mon, 18 Jun 2012)
New Revision: 733
Modified:
branches/mif2/R/mif.R
Log:
an attempt to write mif2 option for mif
Modified: branches/mif2/R/mif.R
===================================================================
--- branches/mif2/R/mif.R 2012-06-18 14:35:01 UTC (rev 732)
+++ branches/mif2/R/mif.R 2012-06-18 18:07:50 UTC (rev 733)
@@ -20,7 +20,10 @@
alpha <- factor^(n-1)
list(alpha=alpha,gamma=alpha^2)
}
-
+mif.cooling2 <- function (factor, n) { # default cooling schedule
+ alpha <- (1+factor)/(n+factor)
+ list(alpha=alpha,gamma=alpha^2)
+}
powerlaw.cooling <- function (init = 1, delta = 0.1, eps = (1-delta)/2, n) {
m <- init
if (n <= m) { # linear cooling regime
@@ -211,76 +214,143 @@
else
pfp <- obj
- for (n in seq_len(Nmif)) { # main loop
+##Branch mif1 and mif2 here
+if (method!="mif2")
+{
+#Nmif = M
+ for (n in seq_len(Nmif)) { # main loop
+
+ ## compute the cooled sigma
+ cool.sched <- try(
+ mif.cooling(cooling.factor,.ndone+n),
+ silent=FALSE
+ )
+ if (inherits(cool.sched,'try-error'))
+ stop("mif error: cooling schedule error",call.=FALSE)
+ sigma.n <- sigma*cool.sched$alpha
+
+ ## initialize the particles' parameter portion...
+ P <- try(
+ particles(tmp.mif,Np=Np[1],center=theta,sd=sigma.n*var.factor),
+ silent=FALSE
+ )
+ if (inherits(P,'try-error'))
+ stop("mif error: error in ",sQuote("particles"),call.=FALSE)
+
+ ## run the particle filter
+ pfp <- try(
+ pfilter.internal(
+ object=obj,
+ params=P,
+ tol=tol,
+ max.fail=max.fail,
+ pred.mean=(n==Nmif),
+ pred.var=((method=="mif")||(n==Nmif)),
+ filter.mean=TRUE,
+ save.states=FALSE,
+ save.params=FALSE,
+ .rw.sd=sigma.n[pars],
+ verbose=verbose,
+ transform=transform
+ ),
+ silent=FALSE
+ )
+ if (inherits(pfp,'try-error'))
+ stop("mif error: error in ",sQuote("pfilter"),call.=FALSE)
+
+ switch(
+ method,
+ mif={ # mif update rule
+ v <- pfp$pred.var[pars,,drop=FALSE] # the prediction variance
+ v1 <- cool.sched$gamma*(1+var.factor^2)*sigma[pars]^2
+ theta.hat <- cbind(theta[pars],pfp$filter.mean[pars,,drop=FALSE])
+ theta[pars] <- theta[pars]+colSums(apply(theta.hat,1,diff)/t(v))*v1
+ },
+ unweighted={ # unweighted (flat) average
+ #theta.hat <- pfp$filter.mean[pars,,drop=FALSE]
+ #theta[pars] <- rowMeans(theta.hat)
+ theta.hat <- pfp$filter.mean[pars,ntimes,drop=FALSE]
+ theta[pars] <- theta.hat
+ },
+ fp={
+ theta.hat <- pfp$filter.mean[pars,ntimes,drop=FALSE]
+ theta[pars] <- theta.hat
+ },
+ #stop("unrecognized method",sQuote(method))
+ )
+
+ ## update the IVPs using fixed-lag smoothing
+ theta[ivps] <- pfp$filter.mean[ivps,ic.lag]
+
+ ## store a record of this iteration
+ conv.rec[n+1,-c(1,2)] <- theta
+ conv.rec[n,c(1,2)] <- c(pfp$loglik,pfp$nfail)
+
+ if (verbose) cat("MIF iteration ",n," of ",Nmif," completed\n")
+
+ }
+}
+else
+{
+ for (n in seq_len(Nmif)) { # main loop
+
+ ## compute the cooled sigma
+ cool.sched <- try(
+ mif.cooling2(cooling.factor,.ndone),#+n,Np=Np[1],n),
+ silent=FALSE
+ )
+ if (inherits(cool.sched,'try-error'))
+ stop("mif error: cooling schedule error",call.=FALSE)
+ sigma.n <- sigma*cool.sched$alpha
+
+ ## initialize the particles' parameter portion...
+ P <- try(
+ particles(tmp.mif,Np=Np[1],center=theta,sd=sigma.n*var.factor),
+ silent=FALSE
+ )
+ if (inherits(P,'try-error'))
+ stop("mif error: error in ",sQuote("particles"),call.=FALSE)
+
+ ## run the particle filter
+ pfp <- try(
+ pfilter.internal(
+ object=obj,
+ params=P,
+ tol=tol,
+ max.fail=max.fail,
+ pred.mean=(n==Nmif),
+ pred.var=((method=="mif")||(n==Nmif)),
+ filter.mean=TRUE,
+ save.states=FALSE,
+ save.params=FALSE,
+ .rw.sd=sigma.n[pars],
+ verbose=verbose,
+ transform=transform
+ ),
+ silent=FALSE
+ )
+ if (inherits(pfp,'try-error'))
+ stop("mif error: error in ",sQuote("pfilter"),call.=FALSE)
+
+ theta.hat <- pfp$filter.mean[pars,ntimes,drop=FALSE]
+ theta[pars] <- theta.hat
+
+
+ ## update the IVPs using fixed-lag smoothing
+ theta[ivps] <- pfp$filter.mean[ivps,ic.lag]
+
+ ## store a record of this iteration
+ conv.rec[n+1,-c(1,2)] <- theta
+ conv.rec[n,c(1,2)] <- c(pfp$loglik,pfp$nfail)
+
+ if (verbose) cat("MIF iteration ",n," of ",Nmif," completed\n")
+
+ }
+
+
+}
- ## compute the cooled sigma
- cool.sched <- try(
- mif.cooling(cooling.factor,.ndone+n),
- silent=FALSE
- )
- if (inherits(cool.sched,'try-error'))
- stop("mif error: cooling schedule error",call.=FALSE)
- sigma.n <- sigma*cool.sched$alpha
- ## initialize the particles' parameter portion...
- P <- try(
- particles(tmp.mif,Np=Np[1],center=theta,sd=sigma.n*var.factor),
- silent=FALSE
- )
- if (inherits(P,'try-error'))
- stop("mif error: error in ",sQuote("particles"),call.=FALSE)
-
- ## run the particle filter
- pfp <- try(
- pfilter.internal(
- object=obj,
- params=P,
- tol=tol,
- max.fail=max.fail,
- pred.mean=(n==Nmif),
- pred.var=((method=="mif")||(n==Nmif)),
- filter.mean=TRUE,
- save.states=FALSE,
- save.params=FALSE,
- .rw.sd=sigma.n[pars],
- verbose=verbose,
- transform=transform
- ),
- silent=FALSE
- )
- if (inherits(pfp,'try-error'))
- stop("mif error: error in ",sQuote("pfilter"),call.=FALSE)
-
- switch(
- method,
- mif={ # mif update rule
- v <- pfp$pred.var[pars,,drop=FALSE] # the prediction variance
- v1 <- cool.sched$gamma*(1+var.factor^2)*sigma[pars]^2
- theta.hat <- cbind(theta[pars],pfp$filter.mean[pars,,drop=FALSE])
- theta[pars] <- theta[pars]+colSums(apply(theta.hat,1,diff)/t(v))*v1
- },
- unweighted={ # unweighted (flat) average
- theta.hat <- pfp$filter.mean[pars,,drop=FALSE]
- theta[pars] <- rowMeans(theta.hat)
- },
- fp={
- theta.hat <- pfp$filter.mean[pars,ntimes,drop=FALSE]
- theta[pars] <- theta.hat
- },
- stop("unrecognized method",sQuote(method))
- )
-
- ## update the IVPs using fixed-lag smoothing
- theta[ivps] <- pfp$filter.mean[ivps,ic.lag]
-
- ## store a record of this iteration
- conv.rec[n+1,-c(1,2)] <- theta
- conv.rec[n,c(1,2)] <- c(pfp$loglik,pfp$nfail)
-
- if (verbose) cat("MIF iteration ",n," of ",Nmif," completed\n")
-
- }
-
## back transform the parameter estimate if necessary
if (transform)
theta <- partrans(pfp,theta,dir="forward")
More information about the pomp-commits
mailing list