[Pomp-commits] r824 - in branches/mif2: . R tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Feb 5 17:22:39 CET 2013
Author: nxdao2000
Date: 2013-02-05 17:22:39 +0100 (Tue, 05 Feb 2013)
New Revision: 824
Modified:
branches/mif2/
branches/mif2/.Rbuildignore
branches/mif2/R/mif.R
branches/mif2/R/pfilter.R
branches/mif2/tests/ou2-mif2.R
Log:
change for different cooling scheme and update for branch only
Property changes on: branches/mif2
___________________________________________________________________
Added: svn:ignore
+ .Rproj.user
.Rhistory
.RData
Modified: branches/mif2/.Rbuildignore
===================================================================
--- branches/mif2/.Rbuildignore 2013-02-05 15:21:39 UTC (rev 823)
+++ branches/mif2/.Rbuildignore 2013-02-05 16:22:39 UTC (rev 824)
@@ -2,3 +2,5 @@
inst/doc/(.+?)\.bst$
inst/doc/(.+?)\.R$
inst/doc/(.+?)\.png$
+^.*\.Rproj$
+^\.Rproj\.user$
Modified: branches/mif2/R/mif.R
===================================================================
--- branches/mif2/R/mif.R 2013-02-05 15:21:39 UTC (rev 823)
+++ branches/mif2/R/mif.R 2013-02-05 16:22:39 UTC (rev 824)
@@ -25,7 +25,8 @@
## frac is the fraction of cooling after 50 iterations
cooling.scalar <- (50*n*frac-1)/(1-frac)
alpha <- (1+cooling.scalar)/(cooling.scalar+nt+n*(m-1))
- list(alpha=alpha)
+
+ list(alpha=alpha,gamma=alpha^2)
}
powerlaw.cooling <- function (init = 1, delta = 0.1, eps = (1-delta)/2, n) {
@@ -172,7 +173,7 @@
stop("mif error: ",sQuote("cooling.fraction")," must be a number between 0 and 1",call.=FALSE)
if (!missing(cooling.factor) && !(is.na(cooling.factor)))
warning(sQuote("cooling.factor")," ignored for method = ",sQuote("mif2"),call.=FALSE)
- cooling.factor <- as.numeric(NA)
+ cooling.factor <- as.numeric(cooling.factor)
if (Np[1]!=Np[ntimes+1])
stop("the first and last values of ",sQuote("Np")," must agree when method = ",sQuote("mif2"))
} else {
@@ -183,7 +184,7 @@
stop("mif error: ",sQuote("cooling.factor")," must be a number between 0 and 1",call.=FALSE)
if (!missing(cooling.fraction) && !(is.na(cooling.fraction)))
warning(sQuote("cooling.fraction")," ignored for method != ",sQuote("mif2"),call.=FALSE)
- cooling.fraction <- as.numeric(NA)
+ cooling.fraction <- as.numeric(cooling.fraction)
}
if (missing(var.factor))
@@ -247,6 +248,8 @@
switch(
method,
mif2=mif2.cooling(frac=cooling.fraction,nt=1,m=.ndone+n,n=ntimes),
+ mif4=mif2.cooling(frac=cooling.fraction,nt=round((.ndone+n)/2),m=.ndone+n,n=ntimes),
+ mif3=mif.cooling(factor=cooling.factor,n=.ndone+n),
mif.cooling(factor=cooling.factor,n=.ndone+n)
),
silent=FALSE
@@ -268,24 +271,24 @@
if (inherits(P,"try-error"))
stop("mif error: error in ",sQuote("particles"),call.=FALSE)
- if ((method=="mif2") && ((n>1) || have.parmat)) {
+ if (((method=="mif2")||(method=="mif3")) && ((n>1) || have.parmat)) {
## use pre-existing particle matrix
P[pars,] <- paramMatrix[pars,]
}
pfp <- try(
- pfilter.internal(
+ pfilter.internal(
object=obj,
params=P,
Np=Np,
tol=tol,
max.fail=max.fail,
pred.mean=(n==Nmif),
- pred.var=((method=="mif")||(n==Nmif)),
+ pred.var=((method=="mif")||(method=="mif4")||(n==Nmif)),
filter.mean=TRUE,
cooling.fraction=cooling.fraction,
cooling.m=.ndone+n,
- .mif2=(method=="mif2"),
+ .mif2=((method=="mif2")||(method=="mif3")),
.rw.sd=sigma.n[pars],
.transform=transform,
save.states=FALSE,
@@ -303,6 +306,9 @@
mif={ # original Ionides et al. (2006) average
theta <- .Call(mif_update,pfp,theta,cool.sched$gamma,var.factor,sigma,pars)
},
+ mif4={ # 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])
},
@@ -313,6 +319,10 @@
paramMatrix <- pfp at paramMatrix
theta[pars] <- rowMeans(paramMatrix[pars,,drop=FALSE])
},
+ mif3={ # "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]
@@ -343,7 +353,7 @@
method=method,
cooling.factor=cooling.factor,
cooling.fraction=cooling.fraction,
- paramMatrix=if (method=="mif2") paramMatrix else array(data=numeric(0),dim=c(0,0))
+ paramMatrix=if ((method=="mif2")||(method=="mif3")) paramMatrix else array(data=numeric(0),dim=c(0,0))
)
}
@@ -358,7 +368,7 @@
particles, rw.sd,
Np, ic.lag, var.factor, cooling.factor,
cooling.fraction,
- method = c("mif","unweighted","fp","mif2"),
+ method = c("mif","unweighted","fp","mif2","mif3","mif4"),
tol = 1e-17, max.fail = Inf,
verbose = getOption("verbose"),
transform = FALSE,
Modified: branches/mif2/R/pfilter.R
===================================================================
--- branches/mif2/R/pfilter.R 2013-02-05 15:21:39 UTC (rev 823)
+++ branches/mif2/R/pfilter.R 2013-02-05 16:22:39 UTC (rev 824)
@@ -192,7 +192,7 @@
for (nt in seq_len(ntimes)) {
- if (mif2) {
+ if ((mif2==T) && (cooling.fraction>0)) {
cool.sched <- try(
mif2.cooling(frac=cooling.fraction,nt=nt,m=cooling.m,n=ntimes),
silent=FALSE
@@ -200,6 +200,7 @@
if (inherits(cool.sched,"try-error"))
stop("pfilter error: cooling schedule error",call.=FALSE)
sigma1 <- sigma*cool.sched$alpha
+ sigma1 <- sigma
} else {
sigma1 <- sigma
}
Modified: branches/mif2/tests/ou2-mif2.R
===================================================================
--- branches/mif2/tests/ou2-mif2.R 2013-02-05 15:21:39 UTC (rev 823)
+++ branches/mif2/tests/ou2-mif2.R 2013-02-05 16:22:39 UTC (rev 824)
@@ -19,9 +19,9 @@
Np=1000,
var.factor=1,
ic.lag=10,
- cooling.factor=0.95,
+ cooling.factor=0,
cooling.fraction=0.05,
- method="mif2",
+ method="mif4",
tol=1e-8
)
@@ -35,13 +35,13 @@
var.factor=1,
ic.lag=10,
cooling.factor=0.95,
- cooling.fraction=0.5,
+ cooling.fraction=0,
max.fail=100,
- method="mif",
+ method="mif3",
tol=1e-8
)
-compare.mif(list(mif1a,mif2a))
+#compare.mif(list(mif1a,mif2a))
set.seed(64857673L)
mif1b <- mif(ou2,Nmif=50,start=guess1,
@@ -53,8 +53,8 @@
Np=1000,
var.factor=1,
ic.lag=10,
- cooling.factor=0.95,
- cooling.fraction=0.05,
+ cooling.factor=0,
+ cooling.fraction=0.15,
method="mif2"
)
mif1b <- continue(mif1b,Nmif=50)
@@ -69,14 +69,16 @@
var.factor=1,
ic.lag=10,
cooling.whatsit=200,
+ cooling.fraction=0,
cooling.factor=0.95,
max.fail=100,
method="mif"
)
mif2b <- continue(mif2b,Nmif=50)
-compare.mif(list(mif1b,mif2b))
+compare.mif(list(mif1a,mif1b,mif2a,mif2b))
+
compare.mif(list(mif1a,mif1b))
compare.mif(list(mif2a,mif2b))
More information about the pomp-commits
mailing list