[Pomp-commits] r93 - in pkg: . R data inst inst/doc man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Apr 14 16:43:04 CEST 2009
Author: kingaa
Date: 2009-04-14 16:43:04 +0200 (Tue, 14 Apr 2009)
New Revision: 93
Modified:
pkg/DESCRIPTION
pkg/R/mif.R
pkg/data/ou2.rda
pkg/inst/ChangeLog
pkg/inst/doc/compiled_code_in_pomp.pdf
pkg/inst/doc/intro_to_pomp.Rnw
pkg/inst/doc/intro_to_pomp.pdf
pkg/man/mif.Rd
pkg/tests/ou2-mif.R
pkg/tests/ou2-mif.Rout.save
pkg/tests/sir.Rout.save
Log:
The parameter 'alg.pars' is now deprecated. The four algorithm parameters Np, ic.lag, var.factor, cooling.factor are now specified directly as required arguments to pomp. The 'alg.pars' argument will be removed in a future release. Use of this argument now results in a warning message.
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2009-04-13 14:27:26 UTC (rev 92)
+++ pkg/DESCRIPTION 2009-04-14 14:43:04 UTC (rev 93)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.22-8
-Date: 2009-04-12
+Version: 0.23-1
+Date: 2009-04-14
Author: Aaron A. King, Edward L. Ionides, Carles Martinez Breto, Steve Ellner, Bruce Kendall
Maintainer: Aaron A. King <kingaa at umich.edu>
Description: Inference methods for partially-observed Markov processes
Modified: pkg/R/mif.R
===================================================================
--- pkg/R/mif.R 2009-04-13 14:27:26 UTC (rev 92)
+++ pkg/R/mif.R 2009-04-14 14:43:04 UTC (rev 93)
@@ -20,6 +20,7 @@
pars = NULL, ivps = NULL,
particles = NULL,
rw.sd = NULL, alg.pars = NULL,
+ 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) {
if (is.null(start)) {
@@ -89,16 +90,40 @@
if (is.null(particles))
particles <- object at particles
- if (is.null(alg.pars))
- alg.pars <- object at alg.pars
- if (!all(c('Np','cooling.factor','ic.lag','var.factor')%in%names(alg.pars)))
- stop(
- "mif error: ",sQuote("alg.pars"),
- " must be a named list with elements ",
- sQuote("Np"),",",sQuote("cooling.factor"),",",sQuote("ic.lag"),
- ",and ",sQuote("var.factor"),
- call.=FALSE
- )
+ if (is.null(alg.pars)) {
+ if (is.null(Np))
+ Np <- object at alg.pars$Np
+ 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
+ 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 ((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 ((length(var.factor)!=1)||(var.factor < 0))
+ stop("mif error: ",sQuote("var.factor")," must be a positive number",call.=FALSE)
+ } else { # use of 'alg.pars' is now deprecated
+ warning("mif warning: use of ",sQuote("alg.pars")," is deprecated; see the documentation for ",sQuote("mif"),".",call.=FALSE)
+ if (!all(c('Np','cooling.factor','ic.lag','var.factor')%in%names(alg.pars)))
+ stop(
+ "mif error: ",sQuote("alg.pars"),
+ " must be a named list with elements ",
+ sQuote("Np"),",",sQuote("cooling.factor"),",",sQuote("ic.lag"),
+ ",and ",sQuote("var.factor"),
+ call.=FALSE
+ )
+ Np <- alg.pars$Np
+ cooling.factor <- alg.pars$cooling.factor
+ var.factor <- alg.pars$var.factor
+ ic.lag <- alg.pars$ic.lag
+ }
if (verbose) {
cat("performing",Nmif,"MIF iteration(s) to estimate parameter(s)",
@@ -108,9 +133,9 @@
cat(" using random-walk with SD\n")
print(rw.sd)
cat(
- "using",alg.pars$Np,"particles, variance factor",alg.pars$var.factor,
- "\ninitial condition smoothing lag",alg.pars$ic.lag,
- "and cooling factor",alg.pars$cooling.factor,"\n"
+ "using",Np,"particles, variance factor",var.factor,
+ "\ninitial condition smoothing lag",ic.lag,
+ "and cooling factor",cooling.factor,"\n"
)
}
@@ -139,7 +164,7 @@
## compute the cooled sigma
cool.sched <- try(
- mif.cooling(alg.pars$cooling.factor,.ndone+n),
+ mif.cooling(cooling.factor,.ndone+n),
silent=FALSE
)
if (inherits(cool.sched,'try-error'))
@@ -148,7 +173,7 @@
## initialize the particles' parameter portion...
P <- try(
- particles(object,Np=alg.pars$Np,center=theta,sd=sigma.n*alg.pars$var.factor),
+ particles(object,Np=Np,center=theta,sd=sigma.n*var.factor),
silent=FALSE
)
if (inherits(P,'try-error'))
@@ -175,7 +200,7 @@
if (weighted) { # MIF update rule
v <- x$pred.var[pars,,drop=FALSE] # the prediction variance
- v1 <- cool.sched$gamma*(1+alg.pars$var.factor^2)*sigma[pars]^2
+ v1 <- cool.sched$gamma*(1+var.factor^2)*sigma[pars]^2
theta.hat <- cbind(theta[pars],x$filter.mean[pars,,drop=FALSE])
theta[pars] <- theta[pars]+apply(apply(theta.hat,1,diff)/t(v),2,sum)*v1
} else { # unweighted (flat) average
@@ -184,7 +209,7 @@
}
## update the IVPs using fixed-lag smoothing
- theta[ivps] <- x$filter.mean[ivps,alg.pars$ic.lag]
+ theta[ivps] <- x$filter.mean[ivps,ic.lag]
## store a record of this iteration
conv.rec[n+1,-c(1,2)] <- theta
@@ -204,7 +229,12 @@
pars=pars,
Nmif=as.integer(Nmif),
particles=particles,
- alg.pars=alg.pars,
+ 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,
@@ -222,7 +252,12 @@
pars=pars,
Nmif=0L,
particles=particles,
- alg.pars=alg.pars,
+ 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
)
@@ -236,6 +271,7 @@
pars, ivps = character(0),
particles,
rw.sd, alg.pars,
+ Np, ic.lag, var.factor, cooling.factor,
weighted = TRUE, tol = 1e-17, warn = TRUE, max.fail = 0,
verbose = FALSE)
{
@@ -273,12 +309,13 @@
call.=FALSE
)
}
- if (missing(alg.pars))
- stop("mif error: ",sQuote("alg.pars")," must be specified",call.=FALSE)
-
+ if (missing(alg.pars)) alg.pars <- NULL # alg.pars is now deprecated
+
mif.internal(object,Nmif=Nmif,start=start,pars=pars,ivps=ivps,particles=particles,
- rw.sd=rw.sd,alg.pars=alg.pars,weighted=weighted,tol=tol,warn=warn,
- max.fail=max.fail,verbose=verbose,.ndone=0)
+ rw.sd=rw.sd,alg.pars=alg.pars,Np=Np,cooling.factor=cooling.factor,
+ var.factor=var.factor,ic.lag=ic.lag,
+ weighted=weighted,tol=tol,warn=warn,max.fail=max.fail,
+ verbose=verbose,.ndone=0)
}
)
Modified: pkg/data/ou2.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2009-04-13 14:27:26 UTC (rev 92)
+++ pkg/inst/ChangeLog 2009-04-14 14:43:04 UTC (rev 93)
@@ -1,5 +1,8 @@
2009-04-13 kingaa
+ * [r92] DESCRIPTION, data/ou2.rda, inst/ChangeLog,
+ inst/doc/compiled_code_in_pomp.pdf, inst/doc/intro_to_pomp.pdf,
+ tests/ou2-mif.Rout.save:
* [r91] inst/ChangeLog:
* [r90] R/mif.R, man/mif.Rd, tests/ou2-mif.R: rewrite mif methods
to use 'mif.internal'.
Modified: pkg/inst/doc/compiled_code_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.Rnw
===================================================================
--- pkg/inst/doc/intro_to_pomp.Rnw 2009-04-13 14:27:26 UTC (rev 92)
+++ pkg/inst/doc/intro_to_pomp.Rnw 2009-04-14 14:43:04 UTC (rev 93)
@@ -350,17 +350,9 @@
In this vignette, we'll use the default (multivariate normal) particle distribution.
Let's jump right in and run MIF to maximize the likelihood over two of the parameters and both initial conditions.
-We'll use 1000 particles, an exponential cooling factor of 0.95, and a fixed-lag smoother with lag 10 for the initial conditions:
+We'll use 1000 particles, an exponential cooling factor of 0.95, and a fixed-lag smoother with lag 10 for the initial conditions.
+Just to make it interesting, we'll start far from the true parameter values.
<<>>=
-alg.pars <- list(
- Np=1000,
- var.factor=1,
- ic.lag=10,
- cooling.factor=0.95
- )
-@
-Just to make it interesting, we'll start far from the true parameter values:
-<<>>=
start.p <- true.p
start.p[c('x1.0','x2.0','alpha.1','alpha.4')] <- c(45,-60,0.8,0.9)
fit <- mif(ou2,Nmif=1,start=start.p,
@@ -369,7 +361,10 @@
x1.0=5,x2.0=5,
alpha.1=0.1,alpha.4=0.1
),
- alg.pars=alg.pars,
+ Np=1000,
+ var.factor=1,
+ ic.lag=10,
+ cooling.factor=0.95,
max.fail=100
)
fit <- continue(fit,Nmif=79,max.fail=100)
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/man/mif.Rd
===================================================================
--- pkg/man/mif.Rd 2009-04-13 14:27:26 UTC (rev 92)
+++ pkg/man/mif.Rd 2009-04-14 14:43:04 UTC (rev 93)
@@ -15,7 +15,7 @@
\usage{
mif(object, \dots)
\S4method{mif}{pomp}(object, Nmif = 1, start, pars, ivps = character(0),
- particles, rw.sd, alg.pars, weighted = TRUE,
+ particles, rw.sd, Np, cooling.factor, var.factor, ic.lag, alg.pars, weighted = TRUE,
tol = 1e-17, warn = TRUE, max.fail = 0, verbose = FALSE)
\S4method{mif}{mif}(object, Nmif, \dots)
\S4method{continue}{mif}(object, Nmif = 1, \dots)
@@ -55,27 +55,26 @@
\code{rw.sd} must be non-negative (zeros are simply ignored),
the name of every positive element of \code{rw.sd} must be in either \code{pars} or \code{ivps}.
}
+ \item{Np}{
+ a positive integer;
+ the number of particles to use in filtering
+ }
+ \item{ic.lag}{
+ a positive integer;
+ the timepoint for fixed-lag smoothing of initial-value parameters
+ }
+ \item{var.factor}{
+ a positive number;
+ the scaling coefficient relating the width of the initial particle distribution to \code{rw.sd}
+ }
+ \item{cooling.factor}{
+ a positive number not greater than 1;
+ the exponential cooling factor, \code{alpha}.
+ }
\item{alg.pars}{
- A named list of algorithm parameters.
- This consists of
- \describe{
- \item{Np}{
- a positive integer;
- the number of particles to use in filtering
- }
- \item{ic.lag}{
- a positive integer;
- the timepoint for fixed-lag smoothing of initial-value parameters
- }
- \item{var.factor}{
- a positive number;
- the scaling coefficient relating the width of the initial particle distribution to \code{rw.sd}
- }
- \item{cooling.factor}{
- a positive number not greater than 1;
- the exponential cooling factor, \code{alpha}.
- }
- }
+ optional; a named list of the algorithm parameters \code{Np}, \code{cooling.factor}, \code{var.factor}, \code{ic.lag}.
+ The use of \code{alg.pars} is now deprecated and generates a warning.
+ It will be removed in the near future.
}
\item{weighted}{
Should a weighted average be used?
Modified: pkg/tests/ou2-mif.R
===================================================================
--- pkg/tests/ou2-mif.R 2009-04-13 14:27:26 UTC (rev 92)
+++ pkg/tests/ou2-mif.R 2009-04-14 14:43:04 UTC (rev 93)
@@ -28,12 +28,10 @@
x1.0=5,x2.0=5,
alpha.1=0.1,alpha.4=0.1
),
- alg.pars=list(
- Np=1000,
- var.factor=1,
- ic.lag=10,
- cooling.factor=0.95
- ),
+ Np=1000,
+ var.factor=1,
+ ic.lag=10,
+ cooling.factor=0.95,
max.fail=100
)
mif.fit <- continue(mif.fit,Nmif=70,max.fail=100)
@@ -59,7 +57,7 @@
pars=c("alpha.1","alpha.4"),
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),
- alg.pars=list(Np=100,cooling.factor=0.95,ic.lag=10,var.factor=1)
+ Np=100,cooling.factor=0.95,ic.lag=10,var.factor=1
)
fit <- mif(
ou2,
@@ -74,7 +72,7 @@
Nmif=2,
ivps=c("x1.0","x2.0"),
rw.sd=c(x1.0=5,x2.0=5,alpha.1=0.1,alpha.4=0.2),
- alg.pars=list(Np=1000,cooling.factor=0.95,ic.lag=10,var.factor=1)
+ Np=1000,cooling.factor=0.95,ic.lag=10,var.factor=1
)
fit <- continue(fit,Nmif=40)
ff <- pfilter(fit,pred.mean=T,filter.mean=T,pred.var=T,max.fail=100)
@@ -84,6 +82,6 @@
s <- coef(fit)
s[2] <- 0.01
fit <- mif(fit,Nmif=10,start=s)
-fit <- mif(ou2,Nmif=3,rw.sd=c(alpha.1=0.1,alpha.4=0.1),alg.pars=list(Np=1000,cooling.factor=0.98,var.factor=1,ic.lag=2))
-fit <- continue(fit,Nmif=5,alg.pars=list(Np=2000,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.factor=0.98,var.factor=1,ic.lag=2)
+fit <- continue(fit,Nmif=5,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=3)
Modified: pkg/tests/ou2-mif.Rout.save
===================================================================
--- pkg/tests/ou2-mif.Rout.save 2009-04-13 14:27:26 UTC (rev 92)
+++ pkg/tests/ou2-mif.Rout.save 2009-04-14 14:43:04 UTC (rev 93)
@@ -58,18 +58,16 @@
+ x1.0=5,x2.0=5,
+ alpha.1=0.1,alpha.4=0.1
+ ),
-+ alg.pars=list(
-+ Np=1000,
-+ var.factor=1,
-+ ic.lag=10,
-+ cooling.factor=0.95
-+ ),
++ Np=1000,
++ var.factor=1,
++ ic.lag=10,
++ cooling.factor=0.95,
+ max.fail=100
+ )
> mif.fit <- continue(mif.fit,Nmif=70,max.fail=100)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 1.054774 mins
+Time difference of 58.27099 secs
> cat("PF estimated log likelihood at MIF MLE\n")
PF estimated log likelihood at MIF MLE
> print(pfilter(mif.fit)$loglik,digits=4)
@@ -98,7 +96,7 @@
+ pars=c("alpha.1","alpha.4"),
+ 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),
-+ alg.pars=list(Np=100,cooling.factor=0.95,ic.lag=10,var.factor=1)
++ Np=100,cooling.factor=0.95,ic.lag=10,var.factor=1
+ )
> fit <- mif(
+ ou2,
@@ -108,12 +106,14 @@
+ rw.sd=c(x1.0=5,x2.0=5,alpha.1=0.1,alpha.4=0.2,alpha.3=0),
+ alg.pars=list(Np=1000,cooling.factor=0.95,ic.lag=10,var.factor=1)
+ )
+Warning message:
+mif warning: use of ‘alg.pars’ is deprecated; see the documentation for ‘mif’.
> fit <- mif(
+ ou2,
+ Nmif=2,
+ ivps=c("x1.0","x2.0"),
+ rw.sd=c(x1.0=5,x2.0=5,alpha.1=0.1,alpha.4=0.2),
-+ alg.pars=list(Np=1000,cooling.factor=0.95,ic.lag=10,var.factor=1)
++ Np=1000,cooling.factor=0.95,ic.lag=10,var.factor=1
+ )
> fit <- continue(fit,Nmif=40)
> ff <- pfilter(fit,pred.mean=T,filter.mean=T,pred.var=T,max.fail=100)
@@ -126,8 +126,8 @@
> s <- coef(fit)
> s[2] <- 0.01
> fit <- mif(fit,Nmif=10,start=s)
-> fit <- mif(ou2,Nmif=3,rw.sd=c(alpha.1=0.1,alpha.4=0.1),alg.pars=list(Np=1000,cooling.factor=0.98,var.factor=1,ic.lag=2))
-> fit <- continue(fit,Nmif=5,alg.pars=list(Np=2000,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.factor=0.98,var.factor=1,ic.lag=2)
+> fit <- continue(fit,Nmif=5,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=3)
Warning message:
mif warning: the variable(s) x2.0 have positive random-walk SDs specified, but are included in neither ‘pars’ nor ‘ivps’. These random walk SDs are ignored.
Modified: pkg/tests/sir.Rout.save
===================================================================
--- pkg/tests/sir.Rout.save 2009-04-13 14:27:26 UTC (rev 92)
+++ pkg/tests/sir.Rout.save 2009-04-14 14:43:04 UTC (rev 93)
@@ -1,5 +1,5 @@
-R version 2.7.1 (2008-06-23)
+R version 2.8.1 (2008-12-22)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
@@ -16,7 +16,8 @@
Type 'q()' to quit R.
> library(pomp)
-Loading required package: odesolve
+Loading required package: deSolve
+Loading required package: subplex
>
> tbasis <- seq(0,25,by=1/52)
> basis <- periodic.bspline.basis(tbasis,nbasis=3)
@@ -147,7 +148,7 @@
> x <- simulate(po,params=log(params),nsim=3)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 3.935402 secs
+Time difference of 2.595228 secs
>
> pdf(file='sir.pdf')
>
@@ -164,7 +165,7 @@
> X3 <- trajectory(po,params=log(params),times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 7.444703 secs
+Time difference of 5.991744 secs
> plot(t3,X3['I',1,],type='l')
>
> f1 <- dprocess(
@@ -255,7 +256,7 @@
> x <- simulate(po,params=log(params),nsim=3)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 0.3896399 secs
+Time difference of 0.2616198 secs
> plot(x[[1]],variables=c("S","I","R","cases","W"))
>
> t3 <- seq(0,20,by=1/52)
@@ -263,7 +264,7 @@
> X4 <- trajectory(po,params=log(params),times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 5.439281 secs
+Time difference of 4.824738 secs
> plot(t3,X4['I',1,],type='l')
>
> f2 <- dprocess(
@@ -314,7 +315,7 @@
R -30746 -26601 -22900 -16156.5 -6714 5662 22453
>
> print(max(abs(f2-f1),na.rm=T),digits=4)
-[1] 7.994e-15
+[1] 8.882e-15
> print(max(abs(g2-g1),na.rm=T),digits=4)
[1] 0
> print(max(abs(h2-h1),na.rm=T),digits=4)
More information about the pomp-commits
mailing list