[Pomp-commits] r658 - in pkg: . data inst inst/data-R inst/doc tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Apr 14 23:12:52 CEST 2012
Author: kingaa
Date: 2012-04-14 23:12:52 +0200 (Sat, 14 Apr 2012)
New Revision: 658
Modified:
pkg/DESCRIPTION
pkg/data/bbs.rda
pkg/data/blowflies.rda
pkg/data/dacca.rda
pkg/data/euler.sir.rda
pkg/data/gillespie.sir.rda
pkg/data/gompertz.rda
pkg/data/ou2.rda
pkg/data/ricker.rda
pkg/data/rw2.rda
pkg/data/verhulst.rda
pkg/inst/NEWS
pkg/inst/data-R/gompertz.R
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/bsmc-ricker-flat-prior.rda
pkg/inst/doc/gompertz-multi-mif.rda
pkg/inst/doc/gompertz-trajmatch.rda
pkg/inst/doc/intro_to_pomp.Rnw
pkg/inst/doc/intro_to_pomp.pdf
pkg/inst/doc/nlf-block-boot.rda
pkg/inst/doc/nlf-boot.rda
pkg/inst/doc/nlf-fit-from-truth.rda
pkg/inst/doc/nlf-fits.rda
pkg/inst/doc/nlf-lag-tests.rda
pkg/inst/doc/nlf-multi-short.rda
pkg/inst/doc/ricker-mif.rda
pkg/inst/doc/ricker-probe-match.rda
pkg/tests/gompertz.R
pkg/tests/gompertz.Rout.save
Log:
- remove name-changes from Gompertz parameter transforms
- bring vignettes up to date
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2012-04-13 12:00:39 UTC (rev 657)
+++ pkg/DESCRIPTION 2012-04-14 21:12:52 UTC (rev 658)
@@ -2,7 +2,7 @@
Type: Package
Title: Statistical inference for partially observed Markov processes
Version: 0.41-5
-Date: 2012-04-12
+Date: 2012-04-14
Revision: $Rev$
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>
Modified: pkg/data/bbs.rda
===================================================================
(Binary files differ)
Modified: pkg/data/blowflies.rda
===================================================================
(Binary files differ)
Modified: pkg/data/dacca.rda
===================================================================
(Binary files differ)
Modified: pkg/data/euler.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/gillespie.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/gompertz.rda
===================================================================
(Binary files differ)
Modified: pkg/data/ou2.rda
===================================================================
(Binary files differ)
Modified: pkg/data/ricker.rda
===================================================================
(Binary files differ)
Modified: pkg/data/rw2.rda
===================================================================
(Binary files differ)
Modified: pkg/data/verhulst.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/NEWS
===================================================================
--- pkg/inst/NEWS 2012-04-13 12:00:39 UTC (rev 657)
+++ pkg/inst/NEWS 2012-04-14 21:12:52 UTC (rev 658)
@@ -2,6 +2,9 @@
0.41-5
o An experimental facility for constructing pomp objects with native C routines is now included.
+ o The 'gompertz' example parameter transformations have been changed.
+ No longer are the names of the parameter vector changed in the transformation.
+
0.41-4
o The 'blowflies', 'euler.sir', 'gillespie.sir', 'bbs', and 'ricker' data()-loadable examples have been changed to make the parameterization simpler and more natural.
Although these changes are not backward compatible, they make for much simpler explanations.
@@ -13,7 +16,7 @@
o Demos have been created to show some examples.
0.41-2
- o Fix bug in 'bbs' example.
+ o A segfault bug in 'bbs' example has been fixed.
o The posterior medians (not means) are now stored in the 'params' slot of the 'bsmcd.pomp' object.
Modified: pkg/inst/data-R/gompertz.R
===================================================================
--- pkg/inst/data-R/gompertz.R 2012-04-13 12:00:39 UTC (rev 657)
+++ pkg/inst/data-R/gompertz.R 2012-04-14 21:12:52 UTC (rev 658)
@@ -15,14 +15,10 @@
statenames=c("X"),
obsnames=c("Y"),
parameter.transform=function(params,...){
- params <- exp(params[c("log.X.0","log.r","log.K","log.tau","log.sigma")])
- names(params) <- c("X.0","r","K","tau","sigma")
- params
+ exp(params)
},
parameter.inv.transform=function(params,...){
- params <- log(params[c("X.0","r","K","tau","sigma")])
- names(params) <- c("log.X.0","log.r","log.K","log.tau","log.sigma")
- params
+ log(params)
}
)
Modified: pkg/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/bsmc-ricker-flat-prior.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/gompertz-multi-mif.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/gompertz-trajmatch.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.Rnw
===================================================================
--- pkg/inst/doc/intro_to_pomp.Rnw 2012-04-13 12:00:39 UTC (rev 657)
+++ pkg/inst/doc/intro_to_pomp.Rnw 2012-04-14 21:12:52 UTC (rev 658)
@@ -159,7 +159,7 @@
As we noted above, for this particular model, it isn't necessary to use plug-and-play methods:
one can obtain exact maximum likelihood estimates of this model's parameters using the Kalman filter.
We will demonstrate this below and use it to check the results of plug-and-play inference.
-For now, let's approach this model as we would a more complex model for which no such exact estimation is available.
+In this document, we'll approach this model as we would a more complex model for which no such exact estimation is available.
\section{Defining a partially observed Markov process in \pkg{pomp}.}
@@ -520,16 +520,10 @@
gompertz <- pomp(
gompertz,
parameter.transform=function(params,...){
- params <- exp(params[c("log.r","log.K","log.tau",
- "log.sigma","log.X.0")])
- names(params) <- c("r","K","tau","sigma","X.0")
- params
+ exp(params)
},
parameter.inv.transform=function(params,...){
- params <- log(params[c("r","K","tau","sigma","X.0")])
- names(params) <- c("log.r","log.K","log.tau",
- "log.sigma","log.X.0")
- params
+ log(params)
}
)
@
@@ -549,14 +543,17 @@
@
On the other hand, given parameter values on the estimation scale, one can set the parameters by
<<eval=T>>=
-coef(gompertz,transform=TRUE) <- c(log.r=log(0.1),log.K=0,
- log.tau=log(0.1),log.sigma=log(0.1),log.X.0=0)
+coef(gompertz,transform=TRUE) <- c(r=log(0.1),K=0,tau=log(0.1),
+ sigma=log(0.1),X.0=0)
@
and get them by
<<eval=T>>=
coef(gompertz,transform=TRUE)
@
-
+We can check that the parameters were set correctly on the natural scale:
+<<eval=T>>=
+coef(gompertz)
+@
To be clear: the parameter-setting expression \code{coef(gompertz) <- theta} assumes that \code{theta} is a named vector of parameters on the natural scale, while \code{coef(gompertz,transform=TRUE) <- theta} assumes that \code{theta} is on the estimation scale.
Note that it's the user's responsibility to ensure that the \code{parameter.transform} and \code{parameter.inv.transform} are actually mutually inverse.
@@ -576,10 +573,13 @@
For your convenience, codes creating this \code{pomp} object are included with the package.
To view them, do:
<<eval=F>>=
-file.show(system.file("examples/gompertz.R",package="pomp"))
+file.show(system.file("demo/gompertz.R",package="pomp"))
file.show(system.file("examples/gompertz.c",package="pomp"))
@
-
+or run the demo
+<<eval=F>>=
+demo(gompertz)
+@
\clearpage
\section{Estimating parameters using iterated filtering: \code{mif}}
@@ -625,9 +625,9 @@
Nmif=100,
start=theta.guess,
transform=TRUE,
- pars=c("log.r","log.sigma","log.tau"),
+ pars=estpars,
rw.sd=c(
- log.r=0.02,log.sigma=0.02,log.tau=0.05
+ r=0.02,sigma=0.02,tau=0.05
),
Np=2000,
var.factor=4,
@@ -691,12 +691,12 @@
<<gompertz-post-mif>>
@
-<<multi-mif-plot,echo=F,eval=F>>=
+<<multi-mif-plot,echo=F,eval=T>>=
op <- par(mfrow=c(4,1),mar=c(3,3,0,0),mgp=c(2,1,0),bty='l')
loglik <- sapply(mf,function(x)conv.rec(x,"loglik"))
-log.r <- sapply(mf,function(x)conv.rec(x,"log.r"))
-log.sigma <- sapply(mf,function(x)conv.rec(x,"log.sigma"))
-log.tau <- sapply(mf,function(x)conv.rec(x,"log.tau"))
+log.r <- sapply(mf,function(x)conv.rec(x,"r"))
+log.sigma <- sapply(mf,function(x)conv.rec(x,"sigma"))
+log.tau <- sapply(mf,function(x)conv.rec(x,"tau"))
matplot(loglik,type='l',lty=1,xlab="",ylab=expression(log~L),xaxt='n',ylim=max(loglik,na.rm=T)+c(-12,3))
matplot(log.r,type='l',lty=1,xlab="",ylab=expression(log~r),xaxt='n')
matplot(log.sigma,type='l',lty=1,xlab="",ylab=expression(log~sigma),xaxt='n')
@@ -779,7 +779,7 @@
new.gompertz,
start=coef(new.gompertz),
transform=TRUE,
- est=c("log.r","log.K","log.tau","log.X.0"),
+ est=c("r","K","tau","X.0"),
method="Nelder-Mead",
maxit=1000,
reltol=1e-8
@@ -798,7 +798,7 @@
Fig.~\ref{fig:trajmatch-plot} shows the results of this fit.
\begin{figure}
-<<trajmatch-plot,fig=T,echo=F,eval=T>>=
+<<trajmatch-plot,fig=T,echo=F,eval=T,height=4,width=6>>=
op <- par(mfrow=c(1,1),mar=c(3,3,0,0),mgp=c(2,1,0),bty='l')
plot(time(tm),obs(tm,"Y"),xlab="time",ylab=expression(X,Y),type='o')
lines(time(tm),states(tm,"X"),lwd=2)
@@ -1743,11 +1743,7 @@
b <- params[c("b1","b2","b3")] # contact-rate coefficients
beta <- b%*%covars # flexible seasonality
beta.sd <- params["beta.sd"] # extrademographic noise intensity
- beta.var <- beta.sd*beta.sd
- if (beta.var > 0)
- dW <- rgamma(n=1,shape=delta.t/beta.var,scale=beta.var)
- else
- dW <- delta.t
+ dW <- rgammawn(n=1,sigma=beta.sd,dt=delta.t) # Gamma white noise
foi <- (beta*x["I"]/N+iota)*dW/delta.t # the force of infection
trans <- c(
rpois(n=1,lambda=mu*N*delta.t), # births are assumed to be Poisson
@@ -1804,16 +1800,16 @@
lty=1,col=c("red","blue","orange"),xlab="time (yr)",
ylab=quote("basis functions"~list(s[1],s[2],s[3])))
bb <- coef(complex.sir,c("b1","b2","b3"))
-plot.window(c(0,2),c(0,1)*exp(max(bb)))
-lines(tbasis,exp(basis%*%bb),col="black",lwd=3,lty=1)
-lines(tbasis,exp(basis%*%bb),col="white",lwd=2.5,lty="23")
+plot.window(c(0,2),c(0,1)*max(bb))
+lines(tbasis,basis%*%bb,col="black",lwd=3,lty=1)
+lines(tbasis,basis%*%bb,col="white",lwd=2.5,lty="23")
axis(side=4)
mtext(
side=4,
line=2,
text=bquote(
- beta==exp(.(b1)*s[1]+.(b2)*s[2]+.(b3)*s[3]),
- where=as.list(exp(coef(complex.sir)))
+ beta==.(b1)*s[1]+.(b2)*s[2]+.(b3)*s[3],
+ where=as.list(coef(complex.sir))
)
)
par(op)
@@ -1821,7 +1817,7 @@
\caption{
Periodic B-spline basis functions can be used to construct flexible periodic functions.
The colored lines show the three basis functions, $s_1$, $s_2$, $s_3$.
- The dashed black line shows the seasonality $\beta(t)$ assumed in \code{complex.sir}.
+ The dashed black line shows the seasonal transmission $\beta(t)$ assumed in \code{complex.sir}.
}
\label{fig:seas-basis-plot}
\end{figure}
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/nlf-block-boot.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/nlf-boot.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/nlf-fit-from-truth.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/nlf-fits.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/nlf-lag-tests.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/nlf-multi-short.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/ricker-mif.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/ricker-probe-match.rda
===================================================================
(Binary files differ)
Modified: pkg/tests/gompertz.R
===================================================================
--- pkg/tests/gompertz.R 2012-04-13 12:00:39 UTC (rev 657)
+++ pkg/tests/gompertz.R 2012-04-14 21:12:52 UTC (rev 658)
@@ -7,7 +7,7 @@
coef(po)
coef(po,transform=TRUE)
coef(po,c("r","X.0"))
-coef(po,c("log.r","log.X.0"),transform=TRUE)
+coef(po,c("r","X.0"),transform=TRUE)
coef(po,c("r","K")) <- c(0.2,1)
coef(po)
coef(po,transform=TRUE)
@@ -18,15 +18,15 @@
Nmif=5,Np=1000,
transform=TRUE,
ic.lag=1,var.factor=1,cooling.factor=0.99,
- rw.sd=c(log.r=0.02,log.K=0.02)
+ rw.sd=c(r=0.02,K=0.02)
)
coef(mf,transform=TRUE)
coef(mf)
conv.rec(mf)
conv.rec(mf,transform=TRUE)
-conv.rec(mf,c("loglik","log.r"))
+conv.rec(mf,c("loglik","r"))
try(conv.rec(mf,c("loglik","r"),transform=FALSE))
-try(conv.rec(mf,c("loglik","log.r"),transform=TRUE))
+try(conv.rec(mf,c("loglik","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)
Modified: pkg/tests/gompertz.Rout.save
===================================================================
--- pkg/tests/gompertz.Rout.save 2012-04-13 12:00:39 UTC (rev 657)
+++ pkg/tests/gompertz.Rout.save 2012-04-14 21:12:52 UTC (rev 658)
@@ -29,21 +29,21 @@
K r sigma tau X.0
1.0 0.1 0.1 0.1 1.0
> coef(po,transform=TRUE)
- log.X.0 log.r log.K log.tau log.sigma
- 0.000 -2.303 0.000 -2.303 -2.303
+ K r sigma tau X.0
+ 0.000 -2.303 -2.303 -2.303 0.000
> coef(po,c("r","X.0"))
r X.0
0.1 1.0
-> coef(po,c("log.r","log.X.0"),transform=TRUE)
- log.r log.X.0
- -2.303 0.000
+> coef(po,c("r","X.0"),transform=TRUE)
+ r X.0
+-2.303 0.000
> coef(po,c("r","K")) <- c(0.2,1)
> coef(po)
K r sigma tau X.0
1.0 0.2 0.1 0.1 1.0
> coef(po,transform=TRUE)
- log.X.0 log.r log.K log.tau log.sigma
- 0.000 -1.609 0.000 -2.303 -2.303
+ K r sigma tau X.0
+ 0.000 -1.609 -2.303 -2.303 0.000
>
> set.seed(93848585L)
> mf <- mif(
@@ -51,57 +51,69 @@
+ Nmif=5,Np=1000,
+ transform=TRUE,
+ ic.lag=1,var.factor=1,cooling.factor=0.99,
-+ rw.sd=c(log.r=0.02,log.K=0.02)
++ rw.sd=c(r=0.02,K=0.02)
+ )
> coef(mf,transform=TRUE)
- log.X.0 log.r log.K log.tau log.sigma
- 0.0000 -1.6195 0.0447 -2.3026 -2.3026
+ K r sigma tau X.0
+ 0.04671 -1.61336 -2.30259 -2.30259 0.00000
> coef(mf)
- X.0 r K tau sigma
-1.000 0.198 1.046 0.100 0.100
+ K r sigma tau X.0
+1.0478 0.1992 0.1000 0.1000 1.0000
> conv.rec(mf)
- loglik nfail log.X.0 log.r log.K log.tau log.sigma
-0 30.55 0 0 -1.609 0.000000 -2.303 -2.303
-1 31.18 0 0 -1.613 0.008437 -2.303 -2.303
-2 30.67 0 0 -1.616 0.017788 -2.303 -2.303
-3 30.55 0 0 -1.615 0.026324 -2.303 -2.303
-4 31.06 0 0 -1.612 0.039947 -2.303 -2.303
-5 NA NA 0 -1.619 0.044705 -2.303 -2.303
+ loglik nfail K r sigma tau X.0
+0 30.35 0 0.00000 -1.609 -2.303 -2.303 0
+1 30.59 0 0.01048 -1.609 -2.303 -2.303 0
+2 30.96 0 0.02378 -1.610 -2.303 -2.303 0
+3 29.78 0 0.03161 -1.607 -2.303 -2.303 0
+4 31.18 0 0.03785 -1.611 -2.303 -2.303 0
+5 NA NA 0.04671 -1.613 -2.303 -2.303 0
> 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
+ K r sigma tau X.0 loglik nfail
+0 1.000 0.2000 0.1 0.1 1 30.35 0
+1 1.011 0.2000 0.1 0.1 1 30.59 0
+2 1.024 0.1999 0.1 0.1 1 30.96 0
+3 1.032 0.2005 0.1 0.1 1 29.78 0
+4 1.039 0.1996 0.1 0.1 1 31.18 0
+5 1.048 0.1992 0.1 0.1 1 NA NA
+> conv.rec(mf,c("loglik","r"))
+ loglik r
+0 30.35 -1.609
+1 30.59 -1.609
+2 30.96 -1.610
+3 29.78 -1.607
+4 31.18 -1.611
+5 NA -1.613
> 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)'
+ loglik r
+0 30.35 -1.609
+1 30.59 -1.609
+2 30.96 -1.610
+3 29.78 -1.607
+4 31.18 -1.611
+5 NA -1.613
+> try(conv.rec(mf,c("loglik","r"),transform=TRUE))
+ loglik r
+0 30.35 0.2000
+1 30.59 0.2000
+2 30.96 0.1999
+3 29.78 0.2005
+4 31.18 0.1996
+5 NA 0.1992
> 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
+0 30.35 0.2000
+1 30.59 0.2000
+2 30.96 0.1999
+3 29.78 0.2005
+4 31.18 0.1996
+5 NA 0.1992
> conv.rec(mf,c("loglik"),transform=TRUE)
0 1 2 3 4 5
-30.55 31.18 30.67 30.55 31.06 NA
+30.35 30.59 30.96 29.78 31.18 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
+1.000 1.011 1.024 1.032 1.039 1.048
>
> proc.time()
user system elapsed
- 3.476 0.032 3.527
+ 1.448 0.048 1.509
More information about the pomp-commits
mailing list