[Pomp-commits] r405 - in pkg: . R man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Nov 5 16:06:37 CET 2010
Author: kingaa
Date: 2010-11-05 16:06:37 +0100 (Fri, 05 Nov 2010)
New Revision: 405
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/aaa.R
pkg/R/probe-match.R
pkg/R/spect-match.R
pkg/R/traj-match.R
pkg/man/traj-match.Rd
pkg/tests/ou2-trajmatch.R
pkg/tests/ou2-trajmatch.Rout.save
Log:
- if eval.only=TRUE in traj.match, spect.match, or probe.match, the returned convergence field is now NA and the message says "no optimization performed"
- traj.match now returns an object of class traj.matched.pomp, a new class that inherits from 'pomp'
- the 'states' slot in the 'traj.matched.pomp' object is now filled with the trajectory
- more optimization methods are now provided in 'traj.match'
- the only class now exported is 'pomp'.
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-10-25 11:02:47 UTC (rev 404)
+++ pkg/DESCRIPTION 2010-11-05 15:06:37 UTC (rev 405)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.35-1
-Date: 2010-10-25
+Version: 0.36-1
+Date: 2010-11-05
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>
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2010-10-25 11:02:47 UTC (rev 404)
+++ pkg/NAMESPACE 2010-11-05 15:06:37 UTC (rev 405)
@@ -33,12 +33,12 @@
importFrom(deSolve,lsoda)
importFrom(subplex,subplex)
-exportClasses('pomp','mif','pmcmc','probed.pomp','probe.matched.pomp','spect.pomp','spect.matched.pomp')
+exportClasses('pomp')
exportMethods(
'plot','show','print','coerce',
'dprocess','rprocess','rmeasure','dmeasure','init.state','skeleton',
- 'data.array','obs','coef','time','time<-','timezero','timezero<-',
+ 'data.array','obs','coef','time','time<-','timezero','timezero<-','$',
'simulate','pfilter',
'particles','mif','continue','coef<-','states','trajectory',
'summary','logLik','window',
Modified: pkg/R/aaa.R
===================================================================
--- pkg/R/aaa.R 2010-10-25 11:02:47 UTC (rev 404)
+++ pkg/R/aaa.R 2010-11-05 15:06:37 UTC (rev 405)
@@ -5,3 +5,4 @@
## cat("This is pomp version ",version,"\n\n",sep="")
## cat("See the NEWS file for important information\n")
##}
+
Modified: pkg/R/probe-match.R
===================================================================
--- pkg/R/probe-match.R 2010-10-25 11:02:47 UTC (rev 404)
+++ pkg/R/probe-match.R 2010-11-05 15:06:37 UTC (rev 405)
@@ -153,7 +153,7 @@
datval=datval,
fail.value=fail.value
)
- conv <- 0
+ conv <- NA
evals <- as.integer(c(1,0))
msg <- paste("no optimization performed")
} else {
Modified: pkg/R/spect-match.R
===================================================================
--- pkg/R/spect-match.R 2010-10-25 11:02:47 UTC (rev 404)
+++ pkg/R/spect-match.R 2010-11-05 15:06:37 UTC (rev 405)
@@ -173,9 +173,9 @@
data.spec=ds,
fail.value=fail.value
)
- conv <- 0
- evals <- as.integer(1)
- msg <- paste(sQuote("spec.mismatch"),"evaluated")
+ conv <- NA
+ evals <- as.integer(c(1,0))
+ msg <- paste("no optimization performed")
} else {
if (method == 'subplex') {
opt <- subplex::subplex(
Modified: pkg/R/traj-match.R
===================================================================
--- pkg/R/traj-match.R 2010-10-25 11:02:47 UTC (rev 404)
+++ pkg/R/traj-match.R 2010-11-05 15:06:37 UTC (rev 405)
@@ -1,41 +1,149 @@
-traj.match <- function (object, start, est, method = "Nelder-Mead", gr = NULL, ...) {
+setClass(
+ "traj.matched.pomp",
+ contains="pomp",
+ representation=representation(
+ evals="integer",
+ convergence="integer",
+ msg="character",
+ value="numeric"
+ )
+ )
+
+setMethod("$",signature(x="traj.matched.pomp"),function(x, name) slot(x,name))
+
+setMethod("logLik",signature(object="traj.matched.pomp"),function(object,...)object at value)
+
+setMethod(
+ "summary",
+ "traj.matched.pomp",
+ function (object, ...) {
+ c(
+ list(
+ params=coef(object),
+ loglik=object at value,
+ eval=object at evals,
+ convergence=object at convergence
+ ),
+ if(length(object at msg)>0) list(msg=object at msg) else NULL
+ )
+ }
+ )
+
+trajectory.nll <- function (par, est, object, params, t0, fail.value = NA, ...) {
+ if (missing(par)) par <- numeric(0)
+ if (missing(est)) est <- integer(0)
+ if (missing(params)) params <- coef(object)
+ if (missing(t0)) t0 <- timezero(object)
+ params[est] <- par
+ x <- trajectory(object,params=params,t0=t0)
+ d <- dmeasure(
+ object,
+ y=data.array(object),
+ x=x,
+ times=time(object),
+ params=as.matrix(params),
+ log=TRUE
+ )
+ -sum(d)
+}
+
+trajectory.ls <- function (par, est, object, params, t0, fail.value = NA, transform, ...) {
+ if (missing(transform)) transform <- identity
+ if (missing(par)) par <- numeric(0)
+ if (missing(est)) est <- integer(0)
+ if (missing(params)) params <- coef(object)
+ if (missing(t0)) t0 <- timezero(object)
+ params[est] <- par
+ x <- trajectory(object,params=params,t0=t0)
+ d <- dmeasure(
+ object,
+ y=data.array(object),
+ x=x,
+ times=time(object),
+ params=as.matrix(params),
+ log=TRUE
+ )
+ -sum(d)
+}
+
+traj.match <- function (object, start, est, method = c("Nelder-Mead","SANN","subplex"),
+ gr = NULL, eval.only = FALSE, ...) {
+
if (!is(object,'pomp'))
stop("traj.match error: ",sQuote("object")," must be a ",sQuote("pomp")," object",call.=FALSE)
- if (is.character(est)) {
- if (!all(est%in%names(start)))
- stop("traj.match error: parameters named in ",sQuote("est")," must exist in ",sQuote("start"),call.=FALSE)
- par.est <- which(names(start)%in%est)
- } else if (is.numeric(est)) {
- est <- as.integer(est)
- if (any((est<1)|(est>length(start))))
- stop("traj.match error: some index in ",sQuote("est")," corresponds to no parameters in ",sQuote("start"),call.=FALSE)
- par.est <- as.integer(est)
- }
+
+ if (missing(start)) start <- coef(object)
+
+ method <- match.arg(method)
+
+ if (!is.character(est)) stop(sQuote("est")," must be a vector of parameter names")
+ if (!all(est%in%names(start)))
+ stop("traj.match error: parameters named in ",sQuote("est")," must exist in ",sQuote("start"),call.=FALSE)
+ par.est <- which(names(start)%in%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,t0=t0)
- d <- dmeasure(
- object,
- y=data.array(object),
- x=x,
- times=time(object),
- params=as.matrix(p),
- log=TRUE
- )
- -sum(d)
- },
- gr=gr,
- method=method,
- control=list(...)
- )
- opt$value <- -opt$value
- start[par.est] <- opt$par
- opt$params <- start
- opt$par <- NULL
- opt
+ obj <- as(object,"pomp")
+
+ obj.fn <- function (x) {
+ p <- start
+ p[par.est] <- x
+ x <- trajectory(obj,params=p,t0=t0)
+ d <- dmeasure(
+ obj,
+ y=data.array(object),
+ x=trajectory(obj,params=p,t0=t0),
+ times=time(object),
+ params=as.matrix(p),
+ log=TRUE
+ )
+ -sum(d)
+ }
+
+ if (eval.only) {
+
+ coef(obj,names(start)) <- unname(start)
+ val <- obj.fn(guess)
+ conv <- NA
+ evals <- c(1,0)
+ msg <- paste("no optimization performed")
+
+ } else {
+
+ if (method=="subplex") {
+
+ opt <- subplex::subplex(
+ par=guess,
+ fn=obj.fn,
+ control=list(...)
+ )
+
+ } else {
+
+ opt <- optim(
+ par=guess,
+ fn=obj.fn,
+ gr=gr,
+ method=method,
+ control=list(...)
+ )
+ }
+
+ coef(obj,names(opt$par)) <- unname(opt$par)
+ obj at states[,] <- trajectory(obj,t0=t0)[,1,]
+ msg <- if (is.null(opt$message)) character(0) else opt$message
+ conv <- opt$convergence
+ evals <- opt$counts
+ val <- opt$value
+
+ }
+
+ ans <- new(
+ "traj.matched.pomp",
+ obj,
+ evals=as.integer(evals),
+ convergence=as.integer(conv),
+ msg=msg,
+ value=as.numeric(-val)
+ )
}
Modified: pkg/man/traj-match.Rd
===================================================================
--- pkg/man/traj-match.Rd 2010-10-25 11:02:47 UTC (rev 404)
+++ pkg/man/traj-match.Rd 2010-11-05 15:06:37 UTC (rev 405)
@@ -1,30 +1,47 @@
\name{traj.match}
\alias{traj.match}
+\alias{logLik,traj.matched.pomp-method}
+\alias{logLik-traj.matched.pomp}
+\alias{$,traj.matched.pomp-method}
+\alias{$-traj.matched.pomp}
+\alias{summary,traj.matched.pomp-method}
+\alias{summary-traj.matched.pomp}
\title{Trajectory matching}
\description{
- Match trajectories.
+ Match trajectories to data.
}
\usage{
-traj.match(object, start, est, method = "Nelder-Mead", gr = NULL, \dots)
+traj.match(object, start, est, method = c("Nelder-Mead", "SANN", "subplex"), gr = NULL,
+ eval.only = FALSE, \dots)
}
\arguments{
\item{object}{A \code{pomp} object.}
\item{start}{Initial guess for parameters.}
\item{est}{Character vector containing the names of parameters to be estimated.}
\item{method}{
- One of the optimization methods recognized by \code{\link{optim}}.
+ Optimization method.
+ Choices are \code{\link[subplex]{subplex}} and any of the methods used by \code{\link{optim}}.
}
\item{gr}{
Passed to \code{\link{optim}}.
}
+ \item{eval.only}{
+ logical; if \code{TRUE}, no optimization is attempted and the quasi-loglikelihood value is evaluated at the \code{start} parameters.
+ }
\item{\dots}{
- Arguments that will be passed to \code{\link{optim}} in the \code{control} list.
+ Arguments that will be passed to \code{\link{optim}} in its \code{control} list.
}
}
\details{
Trajectory matching is accomplished using \code{\link{optim}}.
The \code{\link{trajectory}} method is used for this, which in turn uses the \code{skeleton} slot of the \code{pomp} object.
+ The quantity maximized is the likelihood of the data given the trajectory, as returned by \code{\link{dmeasure}}.
}
+\value{
+ An object of class \code{traj.matched.pomp}.
+ Available methods for objects of this type include \code{summary} and \code{logLik}.
+ The other slots of this object can be accessed via the \code{$} operator.
+}
\examples{
data(ou2)
true.p <- c(
@@ -44,6 +61,13 @@
reltol=1e-8
)
+ summary(res)
+
+ plot(range(time(res)),range(c(obs(res),states(res))),type='n',xlab="time",ylab="x,y")
+ points(y1~time,data=as(res,"data.frame"),col='blue')
+ points(y2~time,data=as(res,"data.frame"),col='red')
+ lines(x1~time,data=as(res,"data.frame"),col='blue')
+ lines(x2~time,data=as(res,"data.frame"),col='red')
}
\seealso{
\code{\link{optim}},
Modified: pkg/tests/ou2-trajmatch.R
===================================================================
--- pkg/tests/ou2-trajmatch.R 2010-10-25 11:02:47 UTC (rev 404)
+++ pkg/tests/ou2-trajmatch.R 2010-11-05 15:06:37 UTC (rev 405)
@@ -16,7 +16,7 @@
maxit=2000,
reltol=1e-8
)
- c(conv=res$convergence,loglik=res$value,res$params)
+ c(conv=res$convergence,loglik=logLik(res),coef(res))
}
)
range(x['conv',])
@@ -30,3 +30,8 @@
),
digits=4
)
+
+summary(traj.match(ou2,est=c('alpha.1','alpha.4','x1.0','x2.0','tau'),method="subplex",maxit=100))
+
+summary(traj.match(ou2,est=c('alpha.1','alpha.4','x1.0','x2.0','tau'),eval.only=T))
+
Modified: pkg/tests/ou2-trajmatch.Rout.save
===================================================================
--- pkg/tests/ou2-trajmatch.Rout.save 2010-10-25 11:02:47 UTC (rev 404)
+++ pkg/tests/ou2-trajmatch.Rout.save 2010-11-05 15:06:37 UTC (rev 405)
@@ -33,7 +33,7 @@
+ maxit=2000,
+ reltol=1e-8
+ )
-+ c(conv=res$convergence,loglik=res$value,res$params)
++ c(conv=res$convergence,loglik=logLik(res),coef(res))
+ }
+ )
> range(x['conv',])
@@ -54,10 +54,49 @@
alpha.2 -0.5 -0.5000 0.00000 0.0000
alpha.3 0.3 0.3000 0.00000 0.0000
alpha.4 0.9 0.8318 -0.06822 0.4922
-sigma.1 3.0 0.0000 -3.00000 0.0000
-sigma.2 -0.5 0.0000 0.50000 0.0000
-sigma.3 2.0 0.0000 -2.00000 0.0000
+sigma.1 3.0 3.0000 0.00000 0.0000
+sigma.2 -0.5 -0.5000 0.00000 0.0000
+sigma.3 2.0 2.0000 0.00000 0.0000
tau 1.0 6.0282 5.02824 1.4802
x1.0 -3.0 -3.3701 -0.37012 2.7226
x2.0 4.0 2.3504 -1.64964 4.4465
>
+> summary(traj.match(ou2,est=c('alpha.1','alpha.4','x1.0','x2.0','tau'),method="subplex",maxit=100))
+$params
+ alpha.1 alpha.2 alpha.3 alpha.4 sigma.1 sigma.2 sigma.3
+ 0.8913424 -0.5000000 0.3000000 0.2309642 3.0000000 -0.5000000 2.0000000
+ tau x1.0 x2.0
+ 5.3538853 -2.6942832 9.5177620
+
+$loglik
+[1] -619.7905
+
+$eval
+[1] 100
+
+$convergence
+[1] -1
+
+$msg
+[1] "number of function evaluations exceeds `maxit'"
+
+>
+> summary(traj.match(ou2,est=c('alpha.1','alpha.4','x1.0','x2.0','tau'),eval.only=T))
+$params
+alpha.1 alpha.2 alpha.3 alpha.4 sigma.1 sigma.2 sigma.3 tau x1.0 x2.0
+ 0.8 -0.5 0.3 0.9 3.0 -0.5 2.0 1.0 -3.0 4.0
+
+$loglik
+[1] -3327.081
+
+$eval
+[1] 1 0
+
+$convergence
+[1] NA
+
+$msg
+[1] "no optimization performed"
+
+>
+>
More information about the pomp-commits
mailing list