[Pomp-commits] r530 - in pkg: . R inst/doc man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Aug 2 19:28:12 CEST 2011
Author: kingaa
Date: 2011-08-02 19:28:12 +0200 (Tue, 02 Aug 2011)
New Revision: 530
Added:
pkg/tests/gompertz.R
pkg/tests/gompertz.Rout.save
Modified:
pkg/DESCRIPTION
pkg/R/mif-methods.R
pkg/R/pomp-methods.R
pkg/inst/doc/intro_to_pomp.Rnw
pkg/inst/doc/intro_to_pomp.pdf
pkg/man/mif-methods.Rd
Log:
- add 'transform' argument to 'conv.rec' method for 'mif' objects
- add new unit test 'tests/gompertz.R' which tests the transformation facilities in 'coef', 'coef<-', and 'conv.rec'
- small changes to 'intro_to_pomp' vignette
- bug fix in 'pomp:::pomp.transform'
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2011-07-28 21:34:42 UTC (rev 529)
+++ pkg/DESCRIPTION 2011-08-02 17:28:12 UTC (rev 530)
@@ -2,7 +2,7 @@
Type: Package
Title: Statistical inference for partially observed Markov processes
Version: 0.39-1
-Date: 2011-07-28
+Date: 2011-08-03
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/R/mif-methods.R
===================================================================
--- pkg/R/mif-methods.R 2011-07-28 21:34:42 UTC (rev 529)
+++ pkg/R/mif-methods.R 2011-08-02 17:28:12 UTC (rev 530)
@@ -11,9 +11,38 @@
setMethod(
'conv.rec',
'mif',
- function (object, pars, ...) {
- if (missing(pars)) pars <- colnames(object at conv.rec)
- object at conv.rec[,pars]
+ function (object, pars, transform = FALSE, ...) {
+ if (transform) {
+ pars.improper <- c("loglik","nfail")
+ pars.proper <- setdiff(colnames(object at conv.rec),pars.improper)
+ retval <- cbind(
+ t(
+ pomp.transform(
+ object,
+ params=t(object at conv.rec[,pars.proper]),
+ dir="inverse"
+ )
+ ),
+ object at conv.rec[,pars.improper]
+ )
+ } else {
+ retval <- object at conv.rec
+ }
+ if (missing(pars))
+ retval
+ else {
+ bad.pars <- setdiff(pars,colnames(retval))
+ if (length(bad.pars)>0)
+ stop(
+ "in ",sQuote("conv.rec"),": name(s) ",
+ paste(sQuote(bad.pars),collapse=","),
+ " correspond to no parameter(s) in ",
+ if (transform) sQuote("conv.rec(object,transform=TRUE)")
+ else sQuote("conv.rec(object,transform=FALSE)"),
+ call.=FALSE
+ )
+ retval[,pars]
+ }
}
)
Modified: pkg/R/pomp-methods.R
===================================================================
--- pkg/R/pomp-methods.R 2011-07-28 21:34:42 UTC (rev 529)
+++ pkg/R/pomp-methods.R 2011-08-02 17:28:12 UTC (rev 530)
@@ -170,11 +170,14 @@
forward=function (x) do.call(object at par.trans,c(list(x),object at userdata)),
inverse=function (x) do.call(object at par.untrans,c(list(x),object at userdata))
)
- if (r > 1)
+ if (r > 1) {
retval <- apply(params,2:r,tfunc)
- else
+ no.names <- is.null(rownames(retval))
+ } else {
retval <- tfunc(params)
- if (is.null(names(retval)))
+ no.names <- is.null(names(retval))
+ }
+ if (no.names)
switch(
dir,
forward=stop(
Modified: pkg/inst/doc/intro_to_pomp.Rnw
===================================================================
--- pkg/inst/doc/intro_to_pomp.Rnw 2011-07-28 21:34:42 UTC (rev 529)
+++ pkg/inst/doc/intro_to_pomp.Rnw 2011-08-02 17:28:12 UTC (rev 530)
@@ -592,7 +592,7 @@
)
@
-The parameter transformations come into play only via the \code{coef} and \code{coef<-} methods for getting and setting parameters.
+The parameter transformations come into play via the \code{coef} and \code{coef<-} methods for getting and setting parameters.
One can set the parameters by doing
<<>>=
coef(gompertz,transform=TRUE) <- theta
@@ -688,6 +688,7 @@
@
<<gompertz-post-mif,eval=F,echo=F>>=
+theta.true <- coef(gompertz)
theta.mif <- apply(sapply(mf,coef),1,mean)
loglik.mif <- replicate(n=10,logLik(pfilter(mf[[1]],params=theta.mif,Np=10000)))
bl <- mean(loglik.mif)
@@ -718,7 +719,7 @@
file=binary.file
)
}
-cbind(
+rbind(
# guess=c(signif(theta.guess[estpars],3),loglik=round(loglik.guess,1)),
mle=c(signif(theta.mif[estpars],3),loglik=round(loglik.mif.est,1),loglik.se=signif(loglik.mif.se,2)),
truth=c(signif(theta.true[estpars],3),loglik=round(loglik.true.est,1),loglik.se=signif(loglik.true.se,2))
@@ -728,6 +729,8 @@
Each of the \Sexpr{length(mf)} \code{mif} runs ends up at a different place.
In this case (but by no means in every case), we can average across these parameter estimates to get an approximate maximum likelihood estimate.
We'll evaluate the likelihood several times at this estimate to get an idea of the Monte Carlo error in our likelihood estimate.
+The particle filter produces an unbiased estimate of the likelihood;
+therefore, we'll average the likelihoods, not the log likelihoods.
<<eval=F>>=
<<gompertz-post-mif>>
@
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/man/mif-methods.Rd
===================================================================
--- pkg/man/mif-methods.Rd 2011-07-28 21:34:42 UTC (rev 529)
+++ pkg/man/mif-methods.Rd 2011-08-02 17:28:12 UTC (rev 530)
@@ -13,7 +13,7 @@
\description{Methods of the "mif" class.}
\usage{
\S4method{logLik}{mif}(object, \dots)
-\S4method{conv.rec}{mif}(object, pars, \dots)
+\S4method{conv.rec}{mif}(object, pars, transform = FALSE, \dots)
\S4method{plot}{mif}(x, y = NULL, \dots)
compare.mif(z)
}
@@ -23,6 +23,11 @@
\item{x}{The \code{mif} object.}
\item{y}{Ignored.}
\item{z}{A \code{mif} object or list of \code{mif} objects.}
+ \item{transform}{
+ optional logical;
+ should the parameter transformations be applied?
+ See \code{\link[=coef-pomp]{coef}} for details.
+ }
\item{\dots}{
Further arguments (either ignored or passed to underlying functions).
}
Added: pkg/tests/gompertz.R
===================================================================
--- pkg/tests/gompertz.R (rev 0)
+++ pkg/tests/gompertz.R 2011-08-02 17:28:12 UTC (rev 530)
@@ -0,0 +1,31 @@
+library(pomp)
+options(digits=4)
+
+data(gompertz)
+
+po <- gompertz
+coef(po)
+coef(po,transform=TRUE)
+coef(po,c("log.r","X.0"))
+coef(po,c("r","X.0"),transform=TRUE)
+coef(po,c("r","K"),transform=TRUE) <- c(0.2,1)
+coef(po)
+coef(po,transform=TRUE)
+
+set.seed(93848585L)
+mf <- mif(
+ po,
+ Nmif=5,Np=1000,
+ ic.lag=1,var.factor=1,cooling.factor=0.99,
+ rw.sd=c(log.r=0.02,log.K=0.02)
+ )
+coef(mf,transform=TRUE)
+conv.rec(mf)
+conv.rec(mf,transform=TRUE)
+conv.rec(mf,c("loglik","log.r"))
+try(conv.rec(mf,c("loglik","r"),transform=FALSE))
+try(conv.rec(mf,c("loglik","log.r"),transform=TRUE))
+conv.rec(mf,c("loglik","r"),transform=TRUE)
+conv.rec(mf,c("loglik"),transform=TRUE)
+conv.rec(mf,c("K"),transform=TRUE)
+
Added: pkg/tests/gompertz.Rout.save
===================================================================
--- pkg/tests/gompertz.Rout.save (rev 0)
+++ pkg/tests/gompertz.Rout.save 2011-08-02 17:28:12 UTC (rev 530)
@@ -0,0 +1,101 @@
+
+R version 2.13.1 (2011-07-08)
+Copyright (C) 2011 The R Foundation for Statistical Computing
+ISBN 3-900051-07-0
+Platform: x86_64-unknown-linux-gnu (64-bit)
+
+R is free software and comes with ABSOLUTELY NO WARRANTY.
+You are welcome to redistribute it under certain conditions.
+Type 'license()' or 'licence()' for distribution details.
+
+R is a collaborative project with many contributors.
+Type 'contributors()' for more information and
+'citation()' on how to cite R or R packages in publications.
+
+Type 'demo()' for some demos, 'help()' for on-line help, or
+'help.start()' for an HTML browser interface to help.
+Type 'q()' to quit R.
+
+> library(pomp)
+Loading required package: mvtnorm
+Loading required package: subplex
+Loading required package: deSolve
+> options(digits=4)
+>
+> data(gompertz)
+>
+> po <- gompertz
+> coef(po)
+ X.0 log.r log.K log.tau log.sigma
+ 1.000 -2.303 0.000 -2.303 -2.303
+> coef(po,transform=TRUE)
+ X.0 r K tau sigma
+ 1.0 0.1 1.0 0.1 0.1
+> coef(po,c("log.r","X.0"))
+ log.r X.0
+-2.303 1.000
+> coef(po,c("r","X.0"),transform=TRUE)
+ r X.0
+0.1 1.0
+> coef(po,c("r","K"),transform=TRUE) <- c(0.2,1)
+> coef(po)
+ X.0 log.r log.K log.tau log.sigma
+ 1.000 -1.609 0.000 -2.303 -2.303
+> coef(po,transform=TRUE)
+ X.0 r K tau sigma
+ 1.0 0.2 1.0 0.1 0.1
+>
+> set.seed(93848585L)
+> mf <- mif(
++ po,
++ Nmif=5,Np=1000,
++ ic.lag=1,var.factor=1,cooling.factor=0.99,
++ rw.sd=c(log.r=0.02,log.K=0.02)
++ )
+> coef(mf,transform=TRUE)
+ X.0 r K tau sigma
+1.000 0.198 1.046 0.100 0.100
+> conv.rec(mf)
+ loglik nfail X.0 log.r log.K log.tau log.sigma
+0 30.55 0 1 -1.609 0.000000 -2.303 -2.303
+1 31.18 0 1 -1.613 0.008437 -2.303 -2.303
+2 30.67 0 1 -1.616 0.017788 -2.303 -2.303
+3 30.55 0 1 -1.615 0.026324 -2.303 -2.303
+4 31.06 0 1 -1.612 0.039947 -2.303 -2.303
+5 NA NA 1 -1.619 0.044705 -2.303 -2.303
+> conv.rec(mf,transform=TRUE)
+ X.0 r K tau sigma loglik nfail
+0 1 0.2000 1.000 0.1 0.1 30.55 0
+1 1 0.1994 1.008 0.1 0.1 31.18 0
+2 1 0.1987 1.018 0.1 0.1 30.67 0
+3 1 0.1989 1.027 0.1 0.1 30.55 0
+4 1 0.1994 1.041 0.1 0.1 31.06 0
+5 1 0.1980 1.046 0.1 0.1 NA NA
+> conv.rec(mf,c("loglik","log.r"))
+ loglik log.r
+0 30.55 -1.609
+1 31.18 -1.613
+2 30.67 -1.616
+3 30.55 -1.615
+4 31.06 -1.612
+5 NA -1.619
+> try(conv.rec(mf,c("loglik","r"),transform=FALSE))
+Error : in 'conv.rec': name(s) 'r' correspond to no parameter(s) in 'conv.rec(object,transform=FALSE)'
+> try(conv.rec(mf,c("loglik","log.r"),transform=TRUE))
+Error : in 'conv.rec': name(s) 'log.r' correspond to no parameter(s) in 'conv.rec(object,transform=TRUE)'
+> conv.rec(mf,c("loglik","r"),transform=TRUE)
+ loglik r
+0 30.55 0.2000
+1 31.18 0.1994
+2 30.67 0.1987
+3 30.55 0.1989
+4 31.06 0.1994
+5 NA 0.1980
+> conv.rec(mf,c("loglik"),transform=TRUE)
+ 0 1 2 3 4 5
+30.55 31.18 30.67 30.55 31.06 NA
+> conv.rec(mf,c("K"),transform=TRUE)
+ 0 1 2 3 4 5
+1.000 1.008 1.018 1.027 1.041 1.046
+>
+>
More information about the pomp-commits
mailing list