[Pomp-commits] r153 - in pkg: . R inst inst/doc man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jun 24 17:19:13 CEST 2009
Author: kingaa
Date: 2009-06-24 17:19:12 +0200 (Wed, 24 Jun 2009)
New Revision: 153
Modified:
pkg/DESCRIPTION
pkg/R/mif.R
pkg/R/pomp-methods.R
pkg/inst/ChangeLog
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/intro_to_pomp.pdf
pkg/man/mif.Rd
pkg/tests/logistic.Rout.save
pkg/tests/ou2-kalman.Rout.save
pkg/tests/ou2-mif.R
pkg/tests/ou2-mif.Rout.save
pkg/tests/ou2-nlf.R
pkg/tests/ou2-nlf.Rout.save
pkg/tests/ou2-simulate.Rout.save
pkg/tests/rw2.R
pkg/tests/rw2.Rout.save
pkg/tests/sir.R
pkg/tests/sir.Rout.save
Log:
better testing
some more error checking in 'mif'
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2009-06-23 19:14:54 UTC (rev 152)
+++ pkg/DESCRIPTION 2009-06-24 15:19:12 UTC (rev 153)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.25-1
-Date: 2009-06-23
+Version: 0.25-2
+Date: 2009-06-24
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-06-23 19:14:54 UTC (rev 152)
+++ pkg/R/mif.R 2009-06-24 15:19:12 UTC (rev 153)
@@ -108,7 +108,11 @@
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)
-
+
+ Nmif <- as.integer(Nmif)
+ if (Nmif<0)
+ stop("mif error: ",sQuote("Nmif")," must be a positive integer",call.=FALSE)
+
if (verbose) {
cat("performing",Nmif,"MIF iteration(s) to estimate parameter(s)",
paste(pars,collapse=", "))
Modified: pkg/R/pomp-methods.R
===================================================================
--- pkg/R/pomp-methods.R 2009-06-23 19:14:54 UTC (rev 152)
+++ pkg/R/pomp-methods.R 2009-06-24 15:19:12 UTC (rev 153)
@@ -121,7 +121,7 @@
'pomp',
function (object) {
print(object)
- cat("zero time, t0 = ",object at t0,"\n")
+ cat("zero time, t0 = ",object at t0,"\n",sep="")
if (length(coef(object))>0) {
cat("parameter(s):\n")
print(coef(object))
@@ -133,9 +133,23 @@
cat("process model density, dprocess = \n")
print(object at dprocess)
cat("measurement model simulator, rmeasure = \n")
- print(object at rmeasure)
+ if (object at rmeasure@use==2) {
+ cat("native function, '",object at rmeasure@native.fun,"'",sep="")
+ if (length(object at rmeasure@PACKAGE)>0)
+ cat(", PACKAGE = '",object at rmeasure@PACKAGE,"'",sep="")
+ cat ("\n")
+ } else {
+ print(object at rmeasure@R.fun)
+ }
cat("measurement model density, dmeasure = \n")
- print(object at dmeasure)
+ if (object at dmeasure@use==2) {
+ cat("native function, '",object at dmeasure@native.fun,"'",sep="")
+ if (length(object at dmeasure@PACKAGE)>0)
+ cat(", PACKAGE = '",object at dmeasure@PACKAGE,"'",sep="")
+ cat ("\n")
+ } else {
+ print(object at dmeasure@R.fun)
+ }
cat("initializer = \n")
print(object at initializer)
cat("userdata = \n")
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2009-06-23 19:14:54 UTC (rev 152)
+++ pkg/inst/ChangeLog 2009-06-24 15:19:12 UTC (rev 153)
@@ -1,3 +1,11 @@
+2009-06-23 kingaa
+
+ * [r152] DESCRIPTION, R/mif.R, data/euler.sir.rda, data/ou2.rda,
+ data/rw2.rda, inst/ChangeLog,
+ inst/doc/advanced_topics_in_pomp.pdf, inst/doc/intro_to_pomp.pdf,
+ man/mif.Rd, tests/ou2-mif.R, tests/ou2-mif.Rout.save: alg.pars
+ has been removed as an argument to 'mif'
+
2009-06-19 kingaa
* [r151] DESCRIPTION, inst/ChangeLog,
Modified: pkg/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/man/mif.Rd
===================================================================
--- pkg/man/mif.Rd 2009-06-23 19:14:54 UTC (rev 152)
+++ pkg/man/mif.Rd 2009-06-24 15:19:12 UTC (rev 153)
@@ -98,14 +98,11 @@
}
\section{Re-running MIF Iterations}{
To re-run a sequence of MIF iterations, one can use the \code{mif} method on a \code{mif} object.
- The call sequence is \code{mif(object)}.
By default, the same parameters used for the original MIF run are re-used (except for \code{weighted}, \code{tol}, \code{warn}, \code{max.fail}, and \code{verbose}, the defaults of which are shown above).
If one does specify additional arguments, these will override the defaults.
}
\section{Continuing MIF Iterations}{
- One can continue a series of MIF iterations from where one left off.
- The call sequence is \code{continue(object, Nmif)}.
- This will perform \code{Nmif} additional MIF iterations on the \code{mif} object \code{object}.
+ One can continue a series of MIF iterations from where one left off using the \code{continue} method.
A call to \code{mif} to perform \code{Nmif=m} iterations followed by a call to \code{continue} to perform \code{Nmif=n} iterations will produce precisely the same effect as a single call to \code{mif} to perform \code{Nmif=m+n} iterations.
By default, all the algorithmic parameters are the same as used in the original call to \code{mif}.
Additional arguments will override the defaults.
Modified: pkg/tests/logistic.Rout.save
===================================================================
--- pkg/tests/logistic.Rout.save 2009-06-23 19:14:54 UTC (rev 152)
+++ pkg/tests/logistic.Rout.save 2009-06-24 15:19:12 UTC (rev 153)
@@ -1,6 +1,6 @@
-R version 2.10.0 Under development (unstable) (2009-05-07 r48492)
-Copyright (C) 2009 The R Foundation for Statistical Computing
+R version 2.8.1 (2008-12-22)
+Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -18,6 +18,7 @@
> library(pomp)
Loading required package: deSolve
Loading required package: subplex
+Loading required package: mvtnorm
>
> po <- pomp(
+ data=rbind(obs=rep(0,1000)),
Modified: pkg/tests/ou2-kalman.Rout.save
===================================================================
--- pkg/tests/ou2-kalman.Rout.save 2009-06-23 19:14:54 UTC (rev 152)
+++ pkg/tests/ou2-kalman.Rout.save 2009-06-24 15:19:12 UTC (rev 153)
@@ -18,6 +18,7 @@
> library(pomp)
Loading required package: deSolve
Loading required package: subplex
+Loading required package: mvtnorm
>
> kalman.filter <- function (y, x0, a, b, sigma, tau) {
+ n <- nrow(y)
@@ -103,7 +104,7 @@
> kalm.fit1 <- optim(p.guess,kalman,object=ou2,params=p.truth,hessian=T)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 39.45996 secs
+Time difference of 37.90149 secs
> tic <- Sys.time()
> print(-kalm.fit1$value,digits=4)
[1] -411.1
Modified: pkg/tests/ou2-mif.R
===================================================================
--- pkg/tests/ou2-mif.R 2009-06-23 19:14:54 UTC (rev 152)
+++ pkg/tests/ou2-mif.R 2009-06-24 15:19:12 UTC (rev 153)
@@ -2,7 +2,7 @@
data(ou2)
-set.seed(64857673)
+set.seed(64857673L)
fit1.pfilter <- pfilter(ou2,Np=1000)
cat("coefficients at `truth'\n")
@@ -11,46 +11,107 @@
print(fit1.pfilter$loglik,digits=4)
p.truth <- coef(ou2)
-guess <- p.truth
-guess[c('x1.0','x2.0','alpha.1','alpha.4')] <- c(45,-60,0.8,0.9)
+guess2 <- guess1 <- p.truth
+guess1[c('x1.0','x2.0','alpha.1','alpha.4')] <- c(45,-60,0.8,0.9)
+guess2[c('x1.0','x2.0','alpha.1','alpha.4')] <- c(60,-45,0.9,0.8)
-fit2.pfilter <- pfilter(ou2,params=guess,Np=1000,max.fail=1000,warn=F)
-cat("coefficients at guess\n")
-print(guess[c('x1.0','x2.0','alpha.1','alpha.4')],digits=4)
-cat("particle filter log likelihood at guess\n")
-print(fit2.pfilter$loglik,digits=4)
+mif1 <- mif(ou2,Nmif=10,start=guess1,
+ 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.1
+ ),
+ Np=1000,
+ var.factor=1,
+ ic.lag=10,
+ cooling.factor=0.95,
+ max.fail=100
+ )
-cat("running MIF\n")
-tic <- Sys.time()
-mif.fit <- mif(ou2,Nmif=10,start=guess,
- 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.1
- ),
- 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)
-cat("PF estimated log likelihood at MIF MLE\n")
-print(pfilter(mif.fit)$loglik,digits=4)
+mif2 <- mif(ou2,Nmif=10,start=guess2,
+ 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.1
+ ),
+ Np=1000,
+ var.factor=1,
+ ic.lag=10,
+ cooling.factor=0.95,
+ max.fail=100
+ )
-cat("coefficients at truth\n")
-print(coef(ou2,c('x1.0','x2.0','alpha.1','alpha.4')),digits=4)
-cat("MIF MLE\n")
-print(coef(mif.fit,c('x1.0','x2.0','alpha.1','alpha.4')),digits=4)
+pdf(file="ou2-mif.pdf")
+plot(mif1)
+compare.mif(mif2)
+try(compare.mif(mif1,mif2))
+compare.mif(list(mif1,mif2))
+dev.off()
-plot(mif.fit)
-compare.mif(mif.fit)
-compare.mif(list(mif.fit,mif.fit))
-
set.seed(33848585L)
+try(
+ mif(
+ ou2,
+ Nmif=1,
+ 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
+ )
+ )
+
+try(
+ mif(
+ ou2,
+ Nmif=1,
+ 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
+ )
+ )
+
+try(
+ mif(
+ ou2,
+ 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
+ )
+ )
+
+try(
+ mif(
+ ou2,
+ Nmif=1,
+ 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,ic.lag=10,var.factor=1
+ )
+ )
+
+try(
+ mif(
+ ou2,
+ 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
+ )
+ )
+
+try(
+ mif(
+ ou2,
+ 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
+ )
+ )
+
fit <- mif(
ou2,
Nmif=0,
@@ -60,20 +121,20 @@
Np=100,cooling.factor=0.95,ic.lag=10,var.factor=1
)
fit <- mif(
- ou2,
+ fit,
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),
Np=1000,cooling.factor=0.95,ic.lag=10,var.factor=1
)
-fit <- continue(fit,Nmif=40)
+fit <- continue(fit)
+fit <- continue(fit,Nmif=2)
ff <- pfilter(fit,pred.mean=T,filter.mean=T,pred.var=T,max.fail=100)
-print(ff$loglik)
fit <- mif(fit,rw.sd=c(x1.0=5,x2.0=5,alpha.1=0.1,alpha.4=0.1))
fit <- continue(fit,Nmif=2,ivps=c("x1.0"),pars=c("alpha.1"))
s <- coef(fit)
s[2] <- 0.01
-fit <- mif(fit,Nmif=10,start=s)
+fit <- mif(fit,Nmif=3,start=s)
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)
+fit <- continue(fit,Nmif=2,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=2)
Modified: pkg/tests/ou2-mif.Rout.save
===================================================================
--- pkg/tests/ou2-mif.Rout.save 2009-06-23 19:14:54 UTC (rev 152)
+++ pkg/tests/ou2-mif.Rout.save 2009-06-24 15:19:12 UTC (rev 153)
@@ -22,7 +22,7 @@
>
> data(ou2)
>
-> set.seed(64857673)
+> set.seed(64857673L)
>
> fit1.pfilter <- pfilter(ou2,Np=1000)
> cat("coefficients at `truth'\n")
@@ -36,61 +36,117 @@
[1] -412.9
>
> p.truth <- coef(ou2)
-> guess <- p.truth
-> guess[c('x1.0','x2.0','alpha.1','alpha.4')] <- c(45,-60,0.8,0.9)
+> guess2 <- guess1 <- p.truth
+> guess1[c('x1.0','x2.0','alpha.1','alpha.4')] <- c(45,-60,0.8,0.9)
+> guess2[c('x1.0','x2.0','alpha.1','alpha.4')] <- c(60,-45,0.9,0.8)
>
-> fit2.pfilter <- pfilter(ou2,params=guess,Np=1000,max.fail=1000,warn=F)
-> cat("coefficients at guess\n")
-coefficients at guess
-> print(guess[c('x1.0','x2.0','alpha.1','alpha.4')],digits=4)
- x1.0 x2.0 alpha.1 alpha.4
- 45.0 -60.0 0.8 0.9
-> cat("particle filter log likelihood at guess\n")
-particle filter log likelihood at guess
-> print(fit2.pfilter$loglik,digits=4)
-[1] -3913
+> mif1 <- mif(ou2,Nmif=10,start=guess1,
++ 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.1
++ ),
++ Np=1000,
++ var.factor=1,
++ ic.lag=10,
++ cooling.factor=0.95,
++ max.fail=100
++ )
>
-> cat("running MIF\n")
-running MIF
-> tic <- Sys.time()
-> mif.fit <- mif(ou2,Nmif=10,start=guess,
-+ 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.1
-+ ),
-+ 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 56.82802 secs
-> cat("PF estimated log likelihood at MIF MLE\n")
-PF estimated log likelihood at MIF MLE
-> print(pfilter(mif.fit)$loglik,digits=4)
-[1] -411.2
+> mif2 <- mif(ou2,Nmif=10,start=guess2,
++ 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.1
++ ),
++ Np=1000,
++ var.factor=1,
++ ic.lag=10,
++ cooling.factor=0.95,
++ max.fail=100
++ )
>
-> cat("coefficients at truth\n")
-coefficients at truth
-> print(coef(ou2,c('x1.0','x2.0','alpha.1','alpha.4')),digits=4)
- x1.0 x2.0 alpha.1 alpha.4
- 50.00 -50.00 0.90 0.99
-> cat("MIF MLE\n")
-MIF MLE
-> print(coef(mif.fit,c('x1.0','x2.0','alpha.1','alpha.4')),digits=4)
- x1.0 x2.0 alpha.1 alpha.4
- 50.8435 -48.7910 0.9035 0.9902
+> pdf(file="ou2-mif.pdf")
+> plot(mif1)
+> compare.mif(mif2)
+> try(compare.mif(mif1,mif2))
+Error in compare.mif(mif1, mif2) :
+ unused argument(s) (<S4 object of class "mif">)
+> compare.mif(list(mif1,mif2))
+> dev.off()
+null device
+ 1
>
-> plot(mif.fit)
-> compare.mif(mif.fit)
-> compare.mif(list(mif.fit,mif.fit))
->
> set.seed(33848585L)
>
+> try(
++ mif(
++ ou2,
++ Nmif=1,
++ 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
++ )
++ )
+Error : mif error: ‘pars’ and ‘ivps’ must be mutually disjoint subsets of ‘names(start)’ and must have a positive random-walk SDs specified in ‘rw.sd’
+>
+> try(
++ mif(
++ ou2,
++ Nmif=1,
++ 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
++ )
++ )
+Error : mif error: ‘pars’ and ‘ivps’ must be mutually disjoint subsets of ‘names(start)’ and must have a positive random-walk SDs specified in ‘rw.sd’
+>
+> try(
++ mif(
++ ou2,
++ 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
++ )
++ )
+Error : mif error: ‘pars’ and ‘ivps’ must be mutually disjoint subsets of ‘names(start)’ and must have a positive random-walk SDs specified in ‘rw.sd’
+>
+> try(
++ mif(
++ ou2,
++ Nmif=1,
++ 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,ic.lag=10,var.factor=1
++ )
++ )
+Error : mif error: ‘cooling.factor’ must be specified
+>
+> try(
++ mif(
++ ou2,
++ 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
++ )
++ )
+Error : mif error: ‘Np’ must be a positive integer
+>
+> try(
++ mif(
++ ou2,
++ 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
++ )
++ )
+Error : mif error: ‘Nmif’ must be a positive integer
+>
> fit <- mif(
+ ou2,
+ Nmif=0,
@@ -100,26 +156,34 @@
+ Np=100,cooling.factor=0.95,ic.lag=10,var.factor=1
+ )
> fit <- mif(
-+ ou2,
++ fit,
+ 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),
+ Np=1000,cooling.factor=0.95,ic.lag=10,var.factor=1
+ )
-> fit <- continue(fit,Nmif=40)
+> fit <- continue(fit)
+> fit <- continue(fit,Nmif=2)
> ff <- pfilter(fit,pred.mean=T,filter.mean=T,pred.var=T,max.fail=100)
-> print(ff$loglik)
-[1] -495.7194
+filtering failure at time t = 3
+filtering failure at time t = 4
+filtering failure at time t = 5
+filtering failure at time t = 6
+filtering failure at time t = 7
+filtering failure at time t = 8
+filtering failure at time t = 9
+filtering failure at time t = 10
+filtering failure at time t = 11
> fit <- mif(fit,rw.sd=c(x1.0=5,x2.0=5,alpha.1=0.1,alpha.4=0.1))
> fit <- continue(fit,Nmif=2,ivps=c("x1.0"),pars=c("alpha.1"))
Warning message:
mif warning: the variable(s) alpha.4, x2.0 have positive random-walk SDs specified, but are included in neither ‘pars’ nor ‘ivps’. These random walk SDs are ignored.
> s <- coef(fit)
> s[2] <- 0.01
-> fit <- mif(fit,Nmif=10,start=s)
+> fit <- mif(fit,Nmif=3,start=s)
> 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)
+> fit <- continue(fit,Nmif=2,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=2)
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/ou2-nlf.R
===================================================================
--- pkg/tests/ou2-nlf.R 2009-06-23 19:14:54 UTC (rev 152)
+++ pkg/tests/ou2-nlf.R 2009-06-24 15:19:12 UTC (rev 153)
@@ -2,14 +2,12 @@
data(ou2)
-set.seed(1066)
+set.seed(1066L)
po <- ou2
coef(po,c("x1.0","x2.0","alpha.1","alpha.4")) <- c(0,0,0.1,0.2)
-po <- simulate(po,times=(1:10000))
-p.truth <- coef(po)
-guess <- p.truth
-## guess[c('x1.0','x2.0','alpha.1','alpha.4')] <- c(45,-60,0.8,0.9)
+po <- simulate(po)
+guess <- p.truth <- coef(po)
m1 <- nlf(
object=po,
@@ -17,10 +15,35 @@
est=c("alpha.1","alpha.4"),
lags=c(1,2),
nconverge=100,
- nasymp=2500,
+ nasymp=1000,
method="Nelder-Mead",
maxit=500,
trace=1,
verbose=TRUE,
lql.frac = 0.025
)
+
+se <- p.truth
+se[] <- NA
+se[names(m1$se)] <- m1$se
+print(cbind(truth=p.truth,fit=m1$params,se=se),digits=3)
+
+po <- simulate(po,times=(1:1000))
+m2 <- nlf(
+ object=po,
+ start=guess,
+ est=c("alpha.1","alpha.4"),
+ lags=c(1,2),
+ nconverge=100,
+ nasymp=10000,
+ method="Nelder-Mead",
+ maxit=500,
+ trace=1,
+ verbose=TRUE,
+ lql.frac = 0.025
+ )
+
+se <- p.truth
+se[] <- NA
+se[names(m2$se)] <- m2$se
+print(cbind(truth=p.truth,fit=m2$params,se=se),digits=3)
Modified: pkg/tests/ou2-nlf.Rout.save
===================================================================
--- pkg/tests/ou2-nlf.Rout.save 2009-06-23 19:14:54 UTC (rev 152)
+++ pkg/tests/ou2-nlf.Rout.save 2009-06-24 15:19:12 UTC (rev 153)
@@ -1,6 +1,6 @@
-R version 2.10.0 Under development (unstable) (2009-05-07 r48492)
-Copyright (C) 2009 The R Foundation for Statistical Computing
+R version 2.8.1 (2008-12-22)
+Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -18,17 +18,16 @@
> library(pomp)
Loading required package: deSolve
Loading required package: subplex
+Loading required package: mvtnorm
>
> data(ou2)
>
-> set.seed(1066)
+> set.seed(1066L)
>
> po <- ou2
> coef(po,c("x1.0","x2.0","alpha.1","alpha.4")) <- c(0,0,0.1,0.2)
-> po <- simulate(po,times=(1:10000))
-> p.truth <- coef(po)
-> guess <- p.truth
-> ## guess[c('x1.0','x2.0','alpha.1','alpha.4')] <- c(45,-60,0.8,0.9)
+> po <- simulate(po)
+> guess <- p.truth <- coef(po)
>
> m1 <- nlf(
+ object=po,
@@ -36,7 +35,7 @@
+ est=c("alpha.1","alpha.4"),
+ lags=c(1,2),
+ nconverge=100,
-+ nasymp=2500,
++ nasymp=1000,
+ method="Nelder-Mead",
+ maxit=500,
+ trace=1,
@@ -44,26 +43,114 @@
+ lql.frac = 0.025
+ )
Nelder-Mead direct search function minimizer
-function value for initial parameters = 40059.679966
- Scaled convergence tolerance is 0.000596936
+function value for initial parameters = 390.960996
+ Scaled convergence tolerance is 5.82577e-06
Stepsize computed as 0.020000
-BUILD 3 40059.679966 40059.013186
-REFLECTION 5 40059.024962 40058.355658
-HI-REDUCTION 7 40059.013186 40058.355658
-HI-REDUCTION 9 40058.480960 40058.355658
-REFLECTION 11 40058.428912 40058.287812
-HI-REDUCTION 13 40058.355658 40058.283261
-HI-REDUCTION 15 40058.287812 40058.270086
-HI-REDUCTION 17 40058.283261 40058.262895
-HI-REDUCTION 19 40058.270086 40058.262895
-HI-REDUCTION 21 40058.263735 40058.260296
-HI-REDUCTION 23 40058.262895 40058.260296
+BUILD 3 391.170142 390.960996
+EXTENSION 5 391.121379 390.848054
+EXTENSION 7 390.960996 390.495673
+LO-REDUCTION 9 390.848054 390.495673
+EXTENSION 11 390.509255 390.032148
+EXTENSION 13 390.495673 389.677144
+REFLECTION 15 390.032148 389.601518
+EXTENSION 17 389.677144 389.069152
+LO-REDUCTION 19 389.601518 389.069152
+LO-REDUCTION 21 389.316744 389.069152
+REFLECTION 23 389.097526 388.964593
+LO-REDUCTION 25 389.069152 388.964593
+LO-REDUCTION 27 388.977691 388.959018
+HI-REDUCTION 29 388.964593 388.955938
+LO-REDUCTION 31 388.959018 388.954600
+HI-REDUCTION 33 388.955938 388.954600
+HI-REDUCTION 35 388.955009 388.953977
+LO-REDUCTION 37 388.954600 388.953960
+HI-REDUCTION 39 388.954009 388.953960
+HI-REDUCTION 41 388.953977 388.953875
+HI-REDUCTION 43 388.953960 388.953849
+LO-REDUCTION 45 388.953875 388.953849
+HI-REDUCTION 47 388.953868 388.953849
Exiting from Nelder Mead minimizer
- 25 function evaluations used
+ 49 function evaluations used
h in NLF = 0.1
-epsilon in NLF = 0.4420236 0.2750263
-Fitted param 1 -4.042142 -4.042142 up in 'nlf'
-Fitted param 1 -4.031832 down in 'NLF'
-Fitted param 2 -4.032517 -4.032517 up in 'nlf'
-Fitted param 2 -4.03145 down in 'NLF'
+epsilon in NLF = 0.3120962 0.3004542
+Fitted param 1 -3.984644 -3.984644 up in ‘nlf’
+Fitted param 1 -3.985982 down in ‘NLF’
+Fitted param 2 -3.993283 -3.993283 up in ‘nlf’
+Fitted param 2 -3.993728 down in ‘NLF’
>
+> se <- p.truth
+> se[] <- NA
+> se[names(m1$se)] <- m1$se
+> print(cbind(truth=p.truth,fit=m1$params,se=se),digits=3)
+ truth fit se
+alpha.1 0.1 -0.1539 0.114
+alpha.2 0.0 0.0000 NA
+alpha.3 0.0 0.0000 NA
+alpha.4 0.2 0.0159 0.110
+sigma.1 1.0 1.0000 NA
+sigma.2 0.0 0.0000 NA
+sigma.3 2.0 2.0000 NA
+tau 1.0 1.0000 NA
+x1.0 0.0 0.0000 NA
+x2.0 0.0 0.0000 NA
+>
+> po <- simulate(po,times=(1:1000))
+> m2 <- nlf(
++ object=po,
++ start=guess,
++ est=c("alpha.1","alpha.4"),
++ lags=c(1,2),
++ nconverge=100,
++ nasymp=10000,
++ method="Nelder-Mead",
++ maxit=500,
++ trace=1,
++ verbose=TRUE,
++ lql.frac = 0.025
++ )
+ Nelder-Mead direct search function minimizer
+function value for initial parameters = 4012.121680
+ Scaled convergence tolerance is 5.97853e-05
+Stepsize computed as 0.020000
+BUILD 3 4012.528145 4012.121680
+REFLECTION 5 4012.384913 4012.067214
+EXTENSION 7 4012.121680 4011.649828
+EXTENSION 9 4012.067214 4011.551849
+REFLECTION 11 4011.649828 4011.355226
+HI-REDUCTION 13 4011.551849 4011.355226
+HI-REDUCTION 15 4011.436440 4011.355226
+LO-REDUCTION 17 4011.403589 4011.351028
+HI-REDUCTION 19 4011.355226 4011.351028
+HI-REDUCTION 21 4011.353532 4011.344389
+HI-REDUCTION 23 4011.351028 4011.344389
+LO-REDUCTION 25 4011.346527 4011.344130
+HI-REDUCTION 27 4011.344389 4011.344130
+HI-REDUCTION 29 4011.344238 4011.343759
+HI-REDUCTION 31 4011.344130 4011.343759
+HI-REDUCTION 33 4011.343867 4011.343759
+LO-REDUCTION 35 4011.343836 4011.343759
+Exiting from Nelder Mead minimizer
+ 37 function evaluations used
+h in NLF = 0.1
+epsilon in NLF = 0.4842402 0.2705555
+Fitted param 1 -4.054405 -4.054405 up in ‘nlf’
+Fitted param 1 -4.052725 down in ‘NLF’
+Fitted param 2 -4.049163 -4.049163 up in ‘nlf’
+Fitted param 2 -4.048027 down in ‘NLF’
+>
+> se <- p.truth
+> se[] <- NA
+> se[names(m2$se)] <- m2$se
+> print(cbind(truth=p.truth,fit=m2$params,se=se),digits=3)
+ truth fit se
+alpha.1 0.1 0.0169 0.0677
+alpha.2 0.0 0.0000 NA
+alpha.3 0.0 0.0000 NA
+alpha.4 0.2 0.1907 0.0490
+sigma.1 1.0 1.0000 NA
+sigma.2 0.0 0.0000 NA
+sigma.3 2.0 2.0000 NA
+tau 1.0 1.0000 NA
+x1.0 0.0 0.0000 NA
+x2.0 0.0 0.0000 NA
+>
Modified: pkg/tests/ou2-simulate.Rout.save
===================================================================
--- pkg/tests/ou2-simulate.Rout.save 2009-06-23 19:14:54 UTC (rev 152)
+++ pkg/tests/ou2-simulate.Rout.save 2009-06-24 15:19:12 UTC (rev 153)
@@ -1,6 +1,6 @@
-R version 2.10.0 Under development (unstable) (2009-05-07 r48492)
-Copyright (C) 2009 The R Foundation for Statistical Computing
+R version 2.8.1 (2008-12-22)
+Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -18,6 +18,7 @@
> library(pomp)
Loading required package: deSolve
Loading required package: subplex
+Loading required package: mvtnorm
>
> data(ou2)
>
@@ -33,7 +34,7 @@
> ou2.sim <- simulate(ou2,params=p,nsim=100,seed=32043858)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 0.3552971 secs
+Time difference of 0.3518229 secs
>
> coef(ou2,c('x1.0','x2.0')) <- c(-50,50)
>
Modified: pkg/tests/rw2.R
===================================================================
--- pkg/tests/rw2.R 2009-06-23 19:14:54 UTC (rev 152)
+++ pkg/tests/rw2.R 2009-06-24 15:19:12 UTC (rev 153)
@@ -58,6 +58,15 @@
if (log) f else exp(f)
}
+bad.initializer <- function (params, t0, ...)
+{
+ ivpnames <- c("x1.0","x2.0")
+ x <- params[ivpnames]
+ x
+}
+
+p <- rbind(s1=c(2,2,3),s2=c(0.1,1,2),tau=c(1,5,0),x1.0=c(0,0,5),x2.0=c(0,0,0))
+
rw2 <- pomp(
rprocess = rw.rprocess,
dprocess = rw.dprocess,
@@ -65,6 +74,7 @@
y1 ~ norm(mean=x1,sd=tau),
y2 ~ norm(mean=x2,sd=tau)
),
+ initializer=bad.initializer,
times=1:100,
data=rbind(
y1=rep(0,100),
@@ -74,7 +84,28 @@
useless=23
)
-p <- rbind(s1=c(2,2,3),s2=c(0.1,1,2),tau=c(1,5,0),x1.0=c(0,0,5),x2.0=c(0,0,0))
+show(rw2)
+
+try(
+ simulate(rw2,params=p)
+ )
+
+rw2 <- pomp(
+ rprocess = rw.rprocess,
+ dprocess = rw.dprocess,
+ measurement.model=list(
+ y1 ~ norm(mean=x1,sd=tau),
+ y2 ~ norm(mean=x2,sd=tau)
+ ),
+ times=1:100,
+ data=rbind(
+ y1=rep(0,100),
+ y2=rep(0,100)
+ ),
+ t0=0,
+ useless=23
+ )
+
examples <- simulate(rw2,params=p)
rw2 <- examples[[1]]
Modified: pkg/tests/rw2.Rout.save
===================================================================
--- pkg/tests/rw2.Rout.save 2009-06-23 19:14:54 UTC (rev 152)
+++ pkg/tests/rw2.Rout.save 2009-06-24 15:19:12 UTC (rev 153)
@@ -18,6 +18,7 @@
> library(pomp)
Loading required package: deSolve
Loading required package: subplex
+Loading required package: mvtnorm
>
> set.seed(45768683)
>
@@ -77,6 +78,15 @@
+ if (log) f else exp(f)
+ }
>
+> bad.initializer <- function (params, t0, ...)
++ {
++ ivpnames <- c("x1.0","x2.0")
++ x <- params[ivpnames]
++ x
++ }
+>
+> p <- rbind(s1=c(2,2,3),s2=c(0.1,1,2),tau=c(1,5,0),x1.0=c(0,0,5),x2.0=c(0,0,0))
+>
> rw2 <- pomp(
+ rprocess = rw.rprocess,
+ dprocess = rw.dprocess,
@@ -84,6 +94,7 @@
+ y1 ~ norm(mean=x1,sd=tau),
+ y2 ~ norm(mean=x2,sd=tau)
+ ),
++ initializer=bad.initializer,
+ times=1:100,
+ data=rbind(
+ y1=rep(0,100),
@@ -93,7 +104,202 @@
+ useless=23
+ )
>
-> p <- rbind(s1=c(2,2,3),s2=c(0.1,1,2),tau=c(1,5,0),x1.0=c(0,0,5),x2.0=c(0,0,0))
+> show(rw2)
+ time y1 y2
+1 1 0 0
+2 2 0 0
+3 3 0 0
+4 4 0 0
+5 5 0 0
+6 6 0 0
+7 7 0 0
+8 8 0 0
+9 9 0 0
+10 10 0 0
+11 11 0 0
+12 12 0 0
+13 13 0 0
+14 14 0 0
+15 15 0 0
+16 16 0 0
+17 17 0 0
+18 18 0 0
+19 19 0 0
+20 20 0 0
+21 21 0 0
+22 22 0 0
+23 23 0 0
+24 24 0 0
+25 25 0 0
+26 26 0 0
+27 27 0 0
+28 28 0 0
+29 29 0 0
+30 30 0 0
+31 31 0 0
+32 32 0 0
+33 33 0 0
+34 34 0 0
+35 35 0 0
+36 36 0 0
+37 37 0 0
+38 38 0 0
+39 39 0 0
+40 40 0 0
+41 41 0 0
+42 42 0 0
+43 43 0 0
+44 44 0 0
+45 45 0 0
+46 46 0 0
+47 47 0 0
+48 48 0 0
+49 49 0 0
+50 50 0 0
+51 51 0 0
+52 52 0 0
+53 53 0 0
+54 54 0 0
+55 55 0 0
+56 56 0 0
+57 57 0 0
+58 58 0 0
+59 59 0 0
+60 60 0 0
+61 61 0 0
+62 62 0 0
+63 63 0 0
+64 64 0 0
+65 65 0 0
+66 66 0 0
+67 67 0 0
+68 68 0 0
+69 69 0 0
+70 70 0 0
+71 71 0 0
+72 72 0 0
+73 73 0 0
+74 74 0 0
+75 75 0 0
+76 76 0 0
+77 77 0 0
+78 78 0 0
+79 79 0 0
+80 80 0 0
+81 81 0 0
+82 82 0 0
+83 83 0 0
+84 84 0 0
+85 85 0 0
+86 86 0 0
+87 87 0 0
+88 88 0 0
+89 89 0 0
+90 90 0 0
+91 91 0 0
+92 92 0 0
+93 93 0 0
+94 94 0 0
+95 95 0 0
+96 96 0 0
+97 97 0 0
+98 98 0 0
+99 99 0 0
+100 100 0 0
+zero time, t0 = 0
+parameter(s) unspecified
+process model simulator, rprocess =
+function (params, xstart, times, ...)
+{
+ nvars <- nrow(xstart)
+ nreps <- ncol(params)
+ ntimes <- length(times)
+ dt <- diff(times)
+ x <- array(0, dim = c(nvars, nreps, ntimes))
+ rownames(x) <- rownames(xstart)
+ noise.sds <- params[c("s1", "s2"), ]
+ x[, , 1] <- xstart
+ for (j in 2:ntimes) {
+ x[c("x1", "x2"), , j] <- rnorm(n = 2 * nreps, mean = x[c("x1",
+ "x2"), , j - 1], sd = noise.sds * dt[j - 1])
+ }
+ x
+}
+process model density, dprocess =
+function (x, times, params, log = FALSE, ...)
+{
+ nreps <- ncol(params)
+ ntimes <- length(times)
+ dt <- diff(times)
+ d <- array(0, dim = c(2, nreps, ntimes - 1))
+ noise.sds <- params[c("s1", "s2"), ]
+ for (j in 2:ntimes) d[, , j - 1] <- dnorm(x[, , j] - x[,
+ , j - 1], mean = 0, sd = noise.sds * dt[j - 1], log = TRUE)
+ d <- apply(d, c(2, 3), sum)
+ if (log)
+ d
+ else exp(d)
+}
+measurement model simulator, rmeasure =
+function (x, t, params, covars, ...)
+{
+ y <- numeric(length = nobs)
+ names(y) <- obsnames
+ for (k in 1:nobs) {
+ y[k] <- eval(rcalls[[k]], envir = as.list(c(x, params,
+ covars, t = t)))
+ }
+ y
+}
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 153
More information about the pomp-commits
mailing list