[Pomp-commits] r361 - in pkg: . R data inst inst/data-R inst/doc man src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Oct 5 18:26:16 CEST 2010
Author: kingaa
Date: 2010-10-05 18:26:16 +0200 (Tue, 05 Oct 2010)
New Revision: 361
Removed:
pkg/R/nlf-lql.R
Modified:
pkg/DESCRIPTION
pkg/R/aaa.R
pkg/R/nlf-objfun.R
pkg/R/nlf.R
pkg/R/pomp-methods.R
pkg/R/simulate-pomp.R
pkg/R/spect.R
pkg/R/traj-match.R
pkg/R/trajectory-pomp.R
pkg/data/euler.sir.rda
pkg/data/gillespie.sir.rda
pkg/data/ou2.rda
pkg/data/ricker.rda
pkg/data/rw2.rda
pkg/data/verhulst.rda
pkg/inst/ChangeLog
pkg/inst/NEWS
pkg/inst/data-R/ou2.R
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/intro_to_pomp.Rnw
pkg/inst/doc/intro_to_pomp.pdf
pkg/inst/doc/ou2-first-mif.rda
pkg/inst/doc/ou2-trajmatch.rda
pkg/man/pomp-class.Rd
pkg/man/simulate-pomp.Rd
pkg/man/trajectory-pomp.Rd
pkg/man/verhulst.Rd
pkg/src/probe.c
pkg/tests/gillespie.R
pkg/tests/gillespie.Rout.save
pkg/tests/logistic.R
pkg/tests/logistic.Rout.save
pkg/tests/ou2-forecast.R
pkg/tests/ou2-forecast.Rout.save
pkg/tests/ou2-nlf.Rout.save
pkg/tests/ou2-probe.Rout.save
pkg/tests/ou2-simulate.Rout.save
pkg/tests/ou2-trajmatch.Rout.save
pkg/tests/ricker-probe.Rout.save
pkg/tests/ricker.R
pkg/tests/ricker.Rout.save
pkg/tests/rw2.Rout.save
pkg/tests/sir.Rout.save
Log:
- The default behaviors of 'simulate' and 'trajectory' with respect to the argument 'times' has changed. There is a new argument, 't0', which specifies the time at which the initial conditions obtain. The default for 't0' is 'timezero(object)'. The default for 'times' is now 'time(object,t0=FALSE)' (it was 'time(object,t0=TRUE)'). This change eliminates an infelicity in which 'simulate' and 'trajectory' returned undesired values at the initial time. To aid users in modifying their codes, a warning is now displayed whenever a potentially affected call is made.
- A banner is now displayed on package attachment. This will be used to report important changes.
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-10-04 19:53:57 UTC (rev 360)
+++ pkg/DESCRIPTION 2010-10-05 16:26:16 UTC (rev 361)
@@ -19,6 +19,6 @@
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
- nlf-funcs.R nlf-guts.R nlf-lql.R nlf-objfun.R nlf.R
+ pmcmc.R pmcmc-methods.R compare-pmcmc.R
+ nlf-funcs.R nlf-guts.R nlf-objfun.R nlf.R
probe.R probe-match.R basic-probes.R spect.R spect-match.R
Modified: pkg/R/aaa.R
===================================================================
--- pkg/R/aaa.R 2010-10-04 19:53:57 UTC (rev 360)
+++ pkg/R/aaa.R 2010-10-05 16:26:16 UTC (rev 361)
@@ -1,5 +1,26 @@
-.onAttach <- function(...) {
- version <- packageDescription(pkg="pomp",fields="Version")
- cat("This is pomp version ",version,"\n",sep="")
- cat("\n",sep="")
+.onAttach <- function (...) {
+ version <- library(help=pomp)$info[[1]]
+ version <- strsplit(version[pmatch("Version",version)]," ")[[1]]
+ version <- version[nchar(version)>0][2]
+ cat("This is pomp version ",version,"\n\n",sep="")
+ cat(
+ paste(
+ "IMPORTANT NOTICE:\n",
+ "The default behaviors of ",sQuote("simulate")," and ",sQuote("trajectory"),
+ " have changed as of release 0.34-1. ",
+ "You can ensure that your code will continue to function as you intend by specifying the values ",
+ "of the ",sQuote("times")," and ",sQuote("t0")," arguments to these functions, thus removing ",
+ "dependence of your code on the defaults. ",
+ "In the meantime, using ",sQuote("simulate")," or ",sQuote("trajectory"),
+ " in such a way as to rely on the default will produce a warning. ",
+ "These warnings will be removed in a future release.\n\n",
+ "See the documentation (",dQuote("pomp?simulate"),", ",dQuote("pomp?trajectory"),
+ ") for a description of the new default behaviors.\n\n",
+ "Subscribe to the pomp-announce list to receive email notification about new releases of pomp ",
+ "including descriptions of feature additions, changes, and fixes.\n\n",
+ sep=""
+ ),
+ fill=TRUE,
+ sep=""
+ )
}
Deleted: pkg/R/nlf-lql.R
===================================================================
--- pkg/R/nlf-lql.R 2010-10-04 19:53:57 UTC (rev 360)
+++ pkg/R/nlf-lql.R 2010-10-05 16:26:16 UTC (rev 361)
@@ -1,62 +0,0 @@
-NLF.LQL <- function (params.fitted, object, params, par.index,
- times, lags, period, tensor, seed = NULL, transform = identity,
- nrbf = 4, verbose = FALSE,
- bootstrap = FALSE, bootsamp = NULL) {
-
-#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-# Computes the vector of log quasi-likelihood values at the observations
-# Note that the log QL itself is returned, not the negative log QL,
-# so a large NEGATIVE value is used to flag bad parameters
-#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-
- FAILED = -99999999999
- params[par.index] <- params.fitted
-
- params <- as.matrix(params)
-
- ## Need to extract number of state variables (nvar) from pomp object
- ## Need to include simulation times in problem specification
- ## Evaluates the NLF objective function given a POMP object.
- ## Version 0.1, 3 Dec. 2007, Bruce E. Kendall & Stephen P. Ellner
- ## Version 0.2, May 2008, Stephen P. Ellner
-
- data.ts <- data.array(object)
-
- y <- try(
- simulate(object,times=times,params=params,seed=seed,obs=TRUE,states=FALSE),
- silent=FALSE
- )
- if (inherits(y,"try-error"))
- stop(sQuote("NLF.LQL")," reports: error in simulation")
- ## Test whether the model time series is valid
- if (!all(is.finite(y))) return(FAILED)
-
- model.ts <- array(
- dim=c(nrow(data.ts),length(times)-1),
- dimnames=list(rownames(data.ts),NULL)
- )
- model.ts[,] <- apply(y[,1,-1,drop=FALSE],c(2,3),transform)
- data.ts[,] <- apply(data.ts,2,transform)
-
- LQL <- try(
- NLF.guts(
- data.mat=data.ts,
- data.times=time(object),
- model.mat=model.ts,
- model.times=times[-1],
- lags=lags,
- period=period,
- tensor=tensor,
- nrbf=nrbf,
- verbose=FALSE,
- bootstrap,
- bootsamp,
- plotfit=FALSE
- ),
- silent=FALSE
- )
- if (inherits(LQL,"try-error"))
- stop(sQuote("NLF.LQL")," reports: error in ",sQuote("NLF.guts"))
- LQL
- }
-
Modified: pkg/R/nlf-objfun.R
===================================================================
--- pkg/R/nlf-objfun.R 2010-10-04 19:53:57 UTC (rev 360)
+++ pkg/R/nlf-objfun.R 2010-10-05 16:26:16 UTC (rev 361)
@@ -1,5 +1,69 @@
+NLF.LQL <- function (params.fitted, object, params, par.index,
+ times, t0, lags, period, tensor, seed = NULL, transform = identity,
+ nrbf = 4, verbose = FALSE,
+ bootstrap = FALSE, bootsamp = NULL) {
+
+###>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+### Computes the vector of log quasi-likelihood values at the observations
+### Note that the log QL itself is returned, not the negative log QL,
+### so a large NEGATIVE value is used to flag bad parameters
+###>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+
+ FAILED = -99999999999
+ params[par.index] <- params.fitted
+
+ params <- as.matrix(params)
+
+ ## Need to extract number of state variables (nvar) from pomp object
+ ## Need to include simulation times in problem specification
+ ## Evaluates the NLF objective function given a POMP object.
+ ## Version 0.1, 3 Dec. 2007, Bruce E. Kendall & Stephen P. Ellner
+ ## Version 0.2, May 2008, Stephen P. Ellner
+
+ data.ts <- data.array(object)
+
+ y <- try(
+ simulate(object,times=times,t0=t0,params=params,seed=seed,obs=TRUE,states=FALSE),
+ silent=FALSE
+ )
+ if (inherits(y,"try-error"))
+ stop(sQuote("NLF.LQL")," reports: error in simulation")
+ ## Test whether the model time series is valid
+ if (!all(is.finite(y))) return(FAILED)
+
+ model.ts <- array(
+ dim=c(nrow(data.ts),length(times)),
+ dimnames=list(rownames(data.ts),NULL)
+ )
+ model.ts[,] <- apply(y[,1,,drop=FALSE],c(2,3),transform)
+ data.ts[,] <- apply(data.ts,2,transform)
+
+ LQL <- try(
+ NLF.guts(
+ data.mat=data.ts,
+ data.times=time(object),
+ model.mat=model.ts,
+ model.times=times,
+ lags=lags,
+ period=period,
+ tensor=tensor,
+ nrbf=nrbf,
+ verbose=FALSE,
+ bootstrap,
+ bootsamp,
+ plotfit=FALSE
+ ),
+ silent=FALSE
+ )
+ if (inherits(LQL,"try-error"))
+ stop(sQuote("NLF.LQL")," reports: error in ",sQuote("NLF.guts"))
+ LQL
+}
+
+
+
nlf.objfun <- function (params.fitted, object, params, par.index,
- times, lags, period, tensor, seed, transform = function(x)x,
+ times, t0, lags, period, tensor, seed, transform = function(x)x,
nrbf = 4, verbose = FALSE, bootstrap = FALSE, bootsamp = NULL)
{
-sum(
@@ -9,6 +73,7 @@
params=params,
par.index=par.index,
times=times,
+ t0=t0,
lags=lags,
period=period,
tensor=tensor,
Modified: pkg/R/nlf.R
===================================================================
--- pkg/R/nlf.R 2010-10-04 19:53:57 UTC (rev 360)
+++ pkg/R/nlf.R 2010-10-05 16:26:16 UTC (rev 361)
@@ -53,21 +53,20 @@
guess <- params[par.index]
dt.tol <- 1e-3
- times <- time(object,t0=TRUE)
- dt <- diff(times[-1])
+ times <- time(object,t0=FALSE)
+ t0 <- timezero(object)
+ dt <- diff(times)
if (diff(range(dt))>dt.tol*mean(dt))
stop(sQuote("nlf")," requires evenly spaced sampling times")
- dt <- times[3]-times[2]
- t0 <- times[1]
+ dt <- times[2]-times[1]
## Vector of times to output the simulation
times <- seq(
from=t0,
- to=t0+(nconverge+nasymp)*dt,
+ length=nconverge+nasymp+1,
by=dt
)
-
if (eval.only) {
val <- nlf.objfun(
params.fitted=guess,
@@ -75,6 +74,7 @@
params=params,
par.index=par.index,
times=times,
+ t0=t0,
lags=lags,
period=period,
tensor=tensor,
@@ -96,6 +96,7 @@
params=params,
par.index=par.index,
times=times,
+ t0=t0,
lags=lags,
period=period,
tensor=tensor,
@@ -117,6 +118,7 @@
params=params,
par.index=par.index,
times=times,
+ t0=t0,
lags=lags,
period=period,
tensor=tensor,
@@ -141,7 +143,8 @@
Jhat <- matrix(0,nfitted,nfitted)
Ihat <- Jhat
f0 <- NLF.LQL(fitted,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor, seed=seed,
+ times=times, t0=t0,
+ lags=lags, period=period, tensor=tensor, seed=seed,
transform=transform, nrbf=4,
verbose=FALSE)
F0 <- mean(f0,na.rm=T)
@@ -161,25 +164,25 @@
guess <- fitted
guess[i] <- fitted[i]-sqrt(2)*h*abs(fitted[i])
Fvals[1] <- mean(NLF.LQL(guess,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor,
+ times=times, t0=t0, lags=lags, period=period, tensor=tensor,
seed=seed, transform=transform,
nrbf=4, verbose=FALSE),na.rm=T)
guess <- fitted
guess[i] <- fitted[i]-h*abs(fitted[i])
Fvals[2] <- mean(NLF.LQL(guess,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor,
+ times=times, t0=t0, lags=lags, period=period, tensor=tensor,
seed=seed, transform=transform, nrbf=4,
verbose=FALSE),na.rm=T)
guess <- fitted
guess[i] <- fitted[i]+h*abs(fitted[i])
Fvals[4] <- mean(NLF.LQL(guess,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor,
+ times=times, t0=t0, lags=lags, period=period, tensor=tensor,
seed=seed, transform=transform, nrbf=4,
verbose=FALSE),na.rm=T)
guess <- fitted
guess[i] <- fitted[i]+sqrt(2)*h*abs(fitted[i])
Fvals[5] <- mean(NLF.LQL(guess,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor,
+ times=times, t0=t0, lags=lags, period=period, tensor=tensor,
seed=seed, transform=transform, nrbf=4,
verbose=FALSE),na.rm=T)
FAILED = - 999999
@@ -197,13 +200,13 @@
guess.up <- fitted
guess.up[i] <- guess.up[i]+eps[i]
f.up <- NLF.LQL(guess.up,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor,
+ times=times, t0=t0, lags=lags, period=period, tensor=tensor,
seed=seed, transform=transform, nrbf=4,
verbose=FALSE)
F.up <- mean(f.up,na.rm=T)
f.up2 <- NLF.LQL(guess.up,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor,
+ times=times, t0=t0, lags=lags, period=period, tensor=tensor,
seed=seed, transform=transform, nrbf=4,
verbose=FALSE)
@@ -213,7 +216,7 @@
guess.down <- fitted
guess.down[i] <- guess.down[i]-eps[i]
f.down <- NLF.LQL(guess.down,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor,
+ times=times, t0=t0, lags=lags, period=period, tensor=tensor,
seed=seed, transform=transform, nrbf=4,
verbose=FALSE)
F.down <- mean(f.down,na.rm=T)
@@ -232,7 +235,7 @@
guess.uu[i] <- guess.uu[i]+eps[i]
guess.uu[j] <- guess.uu[j]+eps[j]
F.uu <- mean(NLF.LQL(guess.uu,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor,
+ times=times, t0=t0, lags=lags, period=period, tensor=tensor,
seed=seed, transform=transform, nrbf=4,
verbose=FALSE),na.rm=T)
@@ -240,7 +243,7 @@
guess.ud[i] <- guess.ud[i]+eps[i]
guess.ud[j] <- guess.ud[j]-eps[j]
F.ud <- mean(NLF.LQL(guess.ud,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor,
+ times=times, t0=t0, lags=lags, period=period, tensor=tensor,
seed=seed, transform=transform, nrbf=4,
verbose=FALSE),na.rm=T)
@@ -248,7 +251,7 @@
guess.du[i] <- guess.du[i]-eps[i]
guess.du[j] <- guess.du[j]+eps[j]
F.du <- mean(NLF.LQL(guess.du,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor,
+ times=times, t0=t0, lags=lags, period=period, tensor=tensor,
seed=seed, transform=transform, nrbf=4,
verbose=FALSE),na.rm=T)
@@ -256,7 +259,7 @@
guess.dd[i] <- guess.dd[i]-eps[i]
guess.dd[j] <- guess.dd[j]-eps[j]
F.dd <- mean(NLF.LQL(guess.dd,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor,
+ times=times, t0=t0, lags=lags, period=period, tensor=tensor,
seed=seed, transform=transform, nrbf=4,
verbose=FALSE),na.rm=T)
Modified: pkg/R/pomp-methods.R
===================================================================
--- pkg/R/pomp-methods.R 2010-10-04 19:53:57 UTC (rev 360)
+++ pkg/R/pomp-methods.R 2010-10-05 16:26:16 UTC (rev 361)
@@ -22,7 +22,7 @@
names(x) <- c("time",rownames(from at data))
if (length(from at states)>0) {
nm <- names(x)
- x <- cbind(x,t(from at states[,-1,drop=FALSE]))
+ x <- cbind(x,t(from at states))
names(x) <- c(nm,rownames(from at states))
}
x
@@ -94,19 +94,17 @@
)
object at data[,object at times%in%tt] <- dd[,tt%in%object at times]
if (length(ss)>0) {
- object at states <- array(
+ object at states <- array(
data=NA,
- dim=c(nrow(ss),length(object at times)+1),
+ dim=c(nrow(ss),length(object at times)),
dimnames=list(rownames(ss),NULL)
)
- if (ncol(ss)>length(tt)) tt <- c(tt0,tt)
- nt <- c(object at t0,object at times)
- for (kt in seq_len(length(nt))) {
- wr <- which(nt[kt]==tt)
- if (length(wr)>0)
- object at states[,kt] <- ss[,wr[1]]
- }
- }
+ for (kt in seq_along(object at times)) {
+ wr <- which(object at times[kt]==tt)
+ if (length(wr)>0)
+ object at states[,kt] <- ss[,wr[1]]
+ }
+ }
object
}
)
Modified: pkg/R/simulate-pomp.R
===================================================================
--- pkg/R/simulate-pomp.R 2010-10-04 19:53:57 UTC (rev 360)
+++ pkg/R/simulate-pomp.R 2010-10-05 16:26:16 UTC (rev 361)
@@ -2,11 +2,39 @@
simulate.internal <- function (object, nsim = 1, seed = NULL, params,
states = FALSE, obs = FALSE,
- times = time(object,t0=TRUE), ...) {
+ times, t0, ...) {
+ warn.condition <- (missing(t0) && ((obs || states) || (!missing(times))))
+ if (warn.condition)
+ warning(
+ "The default behavior of ",sQuote("simulate")," has changed.\n",
+ "See the documentation (",dQuote("pomp?simulate"),") for details\n",
+ "and set the ",sQuote("times")," and ",sQuote("t0"),
+ " arguments appropriately to compensate.\n",
+ call.=FALSE
+ )
+
+ if (missing(times)) {
+ times <- time(object,t0=FALSE)
+ } else {
+ times <- as.numeric(times)
+ }
+
+ if (length(times)==0)
+ stop("if ",sQuote("times")," is empty, there is no work to do",call.=FALSE)
+
+ if (any(diff(times)<0))
+ stop(sQuote("times")," must be a nondecreasing sequence of times",call.=FALSE)
+
+ if (missing(t0)) {
+ t0 <- timezero(object)
+ } else {
+ t0 <- as.numeric(t0)
+ }
+
+ if (t0>times[1])
+ stop("the zero-time ",sQuote("t0")," must occur no later than the first observation",call.=FALSE)
ntimes <- length(times)
- times <- as.numeric(times)
- if (ntimes<1)
- stop("if length of ",sQuote("times")," is less than 1, there is no work to do",call.=FALSE)
+
if (missing(params)) {
if (length(object at params)>0)
params <- object at params
@@ -16,61 +44,67 @@
if (is.null(dim(params)))
params <- matrix(params,ncol=1,dimnames=list(names(params),NULL))
npars <- ncol(params)
+
if (!is.null(seed)) { # set the random seed (be very careful about this)
if (!exists('.Random.seed',envir=.GlobalEnv)) runif(1)
save.seed <- get('.Random.seed',envir=.GlobalEnv)
set.seed(seed)
}
+
if (!is.numeric(nsim)||(length(nsim)>1)||nsim<1)
stop(sQuote("nsim")," must be a positive integer")
nsim <- as.integer(nsim)
- xstart <- init.state(object,params=params,t0=times[1])
+
+ xstart <- init.state(object,params=params,t0=t0)
nreps <- npars*nsim # total number of replicates
## we will do the process model simulations with single calls to the user functions
if (nsim > 1) { # make nsim copies of the IC and parameter matrices
xstart <- array(rep(xstart,nsim),dim=c(nrow(xstart),nreps),dimnames=list(rownames(xstart),NULL))
params <- array(rep(params,nsim),dim=c(nrow(params),nreps),dimnames=list(rownames(params),NULL))
}
- x <- rprocess(object,xstart=xstart,times=times,params=params) # simulate the process model
- ## rprocess returns a rank-3 matrix (nvars x nreps x ntimes)
+
+ x <- rprocess(object,xstart=xstart,times=c(t0,times),params=params) # simulate the process model
+ x <- x[,,-1,drop=FALSE]
+ ## rprocess returns a rank-3 matrix (nvars x nreps x (ntimes+1))
+
if (obs || !states) { # we need only simulate the observation process if obs=T or states=F
y <- rmeasure(object,x=x,times=times,params=params)
## rmeasure returns a rank-3 matrix (nobs x nreps x ntimes)
}
+
if (!is.null(seed)) { # restore the RNG state
assign('.Random.seed',save.seed,envir=.GlobalEnv)
}
+
if (!obs && !states) { # both obs=F and states=F, return a list of pomp objects
nm.y <- rownames(y)
nobs <- nrow(y)
retval <- lapply(
- seq(length=nreps),
+ seq_len(nreps),
function (rep) {
po <- as(object,'pomp') # copy the pomp object
- po at t0 <- times[1]
- po at times <- times[-1]
- po at data <- array(0,dim=c(nobs,ntimes-1),dimnames=list(nm.y,NULL))
- po at data[,] <- y[,rep,-1] # replace the data
- po at params <- params[,rep] # store the parameters
- po at states <- array(dim=dim(x)[c(1,3)],dimnames=list(rownames(x),NULL))
- po at states[,] <- x[,rep,] # store the states
+ po at t0 <- t0
+ po at times <- times
+ po at params <- params[,rep]
+ po at data <- array(0,dim=c(nobs,ntimes),dimnames=list(nm.y,NULL))
+ po at data[,] <- y[,rep,]
+ po at states <- array(dim=c(nrow(x),ntimes),dimnames=list(rownames(x),NULL))
+ po at states[,] <- x[,rep,]
po
}
)
if (nreps==1) {
retval <- retval[[1]]
}
- } else {
- if (!obs && states) { # return states only
- retval <- x
- } else if (obs && !states) { # return only observations
- retval <- y
- } else { # return both states and observations
- retval <- list(
- states = x,
- obs = y
- )
- }
+ } else if (!obs && states) { # return states only
+ retval <- x
+ } else if (obs && !states) { # return only observations
+ retval <- y
+ } else { # return both states and observations
+ retval <- list(
+ states = x,
+ obs = y
+ )
}
retval
}
Modified: pkg/R/spect.R
===================================================================
--- pkg/R/spect.R 2010-10-04 19:53:57 UTC (rev 360)
+++ pkg/R/spect.R 2010-10-05 16:26:16 UTC (rev 361)
@@ -98,9 +98,10 @@
nsim=nsim,
seed=seed,
params=params,
- times=time(object,t0=TRUE),
- obs=TRUE
- )[,,-1,drop=FALSE],
+ obs=TRUE,
+ times=time(object,t0=FALSE),
+ t0=timezero(object)
+ ),
silent=FALSE
)
if (inherits(sims,"try-error"))
Modified: pkg/R/traj-match.R
===================================================================
--- pkg/R/traj-match.R 2010-10-04 19:53:57 UTC (rev 360)
+++ pkg/R/traj-match.R 2010-10-05 16:26:16 UTC (rev 361)
@@ -12,12 +12,13 @@
par.est <- as.integer(est)
}
guess <- start[par.est]
+ t0 <- timezero(object)
opt <- optim(
par=guess,
fn=function (x) {
p <- start
p[par.est] <- x
- x <- trajectory(object,params=p)[,,-1,drop=FALSE]
+ x <- trajectory(object,params=p,t0=t0)
d <- dmeasure(
object,
y=data.array(object),
Modified: pkg/R/trajectory-pomp.R
===================================================================
--- pkg/R/trajectory-pomp.R 2010-10-04 19:53:57 UTC (rev 360)
+++ pkg/R/trajectory-pomp.R 2010-10-05 16:26:16 UTC (rev 361)
@@ -3,6 +3,39 @@
setGeneric('trajectory')
trajectory.internal <- function (object, params, times, t0, ...) {
+
+ warn.condition <- missing(t0)
+ if (warn.condition)
+ warning(
+ "The default behavior of ",sQuote("trajectory")," has changed.\n",
+ "See the documentation (",dQuote("pomp?trajectory"),") for details\n",
+ "and set the ",sQuote("times")," and ",sQuote("t0"),
+ " arguments appropriately to compensate.\n",
+ call.=FALSE
+ )
+
+ if (missing(times)) {
+ times <- time(object,t0=FALSE)
+ } else {
+ times <- as.numeric(times)
+ }
+
+ if (length(times)==0)
+ stop("if ",sQuote("times")," is empty, there is no work to do",call.=FALSE)
+
+ if (any(diff(times)<0))
+ stop(sQuote("times")," must be a nondecreasing sequence of times",call.=FALSE)
+
+ if (missing(t0)) {
+ t0 <- timezero(object)
+ } else {
+ t0 <- as.numeric(t0)
+ }
+
+ if (t0>times[1])
+ stop("the zero-time ",sQuote("t0")," must occur no later than the first observation",call.=FALSE)
+ ntimes <- length(times)
+
if (missing(params)) {
params <- coef(object)
if (length(params)==0) {
@@ -26,12 +59,7 @@
stop("pfilter error: ",sQuote("params")," must have rownames",call.=FALSE)
params <- as.matrix(params)
- if (missing(times))
- times <- time(object,t0=TRUE)
- else
- times <- as.numeric(times)
-
- tm <- times[1]
+ tm <- t0
x0 <- init.state(object,params=params,t0=tm)
nm <- rownames(x0)
dim(x0) <- c(nrow(x0),nrep,1)
@@ -59,7 +87,7 @@
X <- try(
deSolve::lsoda(
y=x0[,j,1],
- times=times,
+ times=c(t0,times),
func=function(t,y,parms){
list(
skeleton(
@@ -84,7 +112,7 @@
stop("trajectory error: error in ",sQuote("lsoda"),call.=FALSE)
if (attr(X,'istate')[[1]]!=2)
warning("abnormal exit from ",sQuote("lsoda"),", istate = ",attr(X,'istate'),call.=FALSE)
- x[,j,] <- t(X[,-1])
+ x[,j,] <- t(X[-1,-1])
}
},
unspecified=stop("deterministic skeleton not specified")
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/ricker.rda
===================================================================
(Binary files differ)
Modified: pkg/data/rw2.rda
===================================================================
(Binary files differ)
Modified: pkg/data/verhulst.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2010-10-04 19:53:57 UTC (rev 360)
+++ pkg/inst/ChangeLog 2010-10-05 16:26:16 UTC (rev 361)
@@ -1,5 +1,55 @@
+2010-10-04 kingaa
+
+ * [r360] R/trajectory-pomp.R: - put guts of trajectory calculation
+ into the named function 'trajectory.internal'
+ * [r359] man/euler.Rd: - remove documentation for removed plugins
+ * [r358] DESCRIPTION, NAMESPACE, R/euler.R: - remove the functions
+ 'euler.simulate', 'onestep.simulate', and 'onestep.density',
+ deprecated since version 0.29-1
+ * [r357] DESCRIPTION, R/aaa.R: - add a banner message on
+ attachment.
+ * [r356] R/simulate-pomp.R: - 'simulate.internal' is now the named
+ function for the guts of 'simulate'
+
+2010-10-01 kingaa
+
+ * [r355] DESCRIPTION, NAMESPACE, R/compare-mif.R,
+ R/compare-pmcmc.R, R/compare.mif.R, R/compare.pmcmc.R,
+ R/init-state-pomp.R, R/init.state-pomp.R, R/mif.R, R/nlf-guts.R,
+ R/pmcmc.R, R/pomp-methods.R, R/probe-match.R, R/probe.R,
+ R/spect-match.R, R/spect.R, data/dacca.rda, data/euler.sir.rda,
+ data/gillespie.sir.rda, data/ou2.rda, data/rw2.rda,
+ data/verhulst.rda, inst/NEWS, man/basic-probes.Rd,
+ man/pomp-methods.Rd, man/probe.Rd, man/spect.Rd,
+ tests/ou2-probe.R, tests/ou2-probe.Rout.save,
+ tests/ricker-probe.R, tests/ricker-probe.Rout.save,
+ tests/ricker-spect.R, tests/ricker.R, tests/ricker.Rout.save: -
+ 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
+
+2010-09-30 kingaa
+
+ * [r354] R/nlf-guts.R, R/nlf-lql.R, R/nlf.R: - cosmetology
+ * [r353] man/basic-probes.Rd: - update the man page to include
+ 'probe.ccf'
+
2010-09-29 kingaa
+ * [r352] DESCRIPTION, NAMESPACE, R/basic-probes.R, inst/ChangeLog,
+ inst/doc/advanced_topics_in_pomp.pdf, inst/doc/intro_to_pomp.pdf,
+ inst/doc/ou2-first-mif.rda, inst/doc/ou2-trajmatch.rda,
+ src/probe_ccf.c, tests/ou2-probe.R: - add in a new probe,
+ 'probe.ccf', which computes the CCF of two variables at various
+ lags
* [r351] R/basic-probes.R, src/mat.c, src/pomp_internal.h,
src/pomp_mat.h, src/probe_marginal.c, src/probe_nlar.c: - remove
'mat.c' in favor of inline functions for regression, defined in
Modified: pkg/inst/NEWS
===================================================================
--- pkg/inst/NEWS 2010-10-04 19:53:57 UTC (rev 360)
+++ pkg/inst/NEWS 2010-10-05 16:26:16 UTC (rev 361)
@@ -1,20 +1,41 @@
-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.34-1:
+
+ o The default behaviors of 'simulate' and 'trajectory' with respect to the argument 'times' has changed.
+ There is a new argument, 't0', which specifies the time at which the initial conditions obtain.
+ The default for 't0' is 'timezero(object)'.
+ The default for 'times' is now 'time(object,t0=FALSE)' (it was 'time(object,t0=TRUE)').
+ This change eliminates an infelicity in which 'simulate' and 'trajectory' returned undesired values at the initial time.
+ To aid users in modifying their codes, a warning is now displayed whenever a potentially affected call is made.
+
+ o A banner is now displayed on package attachment.
+ This will be used to report important changes.
+
+ 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'.
+
+ o The functions 'euler.simulate', 'onestep.simulate', and 'onestep.density', deprecated since version 0.29-1, have been removed.
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.
- o New probes recommended by Simon Wood (Nature 466:1102-1104, 2010) have been implemented. In particular, the new function 'probe.marginal' regresses the marginal distributions of actual and simulated data against a reference distribution, returning as probes the regression coefficients. 'probe.marginal' and 'probe.acf' have been coded in C for speed.
+ o Major improvements in speed in the probe-matching routines have been realized by coding the computationally-intensive portions of these calculations in C.
+ o New probes recommended by Simon Wood (Nature 466:1102-1104, 2010) have been implemented.
+ In particular, the new function 'probe.marginal' regresses the marginal distributions of actual and simulated data against a reference distribution, returning as probes the regression coefficients.
+ 'probe.marginal' and 'probe.acf' have been coded in C for speed.
+
NEW IN VERSION 0.32-1:
+
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 361
More information about the pomp-commits
mailing list