[Pomp-commits] r355 - in pkg: . R data inst man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Oct 1 16:05:22 CEST 2010
Author: kingaa
Date: 2010-10-01 16:05:21 +0200 (Fri, 01 Oct 2010)
New Revision: 355
Added:
pkg/R/compare-mif.R
pkg/R/compare-pmcmc.R
pkg/R/init-state-pomp.R
pkg/tests/ricker-spect.R
Removed:
pkg/R/compare.mif.R
pkg/R/compare.pmcmc.R
pkg/R/init.state-pomp.R
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/mif.R
pkg/R/nlf-guts.R
pkg/R/pmcmc.R
pkg/R/pomp-methods.R
pkg/R/probe-match.R
pkg/R/probe.R
pkg/R/spect-match.R
pkg/R/spect.R
pkg/data/dacca.rda
pkg/data/euler.sir.rda
pkg/data/gillespie.sir.rda
pkg/data/ou2.rda
pkg/data/rw2.rda
pkg/data/verhulst.rda
pkg/inst/NEWS
pkg/man/basic-probes.Rd
pkg/man/pomp-methods.Rd
pkg/man/probe.Rd
pkg/man/spect.Rd
pkg/tests/ou2-probe.R
pkg/tests/ou2-probe.Rout.save
pkg/tests/ricker-probe.R
pkg/tests/ricker-probe.Rout.save
pkg/tests/ricker.R
pkg/tests/ricker.Rout.save
Log:
- fix the behavior of 'coef<-' to make it match simple vector replacement by name. This is not a backward-compatible change but should be stable since it so much simpler, more flexible, and more intuitive.
- drop unneeded 'as.vector' in 'nlf'
- add 'params' argument to 'probe' and 'spect'. This allows the user to override the default 'params=coef(object)'.
- add tests for 'spect' in 'tests/ricker-spect.R'
- drop the slow 'probe.cor' and 'probe.cov' in favor of 'probe.acf' and 'probe.ccf' which are written in C.
- change the names of some source files
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-09-30 19:51:17 UTC (rev 354)
+++ pkg/DESCRIPTION 2010-10-01 14:05:21 UTC (rev 355)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.33-3
-Date: 2010-09-29
+Version: 0.34-1
+Date: 2010-09-30
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>
@@ -15,11 +15,10 @@
LazyData: false
Collate: eulermultinom.R euler.R plugins.R
slice-design.R profile-design.R sobol.R bsplines.R
- pomp-fun.R pomp.R pomp-methods.R rmeasure-pomp.R rprocess-pomp.R
- dmeasure-pomp.R dprocess-pomp.R simulate-pomp.R skeleton-pomp.R
- trajectory-pomp.R plot-pomp.R init.state-pomp.R
+ pomp-fun.R pomp.R pomp-methods.R rmeasure-pomp.R rprocess-pomp.R init-state-pomp.R
+ dmeasure-pomp.R dprocess-pomp.R skeleton-pomp.R simulate-pomp.R trajectory-pomp.R plot-pomp.R
pfilter.R traj-match.R bsmc.R
- mif-class.R particles-mif.R pfilter-mif.R mif.R mif-methods.R compare.mif.R
- pmcmc.R pmcmc-methods.R compare.pmcmc.R
+ mif-class.R particles-mif.R pfilter-mif.R mif.R mif-methods.R compare-mif.R
+ pmcmc.R pmcmc-methods.R compare-pmcmc.R
nlf-funcs.R nlf-guts.R nlf-lql.R nlf-objfun.R nlf.R
probe.R probe-match.R basic-probes.R spect.R spect-match.R
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2010-09-30 19:51:17 UTC (rev 354)
+++ pkg/NAMESPACE 2010-10-01 14:05:21 UTC (rev 355)
@@ -76,8 +76,6 @@
probe.acf,
probe.ccf,
probe.nlar,
- probe.cov,
- probe.cor,
probe.marginal,
probe.match,
spect.match
Copied: pkg/R/compare-mif.R (from rev 351, pkg/R/compare.mif.R)
===================================================================
--- pkg/R/compare-mif.R (rev 0)
+++ pkg/R/compare-mif.R 2010-10-01 14:05:21 UTC (rev 355)
@@ -0,0 +1,103 @@
+compare.mif <- function (z) {
+ ## assumes that x is a list of mifs with identical structure
+ if (!is.list(z)) z <- list(z)
+ if (!all(sapply(z,function(x)is(x,'mif'))))
+ stop("compare.mif error: ",sQuote("z")," must be a mif object or a list of mif objects",call.=FALSE)
+ mar.multi <- c(0,5.1,0,2.1)
+ oma.multi <- c(6,0,5,0)
+ xx <- z[[1]]
+ ivpnames <- xx at ivps
+ estnames <- c(xx at pars,ivpnames)
+ parnames <- names(coef(xx))
+ unestnames <- parnames[-match(estnames,parnames)]
+
+ ## plot filter means
+ filt.diag <- rbind("eff. sample size"=xx at eff.sample.size,filter.mean(xx))
+ filtnames <- rownames(filt.diag)
+ plotnames <- if(length(unestnames)>0) filtnames[-match(unestnames,filtnames)] else filtnames
+ lognames <- filtnames[1] # eff. sample size
+ nplots <- length(plotnames)
+ n.per.page <- min(nplots,10)
+ if(n.per.page<=4) nc <- 1 else nc <- 2
+ nr <- ceiling(n.per.page/nc)
+ oldpar <- par(mar=mar.multi,oma=oma.multi,mfcol=c(nr,nc),ask=dev.interactive(orNone=TRUE))
+ on.exit(par(oldpar))
+ low <- 1
+ hi <- 0
+ time <- time(xx)
+ while (hi<nplots) {
+ hi <- min(low+n.per.page-1,nplots)
+ for (i in seq(from=low,to=hi,by=1)) {
+ n <- i-low+1
+ logplot <- if (plotnames[i]%in%lognames) "y" else ""
+ dat <- sapply(
+ z,
+ function(po, label) {
+ if (label=="eff. sample size")
+ po at eff.sample.size
+ else
+ filter.mean(po,label)
+ },
+ label=plotnames[i]
+ )
+ matplot(
+ y=dat,
+ x=time,
+ axes = FALSE,
+ xlab = "",
+ log=logplot,
+ ylab = "",
+ type = "l"
+ )
+ box()
+ y.side <- 2
+ axis(y.side, xpd = NA)
+ mtext(plotnames[i], y.side, line = 3)
+ do.xax <- (n%%nr==0||n==n.per.page)
+ if (do.xax) axis(1,xpd=NA)
+ if (do.xax) mtext("time",side=1,line=3)
+ }
+ low <- hi+1
+ mtext("Filter diagnostics (last iteration)",3,line=2,outer=TRUE)
+ }
+
+ ## plot mif convergence diagnostics
+ other.diagnostics <- c("loglik", "nfail")
+ plotnames <- c(other.diagnostics,estnames)
+ nplots <- length(plotnames)
+ n.per.page <- min(nplots,10)
+ nc <- if (n.per.page<=4) 1 else 2
+ nr <- ceiling(n.per.page/nc)
+ par(mar=mar.multi,oma=oma.multi,mfcol=c(nr,nc))
+ ## on.exit(par(oldpar))
+ low <- 1
+ hi <- 0
+ iteration <- seq(0,xx at Nmif)
+ while (hi<nplots) {
+ hi <- min(low+n.per.page-1,nplots)
+ for (i in seq(from=low,to=hi,by=1)) {
+ n <- i-low+1
+ dat <- sapply(z,function(po,label) conv.rec(po,label),label=plotnames[i])
+ matplot(
+ y=dat,
+ x=iteration,
+ axes = FALSE,
+ xlab = "",
+ ylab = "",
+ type = "l"
+ )
+ box()
+ y.side <- 2
+ axis(y.side,xpd=NA)
+ mtext(plotnames[i],y.side,line=3)
+ do.xax <- (n%%nr==0||n==n.per.page)
+ if (do.xax) axis(1,xpd=NA)
+ if (do.xax) mtext("MIF iteration",side=1,line=3)
+ }
+ low <- hi+1
+ mtext("MIF convergence diagnostics",3,line=2,outer=TRUE)
+ }
+ invisible(NULL)
+}
+
+
Copied: pkg/R/compare-pmcmc.R (from rev 351, pkg/R/compare.pmcmc.R)
===================================================================
--- pkg/R/compare-pmcmc.R (rev 0)
+++ pkg/R/compare-pmcmc.R 2010-10-01 14:05:21 UTC (rev 355)
@@ -0,0 +1,103 @@
+compare.pmcmc <- function (z) {
+ ## assumes that x is a list of pmcmcs with identical structure
+ if (!is.list(z)) z <- list(z)
+ if (!all(sapply(z,function(x)is(x,'pmcmc'))))
+ stop("compare.pmcmc error: ",sQuote("z")," must be a pmcmc object or a list of pmcmc objects",call.=FALSE)
+ mar.multi <- c(0,5.1,0,2.1)
+ oma.multi <- c(6,0,5,0)
+ xx <- z[[1]]
+ estnames <- xx at pars
+ parnames <- names(coef(xx))
+ unestnames <- parnames[-match(estnames,parnames)]
+
+ ## plot filter means
+ filt.diag <- rbind("eff. sample size"=xx at eff.sample.size,filter.mean(xx))
+ filtnames <- rownames(filt.diag)
+# plotnames <- if(length(unestnames)>0) filtnames[-match(unestnames,filtnames)] else filtnames
+ plotnames <- filtnames
+ lognames <- filtnames[1] # eff. sample size
+ nplots <- length(plotnames)
+ n.per.page <- min(nplots,10)
+ if(n.per.page<=4) nc <- 1 else nc <- 2
+ nr <- ceiling(n.per.page/nc)
+ oldpar <- par(mar=mar.multi,oma=oma.multi,mfcol=c(nr,nc),ask=dev.interactive(orNone=TRUE))
+ on.exit(par(oldpar))
+ low <- 1
+ hi <- 0
+ time <- time(xx)
+ while (hi<nplots) {
+ hi <- min(low+n.per.page-1,nplots)
+ for (i in seq(from=low,to=hi,by=1)) {
+ n <- i-low+1
+ logplot <- if (plotnames[i]%in%lognames) "y" else ""
+ dat <- sapply(
+ z,
+ function(po, label) {
+ if (label=="eff. sample size")
+ po at eff.sample.size
+ else
+ filter.mean(po,label)
+ },
+ label=plotnames[i]
+ )
+ matplot(
+ y=dat,
+ x=time,
+ axes = FALSE,
+ xlab = "",
+ log=logplot,
+ ylab = "",
+ type = "l"
+ )
+ box()
+ y.side <- 2
+ axis(y.side, xpd = NA)
+ mtext(plotnames[i], y.side, line = 3)
+ do.xax <- (n%%nr==0||n==n.per.page)
+ if (do.xax) axis(1,xpd=NA)
+ if (do.xax) mtext("time",side=1,line=3)
+ }
+ low <- hi+1
+ mtext("Filter diagnostics (last iteration)",3,line=2,outer=TRUE)
+ }
+
+ ## plot pmcmc convergence diagnostics
+ other.diagnostics <- c("loglik", "log.prior","nfail")
+ plotnames <- c(other.diagnostics,estnames)
+ nplots <- length(plotnames)
+ n.per.page <- min(nplots,10)
+ nc <- if (n.per.page<=4) 1 else 2
+ nr <- ceiling(n.per.page/nc)
+ par(mar=mar.multi,oma=oma.multi,mfcol=c(nr,nc))
+ ## on.exit(par(oldpar))
+ low <- 1
+ hi <- 0
+ iteration <- seq(0,xx at Nmcmc)
+ while (hi<nplots) {
+ hi <- min(low+n.per.page-1,nplots)
+ for (i in seq(from=low,to=hi,by=1)) {
+ n <- i-low+1
+ dat <- sapply(z,function(po,label) conv.rec(po,label),label=plotnames[i])
+ matplot(
+ y=dat,
+ x=iteration,
+ axes = FALSE,
+ xlab = "",
+ ylab = "",
+ type = "l"
+ )
+ box()
+ y.side <- 2
+ axis(y.side,xpd=NA)
+ mtext(plotnames[i],y.side,line=3)
+ do.xax <- (n%%nr==0||n==n.per.page)
+ if (do.xax) axis(1,xpd=NA)
+ if (do.xax) mtext("PMCMC iteration",side=1,line=3)
+ }
+ low <- hi+1
+ mtext("PMCMC convergence diagnostics",3,line=2,outer=TRUE)
+ }
+ invisible(NULL)
+}
+
+
Deleted: pkg/R/compare.mif.R
===================================================================
--- pkg/R/compare.mif.R 2010-09-30 19:51:17 UTC (rev 354)
+++ pkg/R/compare.mif.R 2010-10-01 14:05:21 UTC (rev 355)
@@ -1,103 +0,0 @@
-compare.mif <- function (z) {
- ## assumes that x is a list of mifs with identical structure
- if (!is.list(z)) z <- list(z)
- if (!all(sapply(z,function(x)is(x,'mif'))))
- stop("compare.mif error: ",sQuote("z")," must be a mif object or a list of mif objects",call.=FALSE)
- mar.multi <- c(0,5.1,0,2.1)
- oma.multi <- c(6,0,5,0)
- xx <- z[[1]]
- ivpnames <- xx at ivps
- estnames <- c(xx at pars,ivpnames)
- parnames <- names(coef(xx))
- unestnames <- parnames[-match(estnames,parnames)]
-
- ## plot filter means
- filt.diag <- rbind("eff. sample size"=xx at eff.sample.size,filter.mean(xx))
- filtnames <- rownames(filt.diag)
- plotnames <- if(length(unestnames)>0) filtnames[-match(unestnames,filtnames)] else filtnames
- lognames <- filtnames[1] # eff. sample size
- nplots <- length(plotnames)
- n.per.page <- min(nplots,10)
- if(n.per.page<=4) nc <- 1 else nc <- 2
- nr <- ceiling(n.per.page/nc)
- oldpar <- par(mar=mar.multi,oma=oma.multi,mfcol=c(nr,nc),ask=dev.interactive(orNone=TRUE))
- on.exit(par(oldpar))
- low <- 1
- hi <- 0
- time <- time(xx)
- while (hi<nplots) {
- hi <- min(low+n.per.page-1,nplots)
- for (i in seq(from=low,to=hi,by=1)) {
- n <- i-low+1
- logplot <- if (plotnames[i]%in%lognames) "y" else ""
- dat <- sapply(
- z,
- function(po, label) {
- if (label=="eff. sample size")
- po at eff.sample.size
- else
- filter.mean(po,label)
- },
- label=plotnames[i]
- )
- matplot(
- y=dat,
- x=time,
- axes = FALSE,
- xlab = "",
- log=logplot,
- ylab = "",
- type = "l"
- )
- box()
- y.side <- 2
- axis(y.side, xpd = NA)
- mtext(plotnames[i], y.side, line = 3)
- do.xax <- (n%%nr==0||n==n.per.page)
- if (do.xax) axis(1,xpd=NA)
- if (do.xax) mtext("time",side=1,line=3)
- }
- low <- hi+1
- mtext("Filter diagnostics (last iteration)",3,line=2,outer=TRUE)
- }
-
- ## plot mif convergence diagnostics
- other.diagnostics <- c("loglik", "nfail")
- plotnames <- c(other.diagnostics,estnames)
- nplots <- length(plotnames)
- n.per.page <- min(nplots,10)
- nc <- if (n.per.page<=4) 1 else 2
- nr <- ceiling(n.per.page/nc)
- par(mar=mar.multi,oma=oma.multi,mfcol=c(nr,nc))
- ## on.exit(par(oldpar))
- low <- 1
- hi <- 0
- iteration <- seq(0,xx at Nmif)
- while (hi<nplots) {
- hi <- min(low+n.per.page-1,nplots)
- for (i in seq(from=low,to=hi,by=1)) {
- n <- i-low+1
- dat <- sapply(z,function(po,label) conv.rec(po,label),label=plotnames[i])
- matplot(
- y=dat,
- x=iteration,
- axes = FALSE,
- xlab = "",
- ylab = "",
- type = "l"
- )
- box()
- y.side <- 2
- axis(y.side,xpd=NA)
- mtext(plotnames[i],y.side,line=3)
- do.xax <- (n%%nr==0||n==n.per.page)
- if (do.xax) axis(1,xpd=NA)
- if (do.xax) mtext("MIF iteration",side=1,line=3)
- }
- low <- hi+1
- mtext("MIF convergence diagnostics",3,line=2,outer=TRUE)
- }
- invisible(NULL)
-}
-
-
Deleted: pkg/R/compare.pmcmc.R
===================================================================
--- pkg/R/compare.pmcmc.R 2010-09-30 19:51:17 UTC (rev 354)
+++ pkg/R/compare.pmcmc.R 2010-10-01 14:05:21 UTC (rev 355)
@@ -1,103 +0,0 @@
-compare.pmcmc <- function (z) {
- ## assumes that x is a list of pmcmcs with identical structure
- if (!is.list(z)) z <- list(z)
- if (!all(sapply(z,function(x)is(x,'pmcmc'))))
- stop("compare.pmcmc error: ",sQuote("z")," must be a pmcmc object or a list of pmcmc objects",call.=FALSE)
- mar.multi <- c(0,5.1,0,2.1)
- oma.multi <- c(6,0,5,0)
- xx <- z[[1]]
- estnames <- xx at pars
- parnames <- names(coef(xx))
- unestnames <- parnames[-match(estnames,parnames)]
-
- ## plot filter means
- filt.diag <- rbind("eff. sample size"=xx at eff.sample.size,filter.mean(xx))
- filtnames <- rownames(filt.diag)
-# plotnames <- if(length(unestnames)>0) filtnames[-match(unestnames,filtnames)] else filtnames
- plotnames <- filtnames
- lognames <- filtnames[1] # eff. sample size
- nplots <- length(plotnames)
- n.per.page <- min(nplots,10)
- if(n.per.page<=4) nc <- 1 else nc <- 2
- nr <- ceiling(n.per.page/nc)
- oldpar <- par(mar=mar.multi,oma=oma.multi,mfcol=c(nr,nc),ask=dev.interactive(orNone=TRUE))
- on.exit(par(oldpar))
- low <- 1
- hi <- 0
- time <- time(xx)
- while (hi<nplots) {
- hi <- min(low+n.per.page-1,nplots)
- for (i in seq(from=low,to=hi,by=1)) {
- n <- i-low+1
- logplot <- if (plotnames[i]%in%lognames) "y" else ""
- dat <- sapply(
- z,
- function(po, label) {
- if (label=="eff. sample size")
- po at eff.sample.size
- else
- filter.mean(po,label)
- },
- label=plotnames[i]
- )
- matplot(
- y=dat,
- x=time,
- axes = FALSE,
- xlab = "",
- log=logplot,
- ylab = "",
- type = "l"
- )
- box()
- y.side <- 2
- axis(y.side, xpd = NA)
- mtext(plotnames[i], y.side, line = 3)
- do.xax <- (n%%nr==0||n==n.per.page)
- if (do.xax) axis(1,xpd=NA)
- if (do.xax) mtext("time",side=1,line=3)
- }
- low <- hi+1
- mtext("Filter diagnostics (last iteration)",3,line=2,outer=TRUE)
- }
-
- ## plot pmcmc convergence diagnostics
- other.diagnostics <- c("loglik", "log.prior","nfail")
- plotnames <- c(other.diagnostics,estnames)
- nplots <- length(plotnames)
- n.per.page <- min(nplots,10)
- nc <- if (n.per.page<=4) 1 else 2
- nr <- ceiling(n.per.page/nc)
- par(mar=mar.multi,oma=oma.multi,mfcol=c(nr,nc))
- ## on.exit(par(oldpar))
- low <- 1
- hi <- 0
- iteration <- seq(0,xx at Nmcmc)
- while (hi<nplots) {
- hi <- min(low+n.per.page-1,nplots)
- for (i in seq(from=low,to=hi,by=1)) {
- n <- i-low+1
- dat <- sapply(z,function(po,label) conv.rec(po,label),label=plotnames[i])
- matplot(
- y=dat,
- x=iteration,
- axes = FALSE,
- xlab = "",
- ylab = "",
- type = "l"
- )
- box()
- y.side <- 2
- axis(y.side,xpd=NA)
- mtext(plotnames[i],y.side,line=3)
- do.xax <- (n%%nr==0||n==n.per.page)
- if (do.xax) axis(1,xpd=NA)
- if (do.xax) mtext("PMCMC iteration",side=1,line=3)
- }
- low <- hi+1
- mtext("PMCMC convergence diagnostics",3,line=2,outer=TRUE)
- }
- invisible(NULL)
-}
-
-
Copied: pkg/R/init-state-pomp.R (from rev 351, pkg/R/init.state-pomp.R)
===================================================================
--- pkg/R/init-state-pomp.R (rev 0)
+++ pkg/R/init-state-pomp.R 2010-10-01 14:05:21 UTC (rev 355)
@@ -0,0 +1,20 @@
+init.state <- function (object, params, t0, ...)
+ stop("function ",sQuote("init.state")," is undefined for objects of class ",sQuote(class(object)))
+setGeneric('init.state')
+
+## initialize the process model
+setMethod(
+ 'init.state',
+ 'pomp',
+ function (object, params, t0, ...) { # the package algorithms will only use these arguments
+ if (missing(t0)) t0 <- object at t0
+ if (missing(params)) params <- coef(object)
+ x <- try(
+ .Call(do_init_state,object,as.matrix(params),t0),
+ silent=FALSE
+ )
+ if (inherits(x,'try-error'))
+ stop("init.state error: error in user ",sQuote("initializer"),call.=FALSE)
+ x
+ }
+ )
Deleted: pkg/R/init.state-pomp.R
===================================================================
--- pkg/R/init.state-pomp.R 2010-09-30 19:51:17 UTC (rev 354)
+++ pkg/R/init.state-pomp.R 2010-10-01 14:05:21 UTC (rev 355)
@@ -1,20 +0,0 @@
-init.state <- function (object, params, t0, ...)
- stop("function ",sQuote("init.state")," is undefined for objects of class ",sQuote(class(object)))
-setGeneric('init.state')
-
-## initialize the process model
-setMethod(
- 'init.state',
- 'pomp',
- function (object, params, t0, ...) { # the package algorithms will only use these arguments
- if (missing(t0)) t0 <- object at t0
- if (missing(params)) params <- coef(object)
- x <- try(
- .Call(do_init_state,object,as.matrix(params),t0),
- silent=FALSE
- )
- if (inherits(x,'try-error'))
- stop("init.state error: error in user ",sQuote("initializer"),call.=FALSE)
- x
- }
- )
Modified: pkg/R/mif.R
===================================================================
--- pkg/R/mif.R 2010-09-30 19:51:17 UTC (rev 354)
+++ pkg/R/mif.R 2010-10-01 14:05:21 UTC (rev 355)
@@ -270,7 +270,7 @@
}
- coef(obj,names(theta)) <- unname(theta)
+ coef(obj) <- theta
if (Nmif>0) {
obj at Nmif <- as.integer(Nmif)
Modified: pkg/R/nlf-guts.R
===================================================================
--- pkg/R/nlf-guts.R 2010-09-30 19:51:17 UTC (rev 354)
+++ pkg/R/nlf-guts.R 2010-10-01 14:05:21 UTC (rev 355)
@@ -162,7 +162,7 @@
for (jvar in seq_len(nvar)) {
model.lm <- lm(model.pred[,jvar]~rbfbasis.model-1)
model.residuals[,jvar] <- residuals(model.lm)
- ck <- as.vector(coef(model.lm))
+ ck <- coef(model.lm)
if (verbose) {
print(ck)
print(summary(model.lm))
Modified: pkg/R/pmcmc.R
===================================================================
--- pkg/R/pmcmc.R 2010-09-30 19:51:17 UTC (rev 354)
+++ pkg/R/pmcmc.R 2010-10-01 14:05:21 UTC (rev 355)
@@ -246,7 +246,7 @@
}
- coef(obj) <- unname(theta)
+ coef(obj) <- theta
if (Nmcmc>0) {
obj at Nmcmc <- as.integer(Nmcmc)
Modified: pkg/R/pomp-methods.R
===================================================================
--- pkg/R/pomp-methods.R 2010-09-30 19:51:17 UTC (rev 354)
+++ pkg/R/pomp-methods.R 2010-10-01 14:05:21 UTC (rev 355)
@@ -169,37 +169,36 @@
"coef<-",
"pomp",
function (object, pars, ..., value) {
- if (length(object at params)==0) {
- if (missing(pars)) {
- pars <- names(value)
- if (is.null(pars))
- stop("in ",sQuote("coef<-"),": ",sQuote("value")," must be a named vector")
- } else {
- if (length(pars)!=length(value))
- stop("in ",sQuote("coef<-"),": ",sQuote("pars")," and ",sQuote("value")," must be of the same length")
- }
- object at params <- as.numeric(value)
+ if (missing(pars)) { ## replace the whole params slot with 'value'
+ pars <- names(value)
+ if (is.null(pars))
+ stop("in ",sQuote("coef<-"),": ",sQuote("value")," must be a named vector")
+ object at params <- numeric(length(pars))
names(object at params) <- pars
- } else {
- if (missing(pars)) {
- pars <- names(object at params)
- if (is.null(pars))
- stop("bad ",sQuote("pomp")," object: slot ",sQuote("params")," should be a named vector")
- } else {
- excl <- !(pars%in%names(object at params))
- if (any(excl)) {
- stop(
- "in ",sQuote("coef<-"),": name(s) ",
- paste(sapply(pars[excl],sQuote),collapse=","),
- " correspond to no parameter(s)"
- )
+ object at params[] <- as.numeric(value)
+ } else { ## replace or append only the parameters named in 'pars'
+ if (!is.null(names(value))) ## we ignore the names of 'value'
+ warning("in ",sQuote("coef<-"),": names of ",sQuote("value")," are being discarded",call.=FALSE)
+ if (length(object at params)==0) { ## no pre-existing 'params' slot
+ object at params <- numeric(length(pars))
+ names(object at params) <- pars
+ object at params[] <- as.numeric(value)
+ } else { ## pre-existing params slot
+ excl <- !(pars%in%names(object at params)) ## new parameters
+ if (any(excl)) { ## append parameters
+ warning(
+ "in ",sQuote("coef<-"),": name(s) ",
+ paste(sQuote(pars[excl]),collapse=","),
+ " are not existing parameter(s);",
+ " they are being concatenated",
+ call.=FALSE
+ )
+ x <- c(object at params,numeric(length(excl)))
+ names(x) <- c(names(object at params),pars[excl])
+ object at params <- x
}
+ object at params[pars] <- as.numeric(value)
}
- if (length(pars)!=length(value))
- stop("in ",sQuote("coef<-"),": ",sQuote("pars")," and ",sQuote("value")," must be of the same length")
- if (!is.null(names(value)))
- warning("in ",sQuote("coef<-"),": names of ",sQuote("value")," are being discarded",call.=FALSE)
- object at params[pars] <- as.numeric(value)
}
object
}
Modified: pkg/R/probe-match.R
===================================================================
--- pkg/R/probe-match.R 2010-09-30 19:51:17 UTC (rev 354)
+++ pkg/R/probe-match.R 2010-10-01 14:05:21 UTC (rev 355)
@@ -162,11 +162,9 @@
msg <- opt$message
}
- coef(object,names(params)) <- unname(params)
-
new(
"probe.matched.pomp",
- probe(object,probes=probes,nsim=nsim,seed=seed),
+ probe(object,probes=probes,params=params,nsim=nsim,seed=seed),
weights=weights,
fail.value=as.numeric(fail.value),
value=val,
Modified: pkg/R/probe.R
===================================================================
--- pkg/R/probe.R 2010-09-30 19:51:17 UTC (rev 354)
+++ pkg/R/probe.R 2010-10-01 14:05:21 UTC (rev 355)
@@ -29,12 +29,16 @@
setMethod(
"probe",
signature(object="pomp"),
- function (object, probes, nsim = 1, seed = NULL, ...) {
+ function (object, probes, params, nsim = 1, seed = NULL, ...) {
+
if (!is.list(probes)) probes <- list(probes)
if (!all(sapply(probes,is.function)))
stop(sQuote("probes")," must be a function or a list of functions")
if (!all(sapply(probes,function(f)length(formals(f))==1)))
stop("each probe must be a function of a single argument")
+
+ if (missing(params)) params <- coef(object)
+
if (is.null(seed)) {
if (exists('.Random.seed',where=.GlobalEnv)) {
seed <- get(".Random.seed",pos=.GlobalEnv)
@@ -46,12 +50,12 @@
## apply probes to model simulations
simval <- .Call(
apply_probe_sim,
- object,
- nsim,
- coef(object),
- seed,
- probes,
- datval
+ object=object,
+ nsim=nsim,
+ params=params,
+ seed=seed,
+ probes=probes,
+ datval=datval
)
nprobes <- length(datval)
@@ -66,6 +70,8 @@
quants[k] <- sum(simval[,k]<datval[k])/nsim
}
+ coef(object) <- params
+
new(
"probed.pomp",
object,
@@ -82,20 +88,27 @@
setMethod(
"probe",
signature(object="probed.pomp"),
- function (object, probes, nsim, seed = NULL, ...) {
+ function (object, probes, params, nsim, seed = NULL, ...) {
+
if (missing(probes)) probes <- object at probes
if (!is.list(probes)) probes <- list(probes)
if (!all(sapply(probes,is.function)))
stop(sQuote("probes")," must be a function or a list of functions")
+
+ if (missing(params)) params <- coef(object)
+
if (is.null(seed)) {
if (exists('.Random.seed',where=.GlobalEnv)) {
seed <- get(".Random.seed",pos=.GlobalEnv)
}
}
+
if (missing(nsim)) nsim <- nrow(object at simvals)
+
probe(
as(object,"pomp"),
probes=probes,
+ params=params,
nsim=nsim,
seed=seed,
...
Modified: pkg/R/spect-match.R
===================================================================
--- pkg/R/spect-match.R 2010-09-30 19:51:17 UTC (rev 354)
+++ pkg/R/spect-match.R 2010-10-01 14:05:21 UTC (rev 355)
@@ -220,12 +220,11 @@
msg <- opt$message
}
- coef(object,names(params)) <- unname(params)
-
new(
"spect.matched.pomp",
spect(
object,
+ params=params,
vars=vars,
kernel.width=kernel.width,
nsim=nsim,
Modified: pkg/R/spect.R
===================================================================
--- pkg/R/spect.R 2010-09-30 19:51:17 UTC (rev 354)
+++ pkg/R/spect.R 2010-10-01 14:05:21 UTC (rev 355)
@@ -91,7 +91,7 @@
list(freq=freq,spec=datspec)
}
-compute.spect.sim <- function (object, vars, params, nsim, seed, transform, detrend, ker) {
+compute.spect.sim <- function (object, params, vars, nsim, seed, transform, detrend, ker) {
sims <- try(
simulate(
object,
@@ -140,10 +140,13 @@
setMethod(
"spect",
signature(object="pomp"),
- function (object, vars, kernel.width, nsim, seed = NULL, transform = identity,
+ function (object, params, vars, kernel.width, nsim, seed = NULL,
+ transform = identity,
detrend = c("none","mean","linear","quadratic"),
...) {
+ if (missing(params)) params <- coef(object)
+
if (missing(vars))
vars <- rownames(object at data)
@@ -172,8 +175,8 @@
datspec <- ds$spec
simspec <- compute.spect.sim(
object,
+ params=params,
vars=vars,
- params=coef(object),
nsim=nsim,
seed=seed,
transform=transform,
@@ -200,6 +203,8 @@
}
pvals[length(vars)+1] <- (nsim+1-sum(totsimdist<totdatdist))/(nsim+1)
+ coef(object) <- params
+
new(
"spect.pomp",
object,
@@ -218,7 +223,8 @@
setMethod(
"spect",
signature(object="spect.pomp"),
- function (object, vars, kernel.width, nsim, seed = NULL, transform, detrend, ...) {
+ function (object, params, vars, kernel.width, nsim, seed = NULL, transform, detrend, ...) {
+ if (missing(params)) params <- coef(object)
if (missing(vars)) probes <- rownames(object at datspec)
if (missing(kernel.width)) kernel.width <- object at kernel.width
if (missing(nsim)) nsim <- nrow(object at simspec)
@@ -231,6 +237,7 @@
if (missing(detrend)) detrend <- object at detrend
spect(
as(object,"pomp"),
+ params=params,
vars=vars,
kernel.width=kernel.width,
nsim=nsim,
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/ou2.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 2010-09-30 19:51:17 UTC (rev 354)
+++ pkg/inst/NEWS 2010-10-01 14:05:21 UTC (rev 355)
@@ -1,3 +1,8 @@
+NEW IN VERSION 0.34-1:
+ o 'coef<-' now behaves like ordinary vector assignment.
+ o 'probe' and 'spect' now take an argument 'params'
+ o 'probe.cov' and 'probe.cor' have been removed in favor of 'probe.acf' and 'probe.ccf'
+
NEW IN VERSION 0.33-1:
o Major improvements in speed in the probe-matching routines have been realized by coding the computationally-intensive portions of these calculations in C.
Modified: pkg/man/basic-probes.Rd
===================================================================
--- pkg/man/basic-probes.Rd 2010-09-30 19:51:17 UTC (rev 354)
+++ pkg/man/basic-probes.Rd 2010-10-01 14:05:21 UTC (rev 355)
@@ -6,8 +6,6 @@
\alias{probe.sd}
\alias{probe.period}
\alias{probe.quantile}
-\alias{probe.cov}
-\alias{probe.cor}
\alias{probe.acf}
\alias{probe.ccf}
\alias{probe.marginal}
@@ -27,10 +25,6 @@
probe.acf(var, lag.max, type = c("covariance", "correlation"),
transform = identity)
probe.ccf(vars, lags, transform = identity)
-probe.cov(vars, lag, method = c("pearson", "kendall", "spearman"),
- transform = identity)
-probe.cor(vars, lag, method = c("pearson", "kendall", "spearman"),
- transform = identity)
probe.period(var, kernel.width, transform = identity)
probe.quantile(var, prob, transform = identity)
}
@@ -57,28 +51,27 @@
\item{lag.max}{
The maximum lag at which the ACF is computed.
}
- \item{lag, lags}{
- In \code{probe.cov} and \code{probe.cor}, the lag between time series.
+ \item{lags}{
In \code{probe.ccf}, a vector of lags between time series.
Positive lags correspond to \code{x} advanced relative to \code{y};
negative lags, to the reverse.
+
+ In \code{probe.nlar}, a vector of lags present in the nonlinear autoregressive model that will be fit to the actual and simulated data.
+ See Details, below, for a precise description.
}
+ \item{powers}{
+ the powers of each term (corresponding to \code{lags}) in the the nonlinear autoregressive model that will be fit to the actual and simulated data.
+ See Details, below, for a precise description.
+ }
\item{type}{
Compute autocorrelation or autocovariance?
}
- \item{method}{
- The type of covariance or correlation to compute: see \code{\link{cov}} and \code{\link{cor}}.
- }
\item{ref}{
empirical reference distribution.
Simulated data will be regressed against the values of \code{ref}, sorted and, optionally, differenced.
The resulting regression coefficients capture information about the shape of the marginal distribution.
A good choice for \code{ref} is the data itself.
}
- \item{lags, powers}{
- the lags and powers, respectively, specifying the nonlinear autoregressive model that will be fit to the actual and simulated data.
- See Details, below, for a precise description.
- }
\item{order}{
order of polynomial regression.
}
Modified: pkg/man/pomp-methods.Rd
===================================================================
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 355
More information about the pomp-commits
mailing list