[Pomp-commits] r895 - in pkg/pomp: . R inst man src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Mar 20 16:07:36 CET 2014
Author: kingaa
Date: 2014-03-20 16:07:33 +0100 (Thu, 20 Mar 2014)
New Revision: 895
Added:
pkg/pomp/R/minim.R
Removed:
pkg/pomp/tests/ou2-icfit.R
pkg/pomp/tests/ou2-icfit.Rout.save
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/NAMESPACE
pkg/pomp/R/generics.R
pkg/pomp/R/pomp.R
pkg/pomp/R/probe-match.R
pkg/pomp/R/probe.R
pkg/pomp/R/simulate-pomp.R
pkg/pomp/R/traj-match.R
pkg/pomp/inst/NEWS
pkg/pomp/man/probe.Rd
pkg/pomp/man/probed-pomp-methods.Rd
pkg/pomp/man/spect.Rd
pkg/pomp/man/traj-match.Rd
pkg/pomp/src/sir.c
pkg/pomp/tests/abc.R
pkg/pomp/tests/abc.Rout.save
pkg/pomp/tests/bbs-trajmatch.R
pkg/pomp/tests/bbs-trajmatch.Rout.save
pkg/pomp/tests/ou2-probe.R
pkg/pomp/tests/ou2-probe.Rout.save
pkg/pomp/tests/ou2-trajmatch.R
pkg/pomp/tests/ou2-trajmatch.Rout.save
pkg/pomp/tests/ricker-probe.R
pkg/pomp/tests/ricker-probe.Rout.save
pkg/pomp/tests/rw2.Rout.save
pkg/pomp/tests/sir.Rout.save
Log:
- barycentric transformation added for initial conditions in SIR example
- probe matching and trajectory matching are now done using a unified optimizer (in 'minim.R')
- the 'gr' and 'eval.only' arguments have been removed from 'probe.match' and 'traj.match'
- 'probe.match.objfun' and 'traj.match.objfun' are now exported S4 methods
- some changes to the underlying structure of 'probe.matched.pomp' and 'traj.matched.pomp' objects
- access to slots in 'probed.pomp' is now available via the '$' operator
- RNG seed is now stored in 'probed.pomp' objects if set
- in 'pomp', when dprior is unspecified, the default (flat, improper) prior is sought in the 'pomp' package
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2014-03-17 22:02:29 UTC (rev 894)
+++ pkg/pomp/DESCRIPTION 2014-03-20 15:07:33 UTC (rev 895)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.48-3
-Date: 2014-03-17
+Version: 0.49-1
+Date: 2014-03-20
Authors at R: c(person(given=c("Aaron","A."),family="King",
role=c("aut","cre"),email="kingaa at umich.edu"),
person(given=c("Edward","L."),family="Ionides",role=c("aut")),
@@ -18,7 +18,7 @@
)
URL: http://pomp.r-forge.r-project.org
Description: Inference methods for partially-observed Markov processes
-Depends: R(>= 2.14.1), stats, graphics, methods, mvtnorm, subplex, deSolve
+Depends: R(>= 2.14.1), stats, graphics, methods, mvtnorm, subplex, nloptr, deSolve
License: GPL(>= 2)
LazyData: true
BuildVignettes: true
@@ -31,7 +31,7 @@
dmeasure-pomp.R dprocess-pomp.R skeleton-pomp.R
dprior-pomp.R rprior-pomp.R
simulate-pomp.R trajectory-pomp.R plot-pomp.R
- pfilter.R pfilter-methods.R traj-match.R bsmc.R
+ pfilter.R pfilter-methods.R minim.R traj-match.R bsmc.R
mif-class.R particles-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-objfun.R nlf.R
Modified: pkg/pomp/NAMESPACE
===================================================================
--- pkg/pomp/NAMESPACE 2014-03-17 22:02:29 UTC (rev 894)
+++ pkg/pomp/NAMESPACE 2014-03-20 15:07:33 UTC (rev 895)
@@ -37,6 +37,7 @@
importFrom(mvtnorm,dmvnorm,rmvnorm)
importFrom(subplex,subplex)
importFrom(deSolve,ode)
+importFrom(nloptr,nloptr)
exportClasses(
pomp,
@@ -61,6 +62,8 @@
particles,mif,continue,states,trajectory,
pred.mean,pred.var,filter.mean,conv.rec,
bsmc,pmcmc,abc,
+ traj.match.objfun,
+ probe.match.objfun,
spect,probe,probe.match,
traj.match
)
@@ -98,7 +101,6 @@
probe.marginal,
sannbox,
spect.match,
- traj.match.objfun,
pompBuilder,
pompExample
)
Modified: pkg/pomp/R/generics.R
===================================================================
--- pkg/pomp/R/generics.R 2014-03-17 22:02:29 UTC (rev 894)
+++ pkg/pomp/R/generics.R 2014-03-20 15:07:33 UTC (rev 895)
@@ -48,6 +48,7 @@
## deterministic trajectory computation
setGeneric("trajectory",function(object,...)standardGeneric("trajectory"))
## trajectory matching
+setGeneric("traj.match.objfun",function(object,...)standardGeneric("traj.match.objfun"))
setGeneric("traj.match",function(object,...)standardGeneric("traj.match"))
## ABC algorithm functions
@@ -70,6 +71,7 @@
## synthetic likelihood
setGeneric("probe",function(object,probes,...)standardGeneric("probe"))
## probe matching
+setGeneric("probe.match.objfun",function(object,...)standardGeneric("probe.match.objfun"))
setGeneric("probe.match",function(object,...)standardGeneric("probe.match"))
## power spectrum
Added: pkg/pomp/R/minim.R
===================================================================
--- pkg/pomp/R/minim.R (rev 0)
+++ pkg/pomp/R/minim.R 2014-03-20 15:07:33 UTC (rev 895)
@@ -0,0 +1,75 @@
+minim.internal <- function(objfun, start, est, object, method, transform, verbose, ...)
+{
+
+ transform <- as.logical(transform)
+ est <- as.character(est)
+
+ if (length(start)<1)
+ stop(sQuote("start")," must be supplied")
+
+ if (transform) {
+ start <- partrans(object,start,dir="inverse")
+ if (is.null(names(start))||(!all(est%in%names(start))))
+ stop(sQuote("est")," must refer to parameters named in ",
+ sQuote("partrans(object,start,dir=\"inverse\")"))
+ guess <- start[est]
+ } else {
+ if (is.null(names(start))||(!all(est%in%names(start))))
+ stop(sQuote("est")," must refer to parameters named in ",
+ sQuote("start"))
+ guess <- start[est]
+ }
+
+ if (length(est)==0) {
+
+ params <- c(A=3)
+ val <- objfun(guess)
+ conv <- NA
+ evals <- as.integer(c(1,0))
+ msg <- "no optimization performed"
+
+ } else {
+
+ opts <- list(...)
+
+ if (method == 'subplex') {
+ opt <- subplex::subplex(par=guess,fn=objfun,control=opts)
+ } else if (method=="sannbox") {
+ opt <- sannbox(par=guess,fn=objfun,control=opts)
+ } else if (method=="nloptr") {
+ opt <- nloptr(x0=guess,eval_f=objfun,opts=opts)
+ } else {
+ opt <- optim(par=guess,fn=objfun,method=method,control=opts)
+ }
+
+ msg <- as.character(opt$message)
+ val <- opt$value
+
+ if (method == "nloptr") {
+
+ start[est] <- unname(opt$solution)
+ conv <- opt$status
+ evals <- opt$iterations
+
+ } else {
+
+ start[est] <- unname(opt$par)
+ conv <- opt$convergence
+ evals <- opt$counts
+
+ }
+ }
+
+ if (transform)
+ start <- partrans(object,start,dir='forward')
+
+ list(
+ params=start,
+ est=est,
+ transform=transform,
+ value=val,
+ convergence=as.integer(conv),
+ evals=as.integer(evals),
+ msg=msg
+ )
+}
Modified: pkg/pomp/R/pomp.R
===================================================================
--- pkg/pomp/R/pomp.R 2014-03-17 22:02:29 UTC (rev 894)
+++ pkg/pomp/R/pomp.R 2014-03-20 15:07:33 UTC (rev 895)
@@ -118,7 +118,7 @@
}
if (missing(dprior)) { ## by default, use flat improper prior
- dprior <- pomp.fun(f="_pomp_default_dprior")
+ dprior <- pomp.fun(f="_pomp_default_dprior",PACKAGE="pomp")
} else {
dprior <- pomp.fun(f=dprior,PACKAGE=PACKAGE,
proto=quote(dprior(params,log,...)))
Modified: pkg/pomp/R/probe-match.R
===================================================================
--- pkg/pomp/R/probe-match.R 2014-03-17 22:02:29 UTC (rev 894)
+++ pkg/pomp/R/probe-match.R 2014-03-20 15:07:33 UTC (rev 895)
@@ -2,10 +2,8 @@
"probe.matched.pomp",
contains="probed.pomp",
representation=representation(
- start="numeric",
transform="logical",
est="character",
- weights="numeric",
fail.value="numeric",
value="numeric",
evals="integer",
@@ -14,6 +12,8 @@
)
)
+setMethod("$",signature=signature(x="probe.matched.pomp"),function(x, name)slot(x,name))
+
setMethod(
"summary",
"probe.matched.pomp",
@@ -31,13 +31,18 @@
}
)
-probe.match.objfun <- function (object, params, est, probes,
- nsim = 1, seed = NULL, fail.value = NA,
- transform = FALSE, ...) {
-
+pmof.internal <- function (object, params, est, probes,
+ nsim, seed = NULL, fail.value = NA,
+ transform = FALSE, ...)
+{
+ object <- as(object,"pomp")
transform <- as.logical(transform)
+ fail.value <- as.numeric(fail.value)
if (missing(est)) est <- character(0)
- if (!is.character(est)) stop(sQuote("est")," must be a vector of parameter names")
+ est <- as.character(est)
+ if (missing(nsim)) stop(sQuote("nsim")," must be specified")
+ nsim <- as.integer(nsim)
+
if (missing(params)) params <- coef(object)
if ((!is.numeric(params))||(is.null(names(params))))
stop(sQuote("params")," must be a named numeric vector")
@@ -58,8 +63,8 @@
if (nprobes > nsim)
stop(sQuote("nsim"),"(=",nsim,") should be (much) larger than the number of probes (=",nprobes,")")
- obj.fun <- function (par) {
-
+ function (par) {
+
params[par.est.idx] <- par
if (transform)
@@ -79,273 +84,158 @@
ll <- .Call(synth_loglik,simval,datval)
if (is.finite(ll)||is.na(fail.value)) -ll else fail.value
}
-
- obj.fun
}
-probe.match.internal <- function(object, start, est,
- probes, weights,
- nsim, seed,
- method, verbose,
- eval.only, fail.value, transform, ...) {
+setMethod(
+ "probe.match.objfun",
+ signature=signature(object="pomp"),
+ function (object, params, est, probes,
+ nsim, seed = NULL, fail.value = NA,
+ transform = FALSE, ...)
+ pmof.internal(
+ object=object,
+ params=params,
+ est=est,
+ probes=probes,
+ nsim=nsim,
+ seed=seed,
+ fail.value=fail.value,
+ transform=transform,
+ ...
+ )
+ )
+
+setMethod(
+ "probe.match.objfun",
+ signature=signature(object="probed.pomp"),
+ function (object, probes, nsim, seed, ...) {
- transform <- as.logical(transform)
-
- if (eval.only) {
- est <- character(0)
- guess <- numeric(0)
- transform <- FALSE
- } else {
- if (!is.character(est)) stop(sQuote("est")," must be a vector of parameter names")
- if (length(start)<1)
- stop(sQuote("start")," must be supplied if ",sQuote("object")," contains no parameters")
- if (transform) {
- tstart <- partrans(object,start,dir="inverse")
- if (is.null(names(tstart))||(!all(est%in%names(tstart))))
- stop(sQuote("est")," must refer to parameters named in ",sQuote("partrans(object,start,dir=\"inverse\")"))
- guess <- tstart[est]
- } else {
- if (is.null(names(start))||(!all(est%in%names(start))))
- stop(sQuote("est")," must refer to parameters named in ",sQuote("start"))
- guess <- start[est]
- }
- }
-
- obj <- as(object,"pomp")
- coef(obj) <- start
-
- obj.fn <- probe.match.objfun(
- obj,
- est=est,
+ if (missing(probes)) probes <- object at probes
+ if (missing(nsim)) nsim <- nrow(object at simvals)
+ if (missing(seed)) seed <- object at seed
+
+ probe.match.objfun(
+ object=as(object,"pomp"),
probes=probes,
nsim=nsim,
seed=seed,
- fail.value=fail.value,
- transform=transform
+ ...
)
-
-
- if (eval.only) {
-
- val <- obj.fn(guess)
- conv <- NA
- evals <- as.integer(c(1,0))
- msg <- "no optimization performed"
-
- } else {
-
- if (method == 'subplex') {
- opt <- subplex::subplex(par=guess,fn=obj.fn,control=list(...))
- } else if (method=="sannbox") {
- opt <- sannbox(par=guess,fn=obj.fn,control=list(...))
- } else {
- opt <- optim(par=guess,fn=obj.fn,method=method,control=list(...))
- }
-
- if (!is.null(names(opt$par)) && !all(est==names(opt$par)))
- stop("mismatch between parameter names returned by optimizer and ",sQuote("est"))
- coef(obj,est,transform=transform) <- unname(opt$par)
- msg <- if (is.null(opt$message)) character(0) else opt$message
- conv <- opt$convergence
- val <- opt$value
- evals <- opt$counts
- }
-
- new(
- "probe.matched.pomp",
- probe(
- obj,
- probes=probes,
- nsim=nsim,
- seed=seed
- ),
- start=start,
- transform=transform,
- est=as.character(est),
- value=val,
- convergence=as.integer(conv),
- evals=as.integer(evals),
- msg=as.character(msg)
- )
-}
-
+ }
+ )
+
setMethod(
"probe.match",
signature=signature(object="pomp"),
function(object, start, est = character(0),
- probes, weights,
- nsim, seed = NULL,
- method = c("subplex","Nelder-Mead","SANN","BFGS","sannbox"),
+ probes, nsim, seed = NULL,
+ method = c("subplex","Nelder-Mead","SANN","BFGS",
+ "sannbox","nloptr"),
verbose = getOption("verbose"),
- eval.only = FALSE, fail.value = NA, transform = FALSE,
+ fail.value = NA,
+ transform = FALSE,
...) {
-
- transform <- as.logical(transform)
+
if (missing(start)) start <- coef(object)
if (missing(probes)) stop(sQuote("probes")," must be supplied")
if (missing(nsim)) stop(sQuote("nsim")," must be supplied")
- if (missing(weights)) weights <- 1
method <- match.arg(method)
+ est <- as.character(est)
+ transform <- as.logical(transform)
+ fail.value <- as.numeric(fail.value)
- probe.match.internal(
- object=object,
- start=start,
- est=est,
- probes=probes,
- weights=weights,
- nsim=nsim,
- seed=seed,
- method=method,
- verbose=verbose,
- eval.only=eval.only,
- fail.value=fail.value,
- transform=transform,
- ...
- )
+ m <- minim.internal(
+ objfun=probe.match.objfun(
+ object=object,
+ params=start,
+ est=est,
+ probes=probes,
+ nsim=nsim,
+ seed=seed,
+ fail.value=fail.value,
+ transform=transform
+ ),
+ start=start,
+ est=est,
+ object=object,
+ method=method,
+ transform=transform,
+ verbose=verbose,
+ ...
+ )
+
+ coef(object) <- m$params
+
+ new(
+ "probe.matched.pomp",
+ probe(
+ object,
+ probes=probes,
+ nsim=nsim,
+ seed=seed
+ ),
+ transform=transform,
+ est=est,
+ fail.value=fail.value,
+ value=m$value,
+ evals=m$evals,
+ convergence=m$convergence,
+ msg=m$msg
+ )
}
)
setMethod(
"probe.match",
signature=signature(object="probed.pomp"),
- function(object, start, est = character(0),
- probes, weights,
- nsim, seed = NULL,
- method = c("subplex","Nelder-Mead","SANN","BFGS","sannbox"),
- verbose = getOption("verbose"),
- eval.only = FALSE, fail.value = NA, transform = FALSE, ...) {
-
- transform <- as.logical(transform)
+ function(object, probes, nsim, seed, ...,
+ verbose = getOption("verbose"))
+ {
- if (missing(start)) start <- coef(object)
if (missing(probes)) probes <- object at probes
if (missing(nsim)) nsim <- nrow(object at simvals)
- if (missing(weights)) weights <- 1
-
- method <- match.arg(method)
+ if (missing(seed)) seed <- object at seed
- probe.match.internal(
- object=object,
- start=start,
- est=est,
- probes=probes,
- weights=weights,
- nsim=nsim,
- seed=seed,
- method=method,
- verbose=verbose,
- eval.only=eval.only,
- fail.value=fail.value,
- transform=transform,
- ...
- )
+ f <- selectMethod("probe.match","pomp")
+
+ f(
+ object=object,
+ probes=probes,
+ nsim=nsim,
+ seed=seed,
+ verbose=verbose,
+ ...
+ )
}
)
setMethod(
"probe.match",
signature=signature(object="probe.matched.pomp"),
- function(object, start, est,
- probes, weights,
- nsim, seed = NULL,
- method = c("subplex","Nelder-Mead","SANN","BFGS","sannbox"),
- verbose = getOption("verbose"),
- eval.only = FALSE, fail.value, transform, ...) {
-
- if (missing(start)) start <- coef(object)
+ function(object, est, probes, nsim, seed, transform,
+ fail.value, ..., verbose = getOption("verbose"))
+ {
+
if (missing(est)) est <- object at est
+ if (missing(probes)) probes <- object at probes
+ if (missing(nsim)) nsim <- nrow(object at simvals)
+ if (missing(seed)) seed <- object at seed
if (missing(transform)) transform <- object at transform
-
- if (missing(probes))
- probes <- object at probes
-
- if (missing(nsim))
- nsim <- nrow(object at simvals)
-
- if (missing(weights)) weights <- 1
-
if (missing(fail.value)) fail.value <- object at fail.value
-
- method <- match.arg(method)
- probe.match.internal(
- object=object,
- start=start,
- est=est,
- probes=probes,
- weights=weights,
- nsim=nsim,
- seed=seed,
- method=method,
- verbose=verbose,
- eval.only=eval.only,
- fail.value=fail.value,
- transform=transform,
- ...
- )
+ f <- selectMethod("probe.match","pomp")
+
+ f(
+ object=object,
+ est=est,
+ probes=probes,
+ nsim=nsim,
+ seed=seed,
+ transform=transform,
+ fail.value=fail.value,
+ verbose=verbose,
+ ...
+ )
}
)
-
-probe.mismatch <- function (par, est, object, probes, params,
- nsim = 1, seed = NULL,
- weights, datval,
- fail.value = NA) {
- if (missing(par)) par <- numeric(0)
- if (missing(est)) est <- integer(0)
- if (missing(params)) params <- coef(object)
-
- params[est] <- par
-
- ## apply probes to model simulations
- simval <- .Call(
- apply_probe_sim,
- object=object,
- nsim=nsim,
- params=params,
- seed=seed,
- probes=probes,
- datval=datval
- )
-
- ## compute a measure of the discrepancies between simulations and data
- sim.means <- colMeans(simval)
- simval <- sweep(simval,2,sim.means)
- discrep <- ((datval-sim.means)^2)/colMeans(simval^2)
- if ((length(weights)>1) && (length(weights)!=length(discrep)))
- stop(length(discrep)," probes have been computed, but ",length(weights)," have been supplied")
- if (!all(is.finite(discrep))) {
- mismatch <- fail.value
- } else if (length(weights)>1) {
- mismatch <- sum(discrep*weights)/sum(weights)
- } else {
- mismatch <- sum(discrep)
- }
-
- mismatch
-}
-
-neg.synth.loglik <- function (par, est, object, probes, params,
- nsim = 1, seed = NULL,
- weights, datval,
- fail.value = NA) {
- if (missing(par)) par <- numeric(0)
- if (missing(est)) est <- integer(0)
- if (missing(params)) params <- coef(object)
-
- params[est] <- par
-
- ## apply probes to model simulations
- simval <- .Call(
- apply_probe_sim,
- object=object,
- nsim=nsim,
- params=params,
- seed=seed,
- probes=probes,
- datval=datval
- )
-
- ll <- .Call(synth_loglik,simval,datval)
- if (is.finite(ll)||is.na(fail.value)) -ll else fail.value
-}
Modified: pkg/pomp/R/probe.R
===================================================================
--- pkg/pomp/R/probe.R 2014-03-17 22:02:29 UTC (rev 894)
+++ pkg/pomp/R/probe.R 2014-03-20 15:07:33 UTC (rev 895)
@@ -7,16 +7,20 @@
simvals="array",
quantiles="numeric",
pvals="numeric",
- synth.loglik="numeric"
+ synth.loglik="numeric",
+ seed="integer"
)
)
probe.internal <- 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")
+
+ seed <- as.integer(seed)
if (missing(params)) params <- coef(object)
@@ -25,6 +29,7 @@
nprobes <- length(datval)
if (nprobes > nsim)
stop(sQuote("nsim"),"(=",nsim,") should be (much) larger than the number of probes (=",nprobes,")")
+
## apply probes to model simulations
simval <- .Call(
apply_probe_sim,
@@ -59,29 +64,41 @@
simvals=simval,
quantiles=quants,
pvals=pvals,
- synth.loglik=ll
+ synth.loglik=ll,
+ seed=seed
)
}
-setMethod("probe",signature(object="pomp"),
- function (object, probes, params, nsim = 1, seed = NULL, ...) {
- probe.internal(object=object,probes=probes,params=params,
- nsim=nsim,seed=seed,...)
+setMethod(
+ "probe",
+ signature=signature(object="pomp"),
+ function (object, probes, params, nsim = 1, seed = NULL, ...)
+ {
+ probe.internal(
+ object=object,
+ probes=probes,
+ params=params,
+ nsim=nsim,
+ seed=seed,
+ ...
+ )
}
)
setMethod(
"probe",
- signature(object="probed.pomp"),
- function (object, probes, params, nsim, ...) {
+ signature=signature(object="probed.pomp"),
+ function (object, probes, params, nsim, seed, ...) {
if (missing(probes)) probes <- object at probes
if (missing(nsim)) nsim <- nrow(object at simvals)
+ if (missing(seed)) seed <- object at seed
probe(
- as(object,"pomp"),
+ object=as(object,"pomp"),
probes=probes,
nsim=nsim,
+ seed=seed,
...
)
}
@@ -187,3 +204,4 @@
)
setMethod("logLik",signature(object="probed.pomp"),function(object,...)object at synth.loglik)
+setMethod("$",signature=signature(x="probed.pomp"),function(x, name)slot(x,name))
Modified: pkg/pomp/R/simulate-pomp.R
===================================================================
--- pkg/pomp/R/simulate-pomp.R 2014-03-17 22:02:29 UTC (rev 894)
+++ pkg/pomp/R/simulate-pomp.R 2014-03-20 15:07:33 UTC (rev 895)
@@ -27,7 +27,8 @@
params <- as.matrix(params)
- if (!is.null(seed)) { # set the random seed (be very careful about this)
+ ## set the random seed (be very careful about this)
+ if (!is.null(seed) && length(seed)>0) {
if (!exists('.Random.seed',envir=.GlobalEnv)) runif(1)
save.seed <- get('.Random.seed',envir=.GlobalEnv)
set.seed(seed)
@@ -55,7 +56,8 @@
if (inherits(retval,'try-error'))
stop(sQuote("simulate")," error",call.=FALSE)
- if (!is.null(seed)) { # restore the RNG state
+ ## restore the RNG state
+ if (!is.null(seed) && length(seed)>0) {
assign('.Random.seed',save.seed,envir=.GlobalEnv)
}
@@ -113,11 +115,22 @@
retval
}
-setMethod("simulate",signature=signature(object="pomp"),
+setMethod(
+ "simulate",
+ signature=signature(object="pomp"),
definition=function (object, nsim = 1, seed = NULL, params,
states = FALSE, obs = FALSE,
times, t0, as.data.frame = FALSE, ...)
- simulate.internal(object=object,nsim=nsim,seed=seed,params=params,
- states=states,obs=obs,times=times,t0=t0,
- as.data.frame=as.data.frame,...)
+ simulate.internal(
+ object=object,
+ nsim=nsim,
+ seed=seed,
+ params=params,
+ states=states,
+ obs=obs,
+ times=times,
+ t0=t0,
+ as.data.frame=as.data.frame,
+ ...
+ )
)
Modified: pkg/pomp/R/traj-match.R
===================================================================
--- pkg/pomp/R/traj-match.R 2014-03-17 22:02:29 UTC (rev 894)
+++ pkg/pomp/R/traj-match.R 2014-03-20 15:07:33 UTC (rev 895)
@@ -2,7 +2,6 @@
"traj.matched.pomp",
contains="pomp",
representation=representation(
- start="numeric",
transform="logical",
est="character",
evals="integer",
@@ -32,11 +31,13 @@
}
)
-traj.match.objfun <- function (object, params, est, transform = FALSE, ...) {
+tmof.internal <- function (object, params, est, transform, ...) {
- transform <- as.logical(transform)
+ object <- as(object,"pomp")
if (missing(est)) est <- character(0)
- if (!is.character(est)) stop(sQuote("est")," must be a vector of parameter names")
+ est <- as.character(est)
+ transform <- as.logical(transform)
+
if (missing(params)) params <- coef(object)
if ((!is.numeric(params))||(is.null(names(params))))
stop(sQuote("params")," must be a named numeric vector")
@@ -46,156 +47,105 @@
if (any(is.na(par.est.idx)))
stop("parameter(s): ",sQuote(est[is.na(par.est.idx)])," not found in ",sQuote("params"))
- obj.fn <- function (par) {
+ function (par) {
params[par.est.idx] <- par
- if (transform) {
+ if (transform)
tparams <- partrans(object,params,dir="forward")
- d <- dmeasure(
+ d <- dmeasure(
+ object,
+ y=object at data,
+ x=trajectory(
object,
- y=object at data,
- x=trajectory(object,params=tparams,...),
- times=time(object),
- params=tparams,
- log=TRUE
- )
- } else {
- d <- dmeasure(
- object,
- y=object at data,
- x=trajectory(object,params=params,...),
- times=time(object),
- params=params,
- log=TRUE
- )
- }
+ params=if (transform) tparams else params,
+ ...
+ ),
+ times=time(object),
+ params=if (transform) tparams else params,
+ log=TRUE
+ )
-sum(d)
}
-
- obj.fn
}
-traj.match.internal <- function (object, start, est, method, gr, eval.only, transform, ...) {
-
- transform <- as.logical(transform)
+setMethod(
+ "traj.match.objfun",
+ signature=signature(object="pomp"),
+ function (object, params, est, transform = FALSE, ...)
+ tmof.internal(
+ object=object,
+ params=params,
+ est=est,
+ transform=transform,
+ ...
+ )
+ )
- if (eval.only) {
- est <- character(0)
- guess <- numeric(0)
- transform <- FALSE
- } else {
- if (!is.character(est)) stop(sQuote("est")," must be a vector of parameter names")
- if (length(start)<1)
- stop(sQuote("start")," must be supplied if ",sQuote("object")," contains no parameters")
- if (transform) {
- tstart <- partrans(object,start,dir="inverse")
- if (is.null(names(tstart))||(!all(est%in%names(tstart))))
- stop(sQuote("est")," must refer to parameters named in ",sQuote("partrans(object,start,dir=\"inverse\")"))
- guess <- tstart[est]
- } else {
- if (is.null(names(start))||(!all(est%in%names(start))))
- stop(sQuote("est")," must refer to parameters named in ",sQuote("start"))
- guess <- start[est]
- }
- }
-
- obj <- as(object,"pomp")
- coef(obj) <- start
-
- obj.fn <- traj.match.objfun(obj,est=est,transform=transform)
-
- if (eval.only) {
-
- val <- obj.fn(guess)
- conv <- NA
- evals <- c(1,0)
- msg <- "no optimization performed"
-
- } else {
-
- if (method=="subplex") {
- opt <- subplex::subplex(par=guess,fn=obj.fn,control=list(...))
- } else if (method=="sannbox") {
- opt <- sannbox(par=guess,fn=obj.fn,control=list(...))
- } else {
- opt <- optim(par=guess,fn=obj.fn,gr=gr,method=method,control=list(...))
- }
-
- if (!is.null(names(opt$par)) && !all(est==names(opt$par)))
- stop("mismatch between parameter names returned by optimizer and ",sQuote("est"))
- coef(obj,est,transform=transform) <- unname(opt$par)
- msg <- if (is.null(opt$message)) character(0) else opt$message
- conv <- opt$convergence
- evals <- opt$counts
- val <- opt$value
-
- }
-
- ## fill 'states' slot of returned object with the trajectory
- x <- trajectory(obj)
- obj at states <- array(data=x,dim=dim(x)[c(1L,3L)])
- rownames(obj at states) <- rownames(x)
-
- new(
- "traj.matched.pomp",
- obj,
- start=start,
- transform=transform,
- est=as.character(est),
- evals=as.integer(evals),
- convergence=as.integer(conv),
- msg=msg,
- value=as.numeric(-val)
- )
-}
-
-traj.match <- function (object, ...)
- stop("function ",sQuote("traj.match")," is undefined for objects of class ",sQuote(class(object)))
-
setMethod(
"traj.match",
signature=signature(object="pomp"),
- function (object, start, est,
- method = c("Nelder-Mead","subplex","SANN","BFGS","sannbox"),
- gr = NULL, eval.only = FALSE, transform = FALSE, ...) {
- transform <- as.logical(transform)
+ function (object, start, est = character(0),
+ method = c("Nelder-Mead","subplex","SANN","BFGS",
+ "sannbox","nloptr"),
+ transform = FALSE, ...)
+ {
+
if (missing(start)) start <- coef(object)
- if (!eval.only && missing(est))
- stop(sQuote("est")," must be supplied if optimization is to be done")
- if (eval.only) est <- character(0)
+
method <- match.arg(method)
- traj.match.internal(
- object=object,
+ est <- as.character(est)
+ transform <- as.logical(transform)
+
+ m <- minim.internal(
+ objfun=traj.match.objfun(
+ object=object,
+ params=start,
+ est=est,
+ transform=transform
+ ),
start=start,
est=est,
+ object=object,
method=method,
- gr=gr,
- eval.only=eval.only,
transform=transform,
...
)
+
+ ## fill params slot appropriately
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 895
More information about the pomp-commits
mailing list