[Pomp-commits] r982 - in pkg/pomp: R tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jun 25 19:39:12 CEST 2014
Author: kingaa
Date: 2014-06-25 19:39:11 +0200 (Wed, 25 Jun 2014)
New Revision: 982
Added:
pkg/pomp/tests/ou2-abc.R
pkg/pomp/tests/ou2-abc.Rout.save
Removed:
pkg/pomp/tests/abc.R
pkg/pomp/tests/abc.Rout.save
Modified:
pkg/pomp/R/abc-methods.R
pkg/pomp/R/pmcmc-methods.R
pkg/pomp/tests/ou2-pmcmc.R
pkg/pomp/tests/ou2-pmcmc.Rout.save
Log:
- fix bug in 'conv.rec' for 'abc' and 'pmcmc' lists
- rename 'tests/abc.R' -> 'tests/ou2-abc.R'
- add some tests for thinning and windowing 'conv.rec' output
Modified: pkg/pomp/R/abc-methods.R
===================================================================
--- pkg/pomp/R/abc-methods.R 2014-06-25 17:38:58 UTC (rev 981)
+++ pkg/pomp/R/abc-methods.R 2014-06-25 17:39:11 UTC (rev 982)
@@ -1,32 +1,5 @@
## this file contains short definitions of methods for the 'abc' class
-## extract the convergence record as an 'mcmc' object
-setMethod(
- 'conv.rec',
- 'abc',
- function (object, pars, ...) {
- if (missing(pars)) pars <- colnames(object at conv.rec)
- coda::mcmc(object at conv.rec[,pars,drop=FALSE],...)
- }
- )
-
-## plot abc object
-setMethod(
- "plot",
- "abc",
- function (x, y, pars, scatter = FALSE, ...) {
- ## if (missing(pars)) pars <- x at pars
- ## if (scatter) {
- ## pairs(as.matrix(conv.rec(x,pars)))
- ## } else {
- ## plot.ts(conv.rec(x,pars),xlab="iteration",...)
- ## }
- if (!missing(y))
- warning(sQuote("y")," is ignored")
- abc.diagnostics(c(x),pars=pars,scatter=scatter,...)
- }
- )
-
## abcList class
setClass(
'abcList',
@@ -100,21 +73,47 @@
}
)
+## extract the convergence record as an 'mcmc' object
+setMethod(
+ 'conv.rec',
+ 'abc',
+ function (object, pars, ...) {
+ if (missing(pars)) pars <- colnames(object at conv.rec)
+ coda::mcmc(object at conv.rec[,pars,drop=FALSE],...)
+ }
+ )
+
## extract the convergence record as an 'mcmc.list' object
setMethod(
'conv.rec',
signature=signature(object='abcList'),
definition=function (object, ...) {
- coda::mcmc.list(lapply(object,conv.rec,...))
+ f <- selectMethod("conv.rec","abc")
+ coda::mcmc.list(lapply(object,f,...))
}
)
+## plot abc object
setMethod(
"plot",
+ "abc",
+ function (x, y, pars, scatter = FALSE, ...) {
+ if (!missing(y)) {
+ y <- substitute(y)
+ warning(sQuote(y)," is ignored")
+ }
+ abc.diagnostics(c(x),pars=pars,scatter=scatter,...)
+ }
+ )
+
+setMethod(
+ "plot",
signature=signature(x='abcList'),
definition=function (x, y, ...) {
- if (!missing(y))
- warning(sQuote("y")," is ignored")
+ if (!missing(y)) {
+ y <- substitute(y)
+ warning(sQuote(y)," is ignored")
+ }
abc.diagnostics(x,...)
}
)
Modified: pkg/pomp/R/pmcmc-methods.R
===================================================================
--- pkg/pomp/R/pmcmc-methods.R 2014-06-25 17:38:58 UTC (rev 981)
+++ pkg/pomp/R/pmcmc-methods.R 2014-06-25 17:39:11 UTC (rev 982)
@@ -3,28 +3,6 @@
## extract the estimated log likelihood
setMethod('logLik','pmcmc',function(object,...)object at loglik)
-## extract the convergence record as a coda::mcmc object
-setMethod(
- 'conv.rec',
- signature=signature(object='pmcmc'),
- function (object, pars, ...) {
- if (missing(pars)) pars <- colnames(object at conv.rec)
- coda::mcmc(object at conv.rec[,pars,drop=FALSE],...)
- }
- )
-
-## plot pmcmc object
-setMethod(
- "plot",
- signature=signature(x='pmcmc'),
- function (x, y, ...) {
- if (!missing(y))
- warning(sQuote("y")," is ignored")
- pmcmc.diagnostics(list(x))
- }
- )
-
-
## pmcmcList class
setClass(
'pmcmcList',
@@ -98,21 +76,48 @@
}
)
+## extract the convergence record as a coda::mcmc object
+setMethod(
+ 'conv.rec',
+ signature=signature(object='pmcmc'),
+ function (object, pars, ...) {
+ if (missing(pars)) pars <- colnames(object at conv.rec)
+ coda::mcmc(object at conv.rec[,pars,drop=FALSE],...)
+ }
+ )
+
## extract the convergence records as a coda::mcmc.list object
setMethod(
'conv.rec',
signature=signature(object='pmcmcList'),
definition=function (object, ...) {
- coda::mcmc.list(lapply(object,conv.rec,...))
+ f <- selectMethod("conv.rec","pmcmc")
+ coda::mcmc.list(lapply(object,f,...))
}
)
+## plot pmcmc object
setMethod(
"plot",
+ signature=signature(x='pmcmc'),
+ function (x, y, ...) {
+ if (!missing(y)) {
+ y <- substitute(y)
+ warning(sQuote(y)," is ignored")
+ }
+ pmcmc.diagnostics(list(x))
+ }
+ )
+
+
+setMethod(
+ "plot",
signature=signature(x='pmcmcList'),
definition=function (x, y, ...) {
- if (!missing(y))
- warning(sQuote("y")," is ignored")
+ if (!missing(y)) {
+ y <- substitute(y)
+ warning(sQuote(y)," is ignored")
+ }
pmcmc.diagnostics(x)
}
)
Deleted: pkg/pomp/tests/abc.R
===================================================================
--- pkg/pomp/tests/abc.R 2014-06-25 17:38:58 UTC (rev 981)
+++ pkg/pomp/tests/abc.R 2014-06-25 17:39:11 UTC (rev 982)
@@ -1,103 +0,0 @@
-### OU2 test of abc for pomp
-
-library(pomp)
-pompExample(ou2)
-
-pdf(file='abc.pdf')
-
-set.seed(2079015564L)
-
-probes.good <- list(
- y1.mean=probe.mean(var="y1"),
- y2.mean=probe.mean(var="y2"),
- probe.acf(var="y1",lags=c(0,5)),
- probe.acf(var="y2",lags=c(0,5)),
- probe.ccf(vars=c("y1","y2"),lags=0)
- )
-psim <- probe(ou2,probes=probes.good,nsim=200)
-plot(psim)
-## why do simulations sometimes seem extreme with respect to these probes?
-
-scale.dat <- apply(psim$simvals,2,sd)
-
-po <- ou2
-
-abc1 <- abc(po,
- Nabc=10000,
- start=coef(ou2),
- pars=c("alpha.1","alpha.2"),
- probes=probes.good,
- scale=scale.dat,
- epsilon=1.7,
- rw.sd= c(alpha.1=0.01,alpha.2=0.01)
- )
-
-plot(abc1,scatter=TRUE)
-plot(abc1)
-
-## check how sticky the chain is:
-runs <- rle(as.vector(conv.rec(abc1)[, "alpha.1"]))
-hist(runs$lengths)
-mean(runs$length)
-
-abc2 <- abc(po,
- Nabc=2000,
- pars=c("alpha.1","alpha.2"),
- probes=probes.good,
- scale=scale.dat,
- epsilon=1,
- rw.sd= c(alpha.1=0.01,alpha.2=0.01)
- )
-plot(abc2)
-
-abc3 <- abc(po,
- Nabc=2000,
- probes=probes.good,
- scale=scale.dat,
- epsilon=2,
- rw.sd= c(alpha.1=0.01,alpha.2=0.01)
- )
-abc3 <- continue(abc3,Nabc=3000)
-plot(abc3)
-
-abc4 <- abc(probe(po,probes=probes.good,nsim=200),
- Nabc=2000,
- scale=scale.dat,
- epsilon=2,
- rw.sd= c(alpha.1=0.01,alpha.2=0.01)
- )
-plot(abc4)
-
-abc5 <- abc(abc4,Nabc=1000)
-plot(abc5)
-
-dprior6 <- function (params, log, ...) {
- ll <- sum(
- dnorm(
- x=params[c("alpha.1","alpha.2","alpha.3","alpha.4")],
- mean=c(0.8,-0.5,0.3,0.9),
- sd=5,
- log=TRUE
- )
- )
- if (log) ll else exp(ll)
-}
-
-abc6 <- abc(pomp(po,dprior=dprior6),
- Nabc=2000,
- pars=c("alpha.1","alpha.2"),
- probes=probes.good,
- scale=scale.dat,
- epsilon=1,
- rw.sd= c(alpha.1=0.01,alpha.2=0.01)
- )
-plot(abc6)
-
-try(abc7 <- c(abc2,abc3))
-plot(abc7 <- c(abc2,abc4))
-plot(abc7,scatter=TRUE)
-plot(conv.rec(c(abc2,abc4)))
-plot(conv.rec(c(abc7,abc6)))
-
-dev.off()
-
Deleted: pkg/pomp/tests/abc.Rout.save
===================================================================
--- pkg/pomp/tests/abc.Rout.save 2014-06-25 17:38:58 UTC (rev 981)
+++ pkg/pomp/tests/abc.Rout.save 2014-06-25 17:39:11 UTC (rev 982)
@@ -1,137 +0,0 @@
-
-R version 3.1.0 (2014-04-10) -- "Spring Dance"
-Copyright (C) 2014 The R Foundation for Statistical Computing
-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.
-
-> ### OU2 test of abc for pomp
->
-> library(pomp)
-Loading required package: mvtnorm
-Loading required package: subplex
-Loading required package: nloptr
-Loading required package: deSolve
-Loading required package: coda
-Loading required package: lattice
-> pompExample(ou2)
-newly created pomp object(s):
- ou2
->
-> pdf(file='abc.pdf')
->
-> set.seed(2079015564L)
->
-> probes.good <- list(
-+ y1.mean=probe.mean(var="y1"),
-+ y2.mean=probe.mean(var="y2"),
-+ probe.acf(var="y1",lags=c(0,5)),
-+ probe.acf(var="y2",lags=c(0,5)),
-+ probe.ccf(vars=c("y1","y2"),lags=0)
-+ )
-> psim <- probe(ou2,probes=probes.good,nsim=200)
-> plot(psim)
-> ## why do simulations sometimes seem extreme with respect to these probes?
->
-> scale.dat <- apply(psim$simvals,2,sd)
->
-> po <- ou2
->
-> abc1 <- abc(po,
-+ Nabc=10000,
-+ start=coef(ou2),
-+ pars=c("alpha.1","alpha.2"),
-+ probes=probes.good,
-+ scale=scale.dat,
-+ epsilon=1.7,
-+ rw.sd= c(alpha.1=0.01,alpha.2=0.01)
-+ )
->
-> plot(abc1,scatter=TRUE)
-> plot(abc1)
->
-> ## check how sticky the chain is:
-> runs <- rle(as.vector(conv.rec(abc1)[, "alpha.1"]))
-> hist(runs$lengths)
-> mean(runs$length)
-[1] 6.949965
->
-> abc2 <- abc(po,
-+ Nabc=2000,
-+ pars=c("alpha.1","alpha.2"),
-+ probes=probes.good,
-+ scale=scale.dat,
-+ epsilon=1,
-+ rw.sd= c(alpha.1=0.01,alpha.2=0.01)
-+ )
-> plot(abc2)
->
-> abc3 <- abc(po,
-+ Nabc=2000,
-+ probes=probes.good,
-+ scale=scale.dat,
-+ epsilon=2,
-+ rw.sd= c(alpha.1=0.01,alpha.2=0.01)
-+ )
-> abc3 <- continue(abc3,Nabc=3000)
-> plot(abc3)
->
-> abc4 <- abc(probe(po,probes=probes.good,nsim=200),
-+ Nabc=2000,
-+ scale=scale.dat,
-+ epsilon=2,
-+ rw.sd= c(alpha.1=0.01,alpha.2=0.01)
-+ )
-> plot(abc4)
->
-> abc5 <- abc(abc4,Nabc=1000)
-> plot(abc5)
->
-> dprior6 <- function (params, log, ...) {
-+ ll <- sum(
-+ dnorm(
-+ x=params[c("alpha.1","alpha.2","alpha.3","alpha.4")],
-+ mean=c(0.8,-0.5,0.3,0.9),
-+ sd=5,
-+ log=TRUE
-+ )
-+ )
-+ if (log) ll else exp(ll)
-+ }
->
-> abc6 <- abc(pomp(po,dprior=dprior6),
-+ Nabc=2000,
-+ pars=c("alpha.1","alpha.2"),
-+ probes=probes.good,
-+ scale=scale.dat,
-+ epsilon=1,
-+ rw.sd= c(alpha.1=0.01,alpha.2=0.01)
-+ )
-> plot(abc6)
->
-> try(abc7 <- c(abc2,abc3))
-Error in validObject(.Object) :
- invalid class "abcList" object: error in 'c': to be combined, 'abc' objects must have chains of equal length
-> plot(abc7 <- c(abc2,abc4))
-> plot(abc7,scatter=TRUE)
-> plot(conv.rec(c(abc2,abc4)))
-> plot(conv.rec(c(abc7,abc6)))
->
-> dev.off()
-null device
- 1
->
->
-> proc.time()
- user system elapsed
- 11.122 0.067 11.192
Copied: pkg/pomp/tests/ou2-abc.R (from rev 981, pkg/pomp/tests/abc.R)
===================================================================
--- pkg/pomp/tests/ou2-abc.R (rev 0)
+++ pkg/pomp/tests/ou2-abc.R 2014-06-25 17:39:11 UTC (rev 982)
@@ -0,0 +1,103 @@
+### OU2 test of abc for pomp
+
+library(pomp)
+pompExample(ou2)
+
+pdf(file='abc.pdf')
+
+set.seed(2079015564L)
+
+probes.good <- list(
+ y1.mean=probe.mean(var="y1"),
+ y2.mean=probe.mean(var="y2"),
+ probe.acf(var="y1",lags=c(0,5)),
+ probe.acf(var="y2",lags=c(0,5)),
+ probe.ccf(vars=c("y1","y2"),lags=0)
+ )
+psim <- probe(ou2,probes=probes.good,nsim=200)
+plot(psim)
+## why do simulations sometimes seem extreme with respect to these probes?
+
+scale.dat <- apply(psim$simvals,2,sd)
+
+po <- ou2
+
+abc1 <- abc(po,
+ Nabc=10000,
+ start=coef(ou2),
+ pars=c("alpha.1","alpha.2"),
+ probes=probes.good,
+ scale=scale.dat,
+ epsilon=1.7,
+ rw.sd= c(alpha.1=0.01,alpha.2=0.01)
+ )
+
+plot(abc1,scatter=TRUE)
+plot(abc1)
+
+## check how sticky the chain is:
+runs <- rle(as.vector(conv.rec(abc1)[, "alpha.1"]))
+hist(runs$lengths)
+mean(runs$length)
+
+abc2 <- abc(po,
+ Nabc=2000,
+ pars=c("alpha.1","alpha.2"),
+ probes=probes.good,
+ scale=scale.dat,
+ epsilon=1,
+ rw.sd= c(alpha.1=0.01,alpha.2=0.01)
+ )
+plot(abc2)
+
+abc3 <- abc(po,
+ Nabc=2000,
+ probes=probes.good,
+ scale=scale.dat,
+ epsilon=2,
+ rw.sd= c(alpha.1=0.01,alpha.2=0.01)
+ )
+abc3 <- continue(abc3,Nabc=3000)
+plot(abc3)
+
+abc4 <- abc(probe(po,probes=probes.good,nsim=200),
+ Nabc=2000,
+ scale=scale.dat,
+ epsilon=2,
+ rw.sd= c(alpha.1=0.01,alpha.2=0.01)
+ )
+plot(abc4)
+
+abc5 <- abc(abc4,Nabc=1000)
+plot(abc5)
+
+dprior6 <- function (params, log, ...) {
+ ll <- sum(
+ dnorm(
+ x=params[c("alpha.1","alpha.2","alpha.3","alpha.4")],
+ mean=c(0.8,-0.5,0.3,0.9),
+ sd=5,
+ log=TRUE
+ )
+ )
+ if (log) ll else exp(ll)
+}
+
+abc6 <- abc(pomp(po,dprior=dprior6),
+ Nabc=2000,
+ pars=c("alpha.1","alpha.2"),
+ probes=probes.good,
+ scale=scale.dat,
+ epsilon=1,
+ rw.sd= c(alpha.1=0.01,alpha.2=0.01)
+ )
+plot(abc6)
+
+try(abc7 <- c(abc2,abc3))
+plot(abc7 <- c(abc2,abc4))
+plot(abc7,scatter=TRUE)
+plot(conv.rec(c(abc2,abc4)))
+plot(conv.rec(c(abc7,abc6),thin=10,start=5000))
+
+dev.off()
+
Copied: pkg/pomp/tests/ou2-abc.Rout.save (from rev 981, pkg/pomp/tests/abc.Rout.save)
===================================================================
--- pkg/pomp/tests/ou2-abc.Rout.save (rev 0)
+++ pkg/pomp/tests/ou2-abc.Rout.save 2014-06-25 17:39:11 UTC (rev 982)
@@ -0,0 +1,137 @@
+
+R version 3.1.0 (2014-04-10) -- "Spring Dance"
+Copyright (C) 2014 The R Foundation for Statistical Computing
+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.
+
+> ### OU2 test of abc for pomp
+>
+> library(pomp)
+Loading required package: mvtnorm
+Loading required package: subplex
+Loading required package: nloptr
+Loading required package: deSolve
+Loading required package: coda
+Loading required package: lattice
+> pompExample(ou2)
+newly created pomp object(s):
+ ou2
+>
+> pdf(file='abc.pdf')
+>
+> set.seed(2079015564L)
+>
+> probes.good <- list(
++ y1.mean=probe.mean(var="y1"),
++ y2.mean=probe.mean(var="y2"),
++ probe.acf(var="y1",lags=c(0,5)),
++ probe.acf(var="y2",lags=c(0,5)),
++ probe.ccf(vars=c("y1","y2"),lags=0)
++ )
+> psim <- probe(ou2,probes=probes.good,nsim=200)
+> plot(psim)
+> ## why do simulations sometimes seem extreme with respect to these probes?
+>
+> scale.dat <- apply(psim$simvals,2,sd)
+>
+> po <- ou2
+>
+> abc1 <- abc(po,
++ Nabc=10000,
++ start=coef(ou2),
++ pars=c("alpha.1","alpha.2"),
++ probes=probes.good,
++ scale=scale.dat,
++ epsilon=1.7,
++ rw.sd= c(alpha.1=0.01,alpha.2=0.01)
++ )
+>
+> plot(abc1,scatter=TRUE)
+> plot(abc1)
+>
+> ## check how sticky the chain is:
+> runs <- rle(as.vector(conv.rec(abc1)[, "alpha.1"]))
+> hist(runs$lengths)
+> mean(runs$length)
+[1] 6.949965
+>
+> abc2 <- abc(po,
++ Nabc=2000,
++ pars=c("alpha.1","alpha.2"),
++ probes=probes.good,
++ scale=scale.dat,
++ epsilon=1,
++ rw.sd= c(alpha.1=0.01,alpha.2=0.01)
++ )
+> plot(abc2)
+>
+> abc3 <- abc(po,
++ Nabc=2000,
++ probes=probes.good,
++ scale=scale.dat,
++ epsilon=2,
++ rw.sd= c(alpha.1=0.01,alpha.2=0.01)
++ )
+> abc3 <- continue(abc3,Nabc=3000)
+> plot(abc3)
+>
+> abc4 <- abc(probe(po,probes=probes.good,nsim=200),
++ Nabc=2000,
++ scale=scale.dat,
++ epsilon=2,
++ rw.sd= c(alpha.1=0.01,alpha.2=0.01)
++ )
+> plot(abc4)
+>
+> abc5 <- abc(abc4,Nabc=1000)
+> plot(abc5)
+>
+> dprior6 <- function (params, log, ...) {
++ ll <- sum(
++ dnorm(
++ x=params[c("alpha.1","alpha.2","alpha.3","alpha.4")],
++ mean=c(0.8,-0.5,0.3,0.9),
++ sd=5,
++ log=TRUE
++ )
++ )
++ if (log) ll else exp(ll)
++ }
+>
+> abc6 <- abc(pomp(po,dprior=dprior6),
++ Nabc=2000,
++ pars=c("alpha.1","alpha.2"),
++ probes=probes.good,
++ scale=scale.dat,
++ epsilon=1,
++ rw.sd= c(alpha.1=0.01,alpha.2=0.01)
++ )
+> plot(abc6)
+>
+> try(abc7 <- c(abc2,abc3))
+Error in validObject(.Object) :
+ invalid class "abcList" object: error in 'c': to be combined, 'abc' objects must have chains of equal length
+> plot(abc7 <- c(abc2,abc4))
+> plot(abc7,scatter=TRUE)
+> plot(conv.rec(c(abc2,abc4)))
+> plot(conv.rec(c(abc7,abc6),thin=10,start=5000))
+>
+> dev.off()
+null device
+ 1
+>
+>
+> proc.time()
+ user system elapsed
+ 9.980 0.064 10.288
Modified: pkg/pomp/tests/ou2-pmcmc.R
===================================================================
--- pkg/pomp/tests/ou2-pmcmc.R 2014-06-25 17:38:58 UTC (rev 981)
+++ pkg/pomp/tests/ou2-pmcmc.R 2014-06-25 17:39:11 UTC (rev 982)
@@ -76,7 +76,7 @@
plot(ff <- c(ff,f5))
plot(conv.rec(c(f2,ff),c("alpha.2","alpha.3")))
plot(conv.rec(ff[2],c("alpha.2")))
-plot(conv.rec(ff[2:3],c("alpha.3")))
+plot(conv.rec(ff[2:3],c("alpha.3"),thin=3,start=2))
plot(conv.rec(ff[[3]],c("alpha.3")))
dev.off()
Modified: pkg/pomp/tests/ou2-pmcmc.Rout.save
===================================================================
--- pkg/pomp/tests/ou2-pmcmc.Rout.save 2014-06-25 17:38:58 UTC (rev 981)
+++ pkg/pomp/tests/ou2-pmcmc.Rout.save 2014-06-25 17:39:11 UTC (rev 982)
@@ -107,7 +107,7 @@
> plot(ff <- c(ff,f5))
> plot(conv.rec(c(f2,ff),c("alpha.2","alpha.3")))
> plot(conv.rec(ff[2],c("alpha.2")))
-> plot(conv.rec(ff[2:3],c("alpha.3")))
+> plot(conv.rec(ff[2:3],c("alpha.3"),thin=3,start=2))
> plot(conv.rec(ff[[3]],c("alpha.3")))
>
> dev.off()
@@ -117,4 +117,4 @@
>
> proc.time()
user system elapsed
- 23.369 0.040 23.847
+ 23.325 0.088 23.646
More information about the pomp-commits
mailing list