[Pomp-commits] r155 - in pkg: . R data man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jul 7 15:35:40 CEST 2009
Author: kingaa
Date: 2009-07-07 15:35:40 +0200 (Tue, 07 Jul 2009)
New Revision: 155
Added:
pkg/R/trajmatch.R
pkg/man/trajmatch.Rd
pkg/tests/ou2-trajmatch.R
pkg/tests/ou2-trajmatch.Rout.save
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/data/euler.sir.rda
pkg/data/ou2.rda
pkg/data/rw2.rda
pkg/tests/ou2-mif.Rout.save
Log:
add trajectory-matching method
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2009-06-24 15:48:22 UTC (rev 154)
+++ pkg/DESCRIPTION 2009-07-07 13:35:40 UTC (rev 155)
@@ -1,11 +1,11 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.25-2
-Date: 2009-06-24
+Version: 0.25-3
+Date: 2009-07-07
Author: Aaron A. King, Edward L. Ionides, Carles Martinez Breto, Steve Ellner, Bruce Kendall
Maintainer: Aaron A. King <kingaa at umich.edu>
Description: Inference methods for partially-observed Markov processes
Depends: R(>= 2.8.1), stats, methods, graphics, deSolve, subplex, mvtnorm
-License: GPL (>= 2)
+License: GPL(>= 2)
LazyLoad: true
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2009-06-24 15:48:22 UTC (rev 154)
+++ pkg/NAMESPACE 2009-07-07 13:35:40 UTC (rev 155)
@@ -42,5 +42,6 @@
bspline.basis,
periodic.bspline.basis,
compare.mif,
- nlf
+ nlf,
+ traj.match
)
Added: pkg/R/trajmatch.R
===================================================================
--- pkg/R/trajmatch.R (rev 0)
+++ pkg/R/trajmatch.R 2009-07-07 13:35:40 UTC (rev 155)
@@ -0,0 +1,40 @@
+traj.match <- function (object, start, est, method = "Nelder-Mead", gr = NULL, ...) {
+ 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)
+ }
+ guess <- start[par.est]
+ opt <- optim(
+ par=guess,
+ fn=function (x) {
+ p <- start
+ p[par.est] <- x
+ x <- simulate(object,nsim=1,params=p,states=TRUE)[,,-1,drop=FALSE]
+ 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
+}
Modified: pkg/data/euler.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/ou2.rda
===================================================================
(Binary files differ)
Modified: pkg/data/rw2.rda
===================================================================
(Binary files differ)
Added: pkg/man/trajmatch.Rd
===================================================================
--- pkg/man/trajmatch.Rd (rev 0)
+++ pkg/man/trajmatch.Rd 2009-07-07 13:35:40 UTC (rev 155)
@@ -0,0 +1,48 @@
+\name{traj.match}
+\alias{traj.match}
+\title{Trajectory matching}
+\description{
+ Match trajectories.
+}
+\usage{
+traj.match(object, start, est, method = "Nelder-Mead", gr = NULL, \dots)
+}
+\arguments{
+ \item{object}{A \code{pomp} object.}
+ \item{start}{Guessed parameters.}
+ \item{est}{Character vector containing the names of parameters to be estimated.}
+ \item{method}{}
+ \item{gr}{}
+ \item{\dots}{Arguments that will be passed to \code{optim} in the
+ \code{control} list.}
+}
+\details{
+ Trajectory matching is accomplished using \code{\link{optim}}. It is
+ assumed that the process model is deterministic.
+}
+\examples{
+ data(ou2)
+ true.p <- c(
+ alpha.1=0.9,alpha.2=0,alpha.3=0,alpha.4=0.99,
+ sigma.1=1,sigma.2=0,sigma.3=2,
+ tau=1,
+ x1.0=50,x2.0=-50
+ )
+ simdata <- simulate(ou2,nsim=1,params=true.p,seed=43553)
+ guess.p <- true.p
+ guess.p[grep('sigma',names(guess.p))] <- 0
+ res <- traj.match(
+ simdata,
+ start=guess.p,
+ est=c('alpha.1','alpha.4','x1.0','x2.0','tau'),
+ maxit=2000,
+ reltol=1e-8
+ )
+}
+\seealso{
+ \code{\link{optim}},
+ \code{\link{pomp}},
+ \code{\link{pomp-class}}
+}
+\keyword{models}
+\keyword{ts}
Modified: pkg/tests/ou2-mif.Rout.save
===================================================================
--- pkg/tests/ou2-mif.Rout.save 2009-06-24 15:48:22 UTC (rev 154)
+++ pkg/tests/ou2-mif.Rout.save 2009-07-07 13:35:40 UTC (rev 155)
@@ -1,6 +1,6 @@
-R version 2.8.1 (2008-12-22)
-Copyright (C) 2008 The R Foundation for Statistical Computing
+R version 2.9.1 (2009-06-26)
+Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -70,8 +70,7 @@
> plot(mif1)
> compare.mif(mif2)
> try(compare.mif(mif1,mif2))
-Error in compare.mif(mif1, mif2) :
- unused argument(s) (<S4 object of class "mif">)
+Error in compare.mif(mif1, mif2) : unused argument(s) (mif2)
> compare.mif(list(mif1,mif2))
> dev.off()
null device
@@ -89,7 +88,7 @@
+ Np=100,cooling.factor=0.95,ic.lag=10,var.factor=1
+ )
+ )
-Error : mif error: ‘pars’ and ‘ivps’ must be mutually disjoint subsets of ‘names(start)’ and must have a positive random-walk SDs specified in ‘rw.sd’
+Error : mif error: 'pars' and 'ivps' must be mutually disjoint subsets of 'names(start)' and must have a positive random-walk SDs specified in 'rw.sd'
>
> try(
+ mif(
@@ -101,7 +100,7 @@
+ Np=100,cooling.factor=0.95,ic.lag=10,var.factor=1
+ )
+ )
-Error : mif error: ‘pars’ and ‘ivps’ must be mutually disjoint subsets of ‘names(start)’ and must have a positive random-walk SDs specified in ‘rw.sd’
+Error : mif error: 'pars' and 'ivps' must be mutually disjoint subsets of 'names(start)' and must have a positive random-walk SDs specified in 'rw.sd'
>
> try(
+ mif(
@@ -112,7 +111,7 @@
+ Np=100,cooling.factor=0.95,ic.lag=10,var.factor=1
+ )
+ )
-Error : mif error: ‘pars’ and ‘ivps’ must be mutually disjoint subsets of ‘names(start)’ and must have a positive random-walk SDs specified in ‘rw.sd’
+Error : mif error: 'pars' and 'ivps' must be mutually disjoint subsets of 'names(start)' and must have a positive random-walk SDs specified in 'rw.sd'
>
> try(
+ mif(
@@ -123,7 +122,7 @@
+ Np=100,ic.lag=10,var.factor=1
+ )
+ )
-Error : mif error: ‘cooling.factor’ must be specified
+Error : mif error: 'cooling.factor' must be specified
>
> try(
+ mif(
@@ -134,7 +133,7 @@
+ Np=-10,cooling.factor=0.95,ic.lag=10,var.factor=1
+ )
+ )
-Error : mif error: ‘Np’ must be a positive integer
+Error : mif error: 'Np' must be a positive integer
>
> try(
+ mif(
@@ -145,7 +144,7 @@
+ Np=11.6,cooling.factor=0.95,ic.lag=10,var.factor=1
+ )
+ )
-Error : mif error: ‘Nmif’ must be a positive integer
+Error : mif error: 'Nmif' must be a positive integer
>
> fit <- mif(
+ ou2,
@@ -177,7 +176,7 @@
> fit <- mif(fit,rw.sd=c(x1.0=5,x2.0=5,alpha.1=0.1,alpha.4=0.1))
> fit <- continue(fit,Nmif=2,ivps=c("x1.0"),pars=c("alpha.1"))
Warning message:
-mif warning: the variable(s) alpha.4, x2.0 have positive random-walk SDs specified, but are included in neither ‘pars’ nor ‘ivps’. These random walk SDs are ignored.
+mif warning: the variable(s) alpha.4, x2.0 have positive random-walk SDs specified, but are included in neither 'pars' nor 'ivps'. These random walk SDs are ignored.
> s <- coef(fit)
> s[2] <- 0.01
> fit <- mif(fit,Nmif=3,start=s)
@@ -185,5 +184,5 @@
> fit <- continue(fit,Nmif=2,Np=2000)
> fit <- continue(fit,ivps=c("x1.0"),rw.sd=c(alpha.1=0.1,alpha.4=0.1,x1.0=5,x2.0=5),Nmif=2)
Warning message:
-mif warning: the variable(s) x2.0 have positive random-walk SDs specified, but are included in neither ‘pars’ nor ‘ivps’. These random walk SDs are ignored.
+mif warning: the variable(s) x2.0 have positive random-walk SDs specified, but are included in neither 'pars' nor 'ivps'. These random walk SDs are ignored.
>
Added: pkg/tests/ou2-trajmatch.R
===================================================================
--- pkg/tests/ou2-trajmatch.R (rev 0)
+++ pkg/tests/ou2-trajmatch.R 2009-07-07 13:35:40 UTC (rev 155)
@@ -0,0 +1,36 @@
+library(pomp)
+
+data(ou2)
+true.p <- c(
+ alpha.1=0.9,alpha.2=0,alpha.3=0,alpha.4=0.99,
+ sigma.1=1,sigma.2=0,sigma.3=2,
+ tau=1,x1.0=50,x2.0=-50
+ )
+simdata <- simulate(ou2,nsim=5,params=true.p,seed=394885)
+guess.p <- true.p
+guess.p[grep('sigma',names(guess.p))] <- 0
+
+x <- sapply(
+ simdata,
+ function (d) {
+ res <- traj.match(
+ d,
+ start=guess.p,
+ est=c('alpha.1','alpha.4','x1.0','x2.0','tau'),
+ maxit=2000,
+ reltol=1e-8
+ )
+ c(conv=res$convergence,loglik=res$value,res$params)
+ }
+ )
+range(x['conv',])
+range(x['loglik',])
+print(
+ cbind(
+ truth=true.p,
+ mean.est=apply(x[names(true.p),],1,mean),
+ bias=apply(x[names(true.p),],1,mean)-true.p,
+ se=apply(x[names(true.p),],1,sd)
+ ),
+ digits=4
+ )
Added: pkg/tests/ou2-trajmatch.Rout.save
===================================================================
--- pkg/tests/ou2-trajmatch.Rout.save (rev 0)
+++ pkg/tests/ou2-trajmatch.Rout.save 2009-07-07 13:35:40 UTC (rev 155)
@@ -0,0 +1,70 @@
+
+R version 2.9.1 (2009-06-26)
+Copyright (C) 2009 The R Foundation for Statistical Computing
+ISBN 3-900051-07-0
+
+R is free software and comes with ABSOLUTELY NO WARRANTY.
+You are welcome to redistribute it under certain conditions.
+Type 'license()' or 'licence()' for distribution details.
+
+R is a collaborative project with many contributors.
+Type 'contributors()' for more information and
+'citation()' on how to cite R or R packages in publications.
+
+Type 'demo()' for some demos, 'help()' for on-line help, or
+'help.start()' for an HTML browser interface to help.
+Type 'q()' to quit R.
+
+> library(pomp)
+Loading required package: deSolve
+Loading required package: subplex
+Loading required package: mvtnorm
+>
+> data(ou2)
+> true.p <- c(
++ alpha.1=0.9,alpha.2=0,alpha.3=0,alpha.4=0.99,
++ sigma.1=1,sigma.2=0,sigma.3=2,
++ tau=1,x1.0=50,x2.0=-50
++ )
+> simdata <- simulate(ou2,nsim=5,params=true.p,seed=394885)
+> guess.p <- true.p
+> guess.p[grep('sigma',names(guess.p))] <- 0
+>
+> x <- sapply(
++ simdata,
++ function (d) {
++ res <- traj.match(
++ d,
++ start=guess.p,
++ est=c('alpha.1','alpha.4','x1.0','x2.0','tau'),
++ maxit=2000,
++ reltol=1e-8
++ )
++ c(conv=res$convergence,loglik=res$value,res$params)
++ }
++ )
+> range(x['conv',])
+[1] 0 0
+> range(x['loglik',])
+[1] -590.7141 -507.0953
+> print(
++ cbind(
++ truth=true.p,
++ mean.est=apply(x[names(true.p),],1,mean),
++ bias=apply(x[names(true.p),],1,mean)-true.p,
++ se=apply(x[names(true.p),],1,sd)
++ ),
++ digits=4
++ )
+ truth mean.est bias se
+alpha.1 0.90 0.9110 0.010953 0.022026
+alpha.2 0.00 0.0000 0.000000 0.000000
+alpha.3 0.00 0.0000 0.000000 0.000000
+alpha.4 0.99 0.9935 0.003476 0.002373
+sigma.1 1.00 0.0000 -1.000000 0.000000
+sigma.2 0.00 0.0000 0.000000 0.000000
+sigma.3 2.00 0.0000 -2.000000 0.000000
+tau 1.00 3.6906 2.690637 0.724948
+x1.0 50.00 47.2476 -2.752419 5.299588
+x2.0 -50.00 -48.8198 1.180172 5.418218
+>
More information about the pomp-commits
mailing list