[Pomp-commits] r822 - in pkg/pomp: . R inst man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Feb 4 13:43:18 CET 2013
Author: kingaa
Date: 2013-02-04 13:43:18 +0100 (Mon, 04 Feb 2013)
New Revision: 822
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/R/mif-class.R
pkg/pomp/R/mif.R
pkg/pomp/R/pfilter.R
pkg/pomp/inst/NEWS
pkg/pomp/man/mif-class.Rd
pkg/pomp/man/mif.Rd
pkg/pomp/tests/bbs-trajmatch.Rout.save
pkg/pomp/tests/bbs.Rout.save
pkg/pomp/tests/blowflies.Rout.save
pkg/pomp/tests/dacca.Rout.save
pkg/pomp/tests/dimchecks.Rout.save
pkg/pomp/tests/fhn.Rout.save
pkg/pomp/tests/filtfail.Rout.save
pkg/pomp/tests/gillespie.Rout.save
pkg/pomp/tests/gompertz.R
pkg/pomp/tests/gompertz.Rout.save
pkg/pomp/tests/logistic.Rout.save
pkg/pomp/tests/ou2-bsmc.Rout.save
pkg/pomp/tests/ou2-forecast.R
pkg/pomp/tests/ou2-forecast.Rout.save
pkg/pomp/tests/ou2-icfit.R
pkg/pomp/tests/ou2-icfit.Rout.save
pkg/pomp/tests/ou2-kalman.Rout.save
pkg/pomp/tests/ou2-mif-fp.R
pkg/pomp/tests/ou2-mif-fp.Rout.save
pkg/pomp/tests/ou2-mif.R
pkg/pomp/tests/ou2-mif.Rout.save
pkg/pomp/tests/ou2-mif2.R
pkg/pomp/tests/ou2-mif2.Rout.save
pkg/pomp/tests/ou2-nlf.Rout.save
pkg/pomp/tests/ou2-pmcmc.Rout.save
pkg/pomp/tests/ou2-probe.Rout.save
pkg/pomp/tests/ou2-procmeas.Rout.save
pkg/pomp/tests/ou2-simulate.Rout.save
pkg/pomp/tests/ou2-trajmatch.Rout.save
pkg/pomp/tests/pfilter.Rout.save
pkg/pomp/tests/pomppomp.Rout.save
pkg/pomp/tests/ricker-bsmc.Rout.save
pkg/pomp/tests/ricker-probe.Rout.save
pkg/pomp/tests/ricker-spect.Rout.save
pkg/pomp/tests/ricker.Rout.save
pkg/pomp/tests/rw2.Rout.save
pkg/pomp/tests/sir.Rout.save
pkg/pomp/tests/skeleton.Rout.save
pkg/pomp/tests/steps.Rout.save
pkg/pomp/tests/synlik.Rout.save
pkg/pomp/tests/verhulst.Rout.save
Log:
- change handling of MIF cooling
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/DESCRIPTION 2013-02-04 12:43:18 UTC (rev 822)
@@ -2,7 +2,7 @@
Type: Package
Title: Statistical inference for partially observed Markov processes
Version: 0.44-1
-Date: 2013-01-15
+Date: 2013-02-04
Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall, Helen Wearing, Matthew J. Ferrari, Michael Lavine, Daniel C. Reuman
Maintainer: Aaron A. King <kingaa at umich.edu>
URL: http://pomp.r-forge.r-project.org
Modified: pkg/pomp/R/mif-class.R
===================================================================
--- pkg/pomp/R/mif-class.R 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/R/mif-class.R 2013-02-04 12:43:18 UTC (rev 822)
@@ -10,7 +10,7 @@
particles = 'function',
var.factor='numeric',
ic.lag='integer',
- cooling.factor='numeric',
+ cooling.type='character',
cooling.fraction='numeric',
method='character',
random.walk.sd = 'numeric',
Modified: pkg/pomp/R/mif.R
===================================================================
--- pkg/pomp/R/mif.R 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/R/mif.R 2013-02-04 12:43:18 UTC (rev 822)
@@ -16,6 +16,43 @@
)
}
+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))
+ )
+}
+
mif.cooling <- function (factor, n) { # default geometric cooling schedule
alpha <- factor^(n-1)
list(alpha=alpha,gamma=alpha^2)
@@ -23,8 +60,8 @@
mif2.cooling <- function (frac, nt, m, n) { # cooling schedule for mif2
## 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))
+ scal <- (50*n*frac-1)/(1-frac)
+ alpha <- (1+scal)/(scal+nt+n*(m-1))
list(alpha=alpha)
}
@@ -44,8 +81,8 @@
start, pars, ivps,
particles,
rw.sd,
- Np, cooling.factor, var.factor, ic.lag,
- cooling.fraction,
+ Np, var.factor, ic.lag,
+ cooling.type, cooling.fraction, cooling.factor,
method,
tol, max.fail,
verbose, transform, .ndone = 0,
@@ -164,28 +201,39 @@
)
}
- if (method=="mif2") {
- if (missing(cooling.fraction) || is.na(cooling.fraction))
- stop("mif error: ",sQuote("cooling.fraction")," must be specified for method = ",sQuote("mif2"),call.=FALSE)
- cooling.fraction <- as.numeric(cooling.fraction)
- if ((length(cooling.fraction)!=1)||(cooling.fraction<0)||(cooling.fraction>1))
- 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)
- if (Np[1]!=Np[ntimes+1])
- stop("the first and last values of ",sQuote("Np")," must agree when method = ",sQuote("mif2"))
- } else {
- if (missing(cooling.factor) || is.na(cooling.factor))
- stop("mif error: ",sQuote("cooling.factor")," must be specified",call.=FALSE)
+ ## the following deals with the deprecated option 'cooling.factor'
+ if (!missing(cooling.factor)) {
+ warning(sQuote("cooling.factor")," is deprecated.\n",
+ "See ?mif for instructions on specifying the cooling schedule.",
+ call.=FALSE)
cooling.factor <- as.numeric(cooling.factor)
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 (!missing(cooling.fraction) && !(is.na(cooling.fraction)))
- warning(sQuote("cooling.fraction")," ignored for method != ",sQuote("mif2"),call.=FALSE)
- cooling.fraction <- as.numeric(NA)
+ if (missing(cooling.fraction)) {
+ cooling.fraction <- cooling.factor^50
+ } else {
+ warning("specification of ",sQuote("cooling.factor"),
+ " is overridden by that of ",sQuote("cooling.fraction"),
+ call.=FALSE)
+ }
}
+
+ if (missing(cooling.fraction))
+ stop("mif error: ",sQuote("cooling.fraction")," must be specified",call.=FALSE)
+ cooling.fraction <- as.numeric(cooling.fraction)
+ if ((length(cooling.fraction)!=1)||(cooling.fraction<0)||(cooling.fraction>1))
+ 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
+ )
+
+ if ((method=="mif2")&&(Np[1]!=Np[ntimes+1]))
+ stop("the first and last values of ",sQuote("Np")," must agree when method = ",sQuote("mif2"))
+
if (missing(var.factor))
stop("mif error: ",sQuote("var.factor")," must be specified",call.=FALSE)
if ((length(var.factor)!=1)||(var.factor < 0))
@@ -243,16 +291,7 @@
for (n in seq_len(Nmif)) { ## iterate the filtering
## get the intensity of artificial noise from the cooling schedule
- cool.sched <- try(
- switch(
- method,
- mif2=mif2.cooling(frac=cooling.fraction,nt=1,m=.ndone+n,n=ntimes),
- mif.cooling(factor=cooling.factor,n=.ndone+n)
- ),
- silent=FALSE
- )
- if (inherits(cool.sched,"try-error"))
- stop("mif error: cooling schedule error",call.=FALSE)
+ cool.sched <- cooling(nt=1,m=.ndone+n)
sigma.n <- sigma*cool.sched$alpha
## initialize the parameter portions of the particles
@@ -283,7 +322,7 @@
pred.mean=(n==Nmif),
pred.var=((method=="mif")||(n==Nmif)),
filter.mean=TRUE,
- cooling.fraction=cooling.fraction,
+ cooling=cooling,
cooling.m=.ndone+n,
.mif2=(method=="mif2"),
.rw.sd=sigma.n[pars],
@@ -341,7 +380,7 @@
tol=tol,
conv.rec=conv.rec,
method=method,
- cooling.factor=cooling.factor,
+ cooling.type=cooling.type,
cooling.fraction=cooling.fraction,
paramMatrix=if (method=="mif2") paramMatrix else array(data=numeric(0),dim=c(0,0))
)
@@ -356,8 +395,9 @@
start,
pars, ivps = character(0),
particles, rw.sd,
- Np, ic.lag, var.factor, cooling.factor,
- cooling.fraction,
+ 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"),
@@ -380,6 +420,8 @@
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
@@ -404,6 +446,7 @@
particles=particles,
rw.sd=rw.sd,
Np=Np,
+ cooling.type=cooling.type,
cooling.factor=cooling.factor,
cooling.fraction=cooling.fraction,
var.factor=var.factor,
@@ -445,8 +488,8 @@
start,
pars, ivps,
particles, rw.sd,
- Np, ic.lag, var.factor, cooling.factor,
- cooling.fraction,
+ Np, ic.lag, var.factor,
+ cooling.type, cooling.fraction,
method,
tol,
transform,
@@ -460,7 +503,7 @@
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.factor)) cooling.factor <- object at cooling.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
@@ -478,7 +521,7 @@
particles=particles,
rw.sd=rw.sd,
Np=Np,
- cooling.factor=cooling.factor,
+ cooling.type=cooling.type,
cooling.fraction=cooling.fraction,
var.factor=var.factor,
ic.lag=ic.lag,
Modified: pkg/pomp/R/pfilter.R
===================================================================
--- pkg/pomp/R/pfilter.R 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/R/pfilter.R 2013-02-04 12:43:18 UTC (rev 822)
@@ -38,7 +38,7 @@
pfilter.internal <- function (object, params, Np,
tol, max.fail,
pred.mean, pred.var, filter.mean,
- cooling.fraction, cooling.m, .mif2 = FALSE,
+ cooling, cooling.m, .mif2 = FALSE,
.rw.sd, seed, verbose,
save.states, save.params,
.transform) {
@@ -184,21 +184,10 @@
else
filt.m <- array(data=numeric(0),dim=c(0,0))
- if (mif2) {
- if (missing(cooling.fraction))
- stop("pfilter error: ",sQuote("cooling.fraction")," must be specified for method mif2",call.=FALSE)
- cooling.fraction <- as.numeric(cooling.fraction)
- }
-
for (nt in seq_len(ntimes)) {
if (mif2) {
- cool.sched <- try(
- mif2.cooling(frac=cooling.fraction,nt=nt,m=cooling.m,n=ntimes),
- silent=FALSE
- )
- if (inherits(cool.sched,"try-error"))
- stop("pfilter error: cooling schedule error",call.=FALSE)
+ cool.sched <- cooling(nt=nt,m=cooling.m)
sigma1 <- sigma*cool.sched$alpha
} else {
sigma1 <- sigma
Modified: pkg/pomp/inst/NEWS
===================================================================
--- pkg/pomp/inst/NEWS 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/inst/NEWS 2013-02-04 12:43:18 UTC (rev 822)
@@ -6,6 +6,12 @@
Before, the default behavior has been to stop with an error on the first filtering failure ('max.fail=0').
Now, the default is 'max.fail=Inf', i.e., an error is never triggered.
+ o The implementation of MIF cooling schedules has been changed to make it more general.
+ The cooling schedule is now specified by a 'type' and a 'fraction'.
+ Currently, supported 'cooling.type's include 'geometric' (the old behavior) and 'hyperbolic', i.e., a 1/(1+n) schedule.
+ The 'cooling.fraction' argument specifies the cooling at 50 iterations.
+ That is, if s is the intensity of the random-walk perturbation to parameters at the first iteration ('rw.sd'), then the intensity at iteration 50 is s*cooling.fraction.
+
o Remove all data()-loadable pomp objects.
To load the prebuilt example pomp objects from previous versions, use the new 'pompExample' function.
E.g., instead of 'data(euler.sir)', do 'pompExample("euler.sir")'.
Modified: pkg/pomp/man/mif-class.Rd
===================================================================
--- pkg/pomp/man/mif-class.Rd 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/man/mif-class.Rd 2013-02-04 12:43:18 UTC (rev 822)
@@ -39,13 +39,9 @@
\item{ic.lag}{
the fixed lag used in the estimation of initial-value parameters (IVPs)
}
- \item{cooling.factor}{
- the exponential cooling factor, where \code{0<cooling.factor<1}.
+ \item{cooling.type, cooling.fraction}{
+ the cooling schedule specifications
}
- \item{cooling.fraction}{
- the fraction that determines the cooling schedule when \code{method='mif2'}.
- \code{0<cooling.fraction<1}.
- }
\item{method}{
character; the mif update method.
}
Modified: pkg/pomp/man/mif.Rd
===================================================================
--- pkg/pomp/man/mif.Rd 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/man/mif.Rd 2013-02-04 12:43:18 UTC (rev 822)
@@ -17,13 +17,15 @@
\usage{
mif(object, \dots)
\S4method{mif}{pomp}(object, Nmif = 1, start, pars, ivps = character(0),
- particles, rw.sd, Np, ic.lag, var.factor, cooling.factor,
- cooling.fraction, method = c("mif","unweighted","fp","mif2"),
+ particles, rw.sd, Np, ic.lag, var.factor,
+ cooling.type, cooling.fraction, cooling.factor,
+ method = c("mif","unweighted","fp","mif2"),
tol = 1e-17, max.fail = Inf,
verbose = getOption("verbose"), transform = FALSE, \dots)
\S4method{mif}{pfilterd.pomp}(object, Nmif = 1, Np, tol, \dots)
\S4method{mif}{mif}(object, Nmif, start, pars, ivps,
- particles, rw.sd, Np, ic.lag, var.factor, cooling.factor, cooling.fraction,
+ particles, rw.sd, Np, ic.lag, var.factor,
+ cooling.type, cooling.fraction,
method, tol, transform, \dots)
\S4method{continue}{mif}(object, Nmif = 1, \dots)
}
@@ -85,13 +87,15 @@
the scaling coefficient relating the width of the starting particle distribution to \code{rw.sd}.
In particular, the width of the distribution of particles at the start of the first MIF iteration will be \code{random.walk.sd*var.factor}.
}
- \item{cooling.factor}{
- a positive number not greater than 1;
- the exponential cooling factor, \code{alpha}.
+ \item{cooling.type, cooling.fraction, cooling.factor}{
+ specifications for the cooling schedule, i.e., the manner in which the intensity of the parameter perturbations is reduced with successive filtering iterations.
+ \code{cooling.type} specifies the nature of the cooling schedule.
+ When \code{cooling.type="geometric"}, on the n-th MIF iteration, the relative perturbation intensity is \code{cooling.fraction^(n/50)}.
+ When \code{cooling.type="hyperbolic"}, on the n-th MIF iteration, the relative perturbation intensity is \code{(s+1)(s+n)}, where \code{(s+1)/(s+50)=cooling.fraction}.
+ \code{cooling.fraction} is the relative magnitude of the parameter perturbations after 50 MIF iterations.
+ \code{cooling.factor} is now deprecated:
+ to achieve the old behavior, use \code{cooling.type="geometric"} and \code{cooling.fraction=(cooling.factor)^50}.
}
- \item{cooling.fraction}{
- [EXPLAIN]
- }
\item{method}{
\code{method} sets the update rule used in the algorithm.
\code{method="mif"} uses the iterated filtering update rule (Ionides 2006, 2011);
Modified: pkg/pomp/tests/bbs-trajmatch.Rout.save
===================================================================
--- pkg/pomp/tests/bbs-trajmatch.Rout.save 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/tests/bbs-trajmatch.Rout.save 2013-02-04 12:43:18 UTC (rev 822)
@@ -70,4 +70,4 @@
>
> proc.time()
user system elapsed
- 2.424 0.060 2.516
+ 2.176 0.032 2.241
Modified: pkg/pomp/tests/bbs.Rout.save
===================================================================
--- pkg/pomp/tests/bbs.Rout.save 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/tests/bbs.Rout.save 2013-02-04 12:43:18 UTC (rev 822)
@@ -60,4 +60,4 @@
>
> proc.time()
user system elapsed
- 2.704 0.100 2.827
+ 2.648 0.048 2.732
Modified: pkg/pomp/tests/blowflies.Rout.save
===================================================================
--- pkg/pomp/tests/blowflies.Rout.save 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/tests/blowflies.Rout.save 2013-02-04 12:43:18 UTC (rev 822)
@@ -72,4 +72,4 @@
>
> proc.time()
user system elapsed
- 1.316 0.064 1.397
+ 1.124 0.056 1.204
Modified: pkg/pomp/tests/dacca.Rout.save
===================================================================
--- pkg/pomp/tests/dacca.Rout.save 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/tests/dacca.Rout.save 2013-02-04 12:43:18 UTC (rev 822)
@@ -52,4 +52,4 @@
>
> proc.time()
user system elapsed
- 3.544 0.096 3.661
+ 3.304 0.044 3.379
Modified: pkg/pomp/tests/dimchecks.Rout.save
===================================================================
--- pkg/pomp/tests/dimchecks.Rout.save 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/tests/dimchecks.Rout.save 2013-02-04 12:43:18 UTC (rev 822)
@@ -164,4 +164,4 @@
>
> proc.time()
user system elapsed
- 0.576 0.032 0.619
+ 0.452 0.040 0.517
Modified: pkg/pomp/tests/fhn.Rout.save
===================================================================
--- pkg/pomp/tests/fhn.Rout.save 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/tests/fhn.Rout.save 2013-02-04 12:43:18 UTC (rev 822)
@@ -93,4 +93,4 @@
>
> proc.time()
user system elapsed
- 1.172 0.064 1.359
+ 0.908 0.060 1.204
Modified: pkg/pomp/tests/filtfail.Rout.save
===================================================================
--- pkg/pomp/tests/filtfail.Rout.save 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/tests/filtfail.Rout.save 2013-02-04 12:43:18 UTC (rev 822)
@@ -121,4 +121,4 @@
>
> proc.time()
user system elapsed
- 0.460 0.056 0.534
+ 0.456 0.052 0.526
Modified: pkg/pomp/tests/gillespie.Rout.save
===================================================================
--- pkg/pomp/tests/gillespie.Rout.save 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/tests/gillespie.Rout.save 2013-02-04 12:43:18 UTC (rev 822)
@@ -133,4 +133,4 @@
>
> proc.time()
user system elapsed
- 4.400 0.056 4.480
+ 2.392 0.040 2.459
Modified: pkg/pomp/tests/gompertz.R
===================================================================
--- pkg/pomp/tests/gompertz.R 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/tests/gompertz.R 2013-02-04 12:43:18 UTC (rev 822)
@@ -20,7 +20,8 @@
start=guess,
Nmif=5,Np=1000,
transform=TRUE,
- ic.lag=1,var.factor=1,cooling.factor=0.99,
+ ic.lag=1,var.factor=1,
+ cooling.fraction=0.99^50,
rw.sd=c(r=0.02,K=0.02)
)
)
@@ -30,7 +31,8 @@
po,
Nmif=5,Np=1000,
transform=TRUE,
- ic.lag=1,var.factor=1,cooling.factor=0.99,
+ ic.lag=1,var.factor=1,
+ cooling.fraction=0.99^50,
rw.sd=c(r=0.02,K=0.02)
)
coef(mf,transform=TRUE)
Modified: pkg/pomp/tests/gompertz.Rout.save
===================================================================
--- pkg/pomp/tests/gompertz.Rout.save 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/tests/gompertz.Rout.save 2013-02-04 12:43:18 UTC (rev 822)
@@ -55,7 +55,8 @@
+ start=guess,
+ Nmif=5,Np=1000,
+ transform=TRUE,
-+ ic.lag=1,var.factor=1,cooling.factor=0.99,
++ ic.lag=1,var.factor=1,
++ cooling.fraction=0.99^50,
+ rw.sd=c(r=0.02,K=0.02)
+ )
+ )
@@ -67,7 +68,8 @@
+ po,
+ Nmif=5,Np=1000,
+ transform=TRUE,
-+ ic.lag=1,var.factor=1,cooling.factor=0.99,
++ ic.lag=1,var.factor=1,
++ cooling.fraction=0.99^50,
+ rw.sd=c(r=0.02,K=0.02)
+ )
> coef(mf,transform=TRUE)
@@ -133,4 +135,4 @@
>
> proc.time()
user system elapsed
- 1.512 0.036 1.573
+ 1.488 0.048 1.562
Modified: pkg/pomp/tests/logistic.Rout.save
===================================================================
--- pkg/pomp/tests/logistic.Rout.save 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/tests/logistic.Rout.save 2013-02-04 12:43:18 UTC (rev 822)
@@ -124,4 +124,4 @@
>
> proc.time()
user system elapsed
- 1.052 0.048 1.219
+ 0.788 0.056 0.968
Modified: pkg/pomp/tests/ou2-bsmc.Rout.save
===================================================================
--- pkg/pomp/tests/ou2-bsmc.Rout.save 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/tests/ou2-bsmc.Rout.save 2013-02-04 12:43:18 UTC (rev 822)
@@ -64,7 +64,7 @@
> post <- smc$post
>
> print(etime <- toc-tic)
-Time difference of 3.497756 secs
+Time difference of 2.994066 secs
>
> print(
+ cbind(
@@ -100,4 +100,4 @@
>
> proc.time()
user system elapsed
- 5.452 0.080 5.564
+ 4.876 0.048 4.960
Modified: pkg/pomp/tests/ou2-forecast.R
===================================================================
--- pkg/pomp/tests/ou2-forecast.R 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/tests/ou2-forecast.R 2013-02-04 12:43:18 UTC (rev 822)
@@ -28,7 +28,7 @@
mse[,k,] <- bias^2+sd^2 ## mean squared error
}
-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 <- mif(ou2,Nmif=3,rw.sd=c(alpha.1=0.1,alpha.4=0.1),Np=1000,cooling.fraction=0.98^50,var.factor=1,ic.lag=2)
pf <- pfilter(fit,save.states=TRUE,save.params=TRUE)
pdf(file="ou2-forecast.pdf")
Modified: pkg/pomp/tests/ou2-forecast.Rout.save
===================================================================
--- pkg/pomp/tests/ou2-forecast.Rout.save 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/tests/ou2-forecast.Rout.save 2013-02-04 12:43:18 UTC (rev 822)
@@ -51,7 +51,7 @@
+ mse[,k,] <- bias^2+sd^2 ## mean squared error
+ }
>
-> 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 <- mif(ou2,Nmif=3,rw.sd=c(alpha.1=0.1,alpha.4=0.1),Np=1000,cooling.fraction=0.98^50,var.factor=1,ic.lag=2)
> pf <- pfilter(fit,save.states=TRUE,save.params=TRUE)
>
> pdf(file="ou2-forecast.pdf")
@@ -63,4 +63,4 @@
>
> proc.time()
user system elapsed
- 1.928 0.060 2.108
+ 1.528 0.044 1.689
Modified: pkg/pomp/tests/ou2-icfit.R
===================================================================
--- pkg/pomp/tests/ou2-icfit.R 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/tests/ou2-icfit.R 2013-02-04 12:43:18 UTC (rev 822)
@@ -20,7 +20,7 @@
x1.0=1,x2.0=1
),
Np=1000,
- cooling.factor=1,
+ cooling.fraction=1,
max.fail=10
)
@@ -35,7 +35,7 @@
x1.0=1,x2.0=1
),
Np=1000,
- cooling.factor=1,
+ cooling.fraction=1,
max.fail=10
)
@@ -50,7 +50,7 @@
x1.0=1,x2.0=1
),
Np=1000,
- cooling.factor=1,
+ cooling.fraction=1,
max.fail=10
)
Modified: pkg/pomp/tests/ou2-icfit.Rout.save
===================================================================
--- pkg/pomp/tests/ou2-icfit.Rout.save 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/tests/ou2-icfit.Rout.save 2013-02-04 12:43:18 UTC (rev 822)
@@ -43,7 +43,7 @@
+ x1.0=1,x2.0=1
+ ),
+ Np=1000,
-+ cooling.factor=1,
++ cooling.fraction=1,
+ max.fail=10
+ )
Warning message:
@@ -60,7 +60,7 @@
+ x1.0=1,x2.0=1
+ ),
+ Np=1000,
-+ cooling.factor=1,
++ cooling.fraction=1,
+ max.fail=10
+ )
Warning message:
@@ -77,7 +77,7 @@
+ x1.0=1,x2.0=1
+ ),
+ Np=1000,
-+ cooling.factor=1,
++ cooling.fraction=1,
+ max.fail=10
+ )
>
@@ -109,4 +109,4 @@
>
> proc.time()
user system elapsed
- 17.629 0.060 17.749
+ 16.897 0.048 17.022
Modified: pkg/pomp/tests/ou2-kalman.Rout.save
===================================================================
--- pkg/pomp/tests/ou2-kalman.Rout.save 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/tests/ou2-kalman.Rout.save 2013-02-04 12:43:18 UTC (rev 822)
@@ -166,7 +166,7 @@
117 function evaluations used
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 5.160968 secs
+Time difference of 5.084913 secs
> tic <- Sys.time()
> print(loglik.mle <- -kalm.fit1$value,digits=4)
[1] -477.2
@@ -190,4 +190,4 @@
>
> proc.time()
user system elapsed
- 5.680 0.044 5.752
+ 5.596 0.044 5.680
Modified: pkg/pomp/tests/ou2-mif-fp.R
===================================================================
--- pkg/pomp/tests/ou2-mif-fp.R 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/tests/ou2-mif-fp.R 2013-02-04 12:43:18 UTC (rev 822)
@@ -14,30 +14,32 @@
mif1 <- mif(ou2,Nmif=100,start=guess1,
pars=c('alpha.2','alpha.3'),ivps=c('x1.0','x2.0'),
rw.sd=c(
- x1.0=5,x2.0=5,
- alpha.2=0.1,alpha.3=0.1
- ),
- Np=1000,
- var.factor=1,
- ic.lag=10,
- cooling.factor=0.95,
- max.fail=100,
- method="fp"
+ x1.0=5,x2.0=5,
+ alpha.2=0.1,alpha.3=0.1
+ ),
+ Np=1000,
+ var.factor=1,
+ ic.lag=10,
+ cooling.type="geometric",
+ cooling.fraction=0.95^50,
+ max.fail=100,
+ method="fp"
)
mif2 <- mif(ou2,Nmif=100,start=guess2,
pars=c('alpha.2','alpha.3'),ivps=c('x1.0','x2.0'),
rw.sd=c(
- x1.0=5,x2.0=5,
- alpha.2=0.1,alpha.3=0.1
- ),
- Np=1000,
- var.factor=1,
- ic.lag=10,
- cooling.factor=0.95,
- max.fail=100,
- method="fp"
- )
+ x1.0=5,x2.0=5,
+ alpha.2=0.1,alpha.3=0.1
+ ),
+ Np=1000,
+ var.factor=1,
+ ic.lag=10,
+ cooling.type="geometric",
+ cooling.fraction=0.95^50,
+ max.fail=100,
+ method="fp"
+ )
compare.mif(list(mif1,mif2))
Modified: pkg/pomp/tests/ou2-mif-fp.Rout.save
===================================================================
--- pkg/pomp/tests/ou2-mif-fp.Rout.save 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/tests/ou2-mif-fp.Rout.save 2013-02-04 12:43:18 UTC (rev 822)
@@ -37,30 +37,32 @@
> mif1 <- mif(ou2,Nmif=100,start=guess1,
+ pars=c('alpha.2','alpha.3'),ivps=c('x1.0','x2.0'),
+ rw.sd=c(
-+ x1.0=5,x2.0=5,
-+ alpha.2=0.1,alpha.3=0.1
-+ ),
-+ Np=1000,
-+ var.factor=1,
-+ ic.lag=10,
-+ cooling.factor=0.95,
-+ max.fail=100,
-+ method="fp"
++ x1.0=5,x2.0=5,
++ alpha.2=0.1,alpha.3=0.1
++ ),
++ Np=1000,
++ var.factor=1,
++ ic.lag=10,
++ cooling.type="geometric",
++ cooling.fraction=0.95^50,
++ max.fail=100,
++ method="fp"
+ )
>
> mif2 <- mif(ou2,Nmif=100,start=guess2,
+ pars=c('alpha.2','alpha.3'),ivps=c('x1.0','x2.0'),
+ rw.sd=c(
-+ x1.0=5,x2.0=5,
-+ alpha.2=0.1,alpha.3=0.1
-+ ),
-+ Np=1000,
-+ var.factor=1,
-+ ic.lag=10,
-+ cooling.factor=0.95,
-+ max.fail=100,
-+ method="fp"
-+ )
++ x1.0=5,x2.0=5,
++ alpha.2=0.1,alpha.3=0.1
++ ),
++ Np=1000,
++ var.factor=1,
++ ic.lag=10,
++ cooling.type="geometric",
++ cooling.fraction=0.95^50,
++ max.fail=100,
++ method="fp"
++ )
>
> compare.mif(list(mif1,mif2))
>
@@ -70,4 +72,4 @@
>
> proc.time()
user system elapsed
- 20.065 0.076 20.346
+ 19.653 0.024 19.956
Modified: pkg/pomp/tests/ou2-mif.R
===================================================================
--- pkg/pomp/tests/ou2-mif.R 2013-01-21 23:14:11 UTC (rev 821)
+++ pkg/pomp/tests/ou2-mif.R 2013-02-04 12:43:18 UTC (rev 822)
@@ -27,6 +27,7 @@
Np=1000,
var.factor=1,
ic.lag=10,
+ cooling.type="geometric",
cooling.factor=0.95,
max.fail=100
)
@@ -41,7 +42,8 @@
Np=1000,
var.factor=1,
ic.lag=10,
- cooling.factor=0.95,
+ cooling.type="geometric",
+ cooling.fraction=0.95^50,
max.fail=100
)
@@ -61,7 +63,8 @@
pars=c("alpha.1","alpha.4","x1.0"),
ivps=c("x1.0","x2.0"),
rw.sd=c(x1.0=5,x2.0=5,alpha.1=0.1,alpha.4=0.2,alpha.3=0),
- Np=100,cooling.factor=0.95,ic.lag=10,var.factor=1
+ Np=100,cooling.type="geometric",cooling.fraction=0.95^50,
+ ic.lag=10,var.factor=1
)
)
@@ -72,7 +75,9 @@
pars=c("alpha.1","alpha.4"),
ivps=c("x1.0","x2.0"),
rw.sd=c(x1.0=5,x2.0=5,alpha.1=0,alpha.4=0.2,alpha.3=0),
- Np=100,cooling.factor=0.95,ic.lag=10,var.factor=1
+ Np=100,
+ cooling.type="geometric",cooling.fraction=0.95^50,
+ ic.lag=10,var.factor=1
)
)
@@ -82,7 +87,8 @@
Nmif=1,
ivps=c("x1.0","x2.0"),
rw.sd=c(alpha.1=0.1,alpha.4=0.2,alpha.3=0),
- Np=100,cooling.factor=0.95,ic.lag=10,var.factor=1
+ Np=100,cooling.type="geometric",cooling.fraction=0.95^50,
+ cooling.factor=0.95,ic.lag=10,var.factor=1
)
)
@@ -102,7 +108,8 @@
Nmif=1,
ivps=c("x1.0","x2.0"),
rw.sd=c(x1.0=5,x2.0=5,alpha.1=0,alpha.4=0.2,alpha.3=0),
- Np=-10,cooling.factor=0.95,ic.lag=10,var.factor=1
+ Np=-10,cooling.type="geometric",cooling.fraction=0.95^50,
+ ic.lag=10,var.factor=1
)
)
@@ -112,7 +119,8 @@
Nmif=-3,
ivps=c("x1.0","x2.0"),
rw.sd=c(x1.0=5,x2.0=5,alpha.1=0,alpha.4=0.2,alpha.3=0),
- Np=11.6,cooling.factor=0.95,ic.lag=10,var.factor=1
+ Np=11.6,cooling.type="geometric",cooling.fraction=0.95^50,
+ ic.lag=10,var.factor=1
)
)
@@ -123,7 +131,8 @@
start=c(alpha.1=0.9,alpha.2=0,alpha.3=0,alpha.4=-Inf,sigma.1=1,sigma.2=0,sigma.3=2,tau=1,x1.0=50,x2.0=-50),
ivps=c("x1.0","x2.0"),
rw.sd=c(x1.0=5,x2.0=5,alpha.1=0,alpha.4=0.2,alpha.3=0),
- Np=11,cooling.factor=0.95,ic.lag=10,var.factor=1
+ Np=11,cooling.type="geometric",cooling.fraction=0.95^50,
+ ic.lag=10,var.factor=1
)
)
@@ -134,7 +143,8 @@
start=c(alpha.1=0.9,alpha.2=0,alpha.3=0,alpha.4=0.99,sigma.1=1,sigma.2=0,sigma.3=2,tau=1,x1.0=50,x2.0=NaN),
ivps=c("x1.0","x2.0"),
rw.sd=c(x1.0=5,x2.0=5,alpha.1=0,alpha.4=0.2,alpha.3=0),
- Np=11,cooling.factor=0.95,ic.lag=10,var.factor=1
+ Np=11,cooling.type="geometric",cooling.fraction=0.95^50,
+ ic.lag=10,var.factor=1
)
)
@@ -144,7 +154,8 @@
pars=c("alpha.2","alpha.3"),
ivps=c("x1.0","x2.0"),
rw.sd=c(x1.0=5,x2.0=5,alpha.2=0.1,alpha.3=0.2,alpha.3=0),
- Np=100,cooling.factor=0.95,ic.lag=10,var.factor=1
+ Np=100,cooling.type="geometric",cooling.fraction=0.95^50,
+ ic.lag=10,var.factor=1
)
fit <- mif(
fit,
@@ -152,7 +163,8 @@
ivps=c("x1.0","x2.0"),
rw.sd=c(x1.0=5,x2.0=5,alpha.2=0.1,alpha.3=0.2),
Np=function(k)if(k<10) 2000 else 500,
- cooling.factor=0.95,ic.lag=10,var.factor=1
+ cooling.type="geometric",cooling.fraction=0.95^50,
+ ic.lag=10,var.factor=1
)
fit <- continue(fit)
fit <- continue(fit,Nmif=2)
@@ -163,7 +175,7 @@
s <- coef(fit)
s[2] <- 0.01
fit <- mif(fit,Nmif=3,start=s)
-fit <- mif(ou2,Nmif=3,rw.sd=c(alpha.2=0.1,alpha.3=0.1),Np=1000,cooling.factor=0.98,var.factor=1,ic.lag=2)
+fit <- mif(ou2,Nmif=3,rw.sd=c(alpha.2=0.1,alpha.3=0.1),Np=1000,cooling.type="geometric",cooling.fraction=0.98^50,var.factor=1,ic.lag=2)
fit <- continue(fit,Nmif=2,Np=2000)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 822
More information about the pomp-commits
mailing list