[Pomp-commits] r789 - in branches/mif2: . R inst/doc man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jan 7 17:07:24 CET 2013
Author: kingaa
Date: 2013-01-07 17:07:24 +0100 (Mon, 07 Jan 2013)
New Revision: 789
Added:
branches/mif2/R/authors.R
Modified:
branches/mif2/DESCRIPTION
branches/mif2/R/aaa.R
branches/mif2/R/builder.R
branches/mif2/R/mif-class.R
branches/mif2/R/mif.R
branches/mif2/R/pfilter.R
branches/mif2/R/pomp-methods.R
branches/mif2/R/pomp.R
branches/mif2/R/traj-match.R
branches/mif2/inst/doc/advanced_topics_in_pomp.pdf
branches/mif2/inst/doc/intro_to_pomp.Rnw
branches/mif2/inst/doc/intro_to_pomp.pdf
branches/mif2/man/pomp.Rd
branches/mif2/man/traj-match.Rd
branches/mif2/src/dmeasure.c
branches/mif2/src/pomp.h
branches/mif2/src/pomp_fun.c
branches/mif2/src/pomp_internal.h
Log:
- merge in changes to main branch of pomp
Modified: branches/mif2/DESCRIPTION
===================================================================
--- branches/mif2/DESCRIPTION 2012-10-17 16:56:16 UTC (rev 788)
+++ branches/mif2/DESCRIPTION 2013-01-07 16:07:24 UTC (rev 789)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.42-6
-Date: 2012-06-01
+Version: 1.0-1
+Date: 2012-10-09
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>
URL: http://pomp.r-forge.r-project.org
Modified: branches/mif2/R/aaa.R
===================================================================
--- branches/mif2/R/aaa.R 2012-10-17 16:56:16 UTC (rev 788)
+++ branches/mif2/R/aaa.R 2013-01-07 16:07:24 UTC (rev 789)
@@ -17,6 +17,8 @@
setGeneric("pred.mean",function(object,...)standardGeneric("pred.mean"))
setGeneric("pred.var",function(object,...)standardGeneric("pred.var"))
setGeneric("filter.mean",function(object,...)standardGeneric("filter.mean"))
+setGeneric("cond.logLik",function(object,...)standardGeneric("cond.logLik"))
+setGeneric("eff.sample.size",function(object,...)standardGeneric("eff.sample.size"))
if (!exists("paste0",where="package:base")) {
paste0 <- function(...) paste(...,sep="")
Added: branches/mif2/R/authors.R
===================================================================
--- branches/mif2/R/authors.R (rev 0)
+++ branches/mif2/R/authors.R 2013-01-07 16:07:24 UTC (rev 789)
@@ -0,0 +1,12 @@
+list(
+ aak=person(given=c("Aaron","A."),family="King",role=c("aut","cre"),email="kingaa at umich.edu"),
+ eli=person(given=c("Edward","L."),family="Ionides",role=c("aut"),email="ionides at umich.edu"),
+ cb=person(given=c("Carles"),family="Breto",role=c("aut")),
+ spe=person(given=c("Stephen","P."),family="Ellner",role=c("aut")),
+ bek=person(given=c("Bruce","E."),family="Kendall",role=c("aut")),
+ mf=person(given=c("Matthew","J."),family="Ferrari",role=c("ctb")),
+ ml=person(given=c("Michael"),family="Lavine",role=c("ctb")),
+ dcr=person(given=c("Daniel","C."),family="Reuman",role=c("ctb")),
+ hw=person(given=c("Helen"),family="Wearing",role=c("ctb")),
+ snw=person(given=c("Simon","N."),family="Wood",role=c("ctb"))
+ ) -> author.list
Modified: branches/mif2/R/builder.R
===================================================================
--- branches/mif2/R/builder.R 2012-10-17 16:56:16 UTC (rev 788)
+++ branches/mif2/R/builder.R 2013-01-07 16:07:24 UTC (rev 789)
@@ -110,9 +110,7 @@
reulermultinom="\tvoid (*reulermultinom)(int,double,double*,double,double*);\nreulermultinom = (void (*)(int,double,double*,double,double*)) R_GetCCallable(\"pomp\",\"reulermultinom\");\n",
get_pomp_userdata="\tconst SEXP (*get_pomp_userdata)(const char *);\npomp_get_userdata = (const SEXP (*)(const char*)) R_GetCCallable(\"pomp\",\"get_pomp_userdata\");\n",
get_pomp_userdata_int="\tconst int * (*get_pomp_userdata_int)(const char *);\npomp_get_userdata_int = (const int *(*)(const char*)) R_GetCCallable(\"pomp\",\"get_pomp_userdata_int\");\n",
- get_pomp_userdata_double="\tconst double * (*get_pomp_userdata_double)(const char *);\npomp_get_userdata_double = (const double *(*)(const char*)) R_GetCCallable(\"pomp\",\"get_pomp_userdata_double\");\n",
- logit="\tdouble logit(double p){return log(p/(1.0-p));}\n\n",
- expit="\tdouble expit(double x){return 1.0/(1.0+exp(-x));}\n\n"
+ get_pomp_userdata_double="\tconst double * (*get_pomp_userdata_double)(const char *);\npomp_get_userdata_double = (const double *(*)(const char*)) R_GetCCallable(\"pomp\",\"get_pomp_userdata_double\");\n"
)
footer <- list(
Modified: branches/mif2/R/mif-class.R
===================================================================
--- branches/mif2/R/mif-class.R 2012-10-17 16:56:16 UTC (rev 788)
+++ branches/mif2/R/mif-class.R 2013-01-07 16:07:24 UTC (rev 789)
@@ -6,10 +6,10 @@
transform = "logical",
ivps = 'character',
pars = 'character',
- Nmif = 'integer',
- option='character',
- cooling.fraction='numeric',
- particles = 'function',
+ Nmif = 'integer',
+ option='character',
+ cooling.fraction='numeric',
+ particles = 'function',
var.factor='numeric',
ic.lag='integer',
cooling.factor='numeric',
Modified: branches/mif2/R/mif.R
===================================================================
--- branches/mif2/R/mif.R 2012-10-17 16:56:16 UTC (rev 788)
+++ branches/mif2/R/mif.R 2013-01-07 16:07:24 UTC (rev 789)
@@ -21,8 +21,8 @@
list(alpha=alpha,gamma=alpha^2)
}
-mif.cooling2 <- function (cooling.fraction, nt, m, n) { # cooling schedule for mif2
- ## cooling.fraction now is the fraction of cooling on 50 iteration
+mif2.cooling <- function (cooling.fraction, nt, m, n) { # cooling schedule for mif2
+ ## cooling.fraction is the fraction of cooling after 50 iterations
cooling.scalar <- (50*n*cooling.fraction-1)/(1-cooling.fraction)
alpha <- (1+cooling.scalar)/(cooling.scalar+nt+n*(m-1))
list(alpha=alpha,gamma=alpha^2)
@@ -46,7 +46,7 @@
rw.sd,
option, cooling.fraction, paramMatrix,
Np, cooling.factor, var.factor, ic.lag,
- method, tol, max.fail,
+ tol, max.fail,
verbose, transform, .ndone) {
transform <- as.logical(transform)
@@ -163,19 +163,11 @@
)
}
- if (missing(option) && missing(method))
- { option <- 'mif'
- warning(sQuote("mif")," warning: use mif as default")
- }
- if (missing(option) && !missing(method) )
- { option <- method
- warning(sQuote("mif")," warning: ",sQuote("method")," flag is deprecated, use ",sQuote("option"))
- }
if (missing(cooling.factor)&&(option=="mif2")) ##Default value for the slot cooling.fraction
- cooling.factor=1
+ cooling.factor <- 1
if (missing(cooling.factor)&&(option!="mif2"))
stop("mif error: ",sQuote("cooling.factor")," must be specified",call.=FALSE)
- if ((option!="mif2")&&((length(cooling.factor)!=1)||(cooling.factor < 0)||(cooling.factor>1)))
+ if ((option!="mif2")&&((length(cooling.factor)!=1)||(cooling.factor<0)||(cooling.factor>1)))
stop("mif error: ",sQuote("cooling.factor")," must be a number between 0 and 1",call.=FALSE)
if (missing(var.factor))
@@ -189,11 +181,11 @@
if (Nmif<0)
stop("mif error: ",sQuote("Nmif")," must be a positive integer",call.=FALSE)
if (option=="mif2" && missing(cooling.fraction))
- stop("mif error: ",sQuote("cooling.fraction")," must be specified for method mif2",call.=FALSE)
+ stop("mif error: ",sQuote("cooling.fraction")," must be specified for method ",sQuote("mif2"),call.=FALSE)
if (option=="mif2")
cooling.fraction <- as.numeric(cooling.fraction)
if (missing(cooling.fraction)&&(option!="mif2")) ##Default value for the slot cooling.fraction
- cooling.fraction=-1
+ cooling.fraction <- NA
theta <- start
@@ -237,84 +229,101 @@
for (n in seq_len(Nmif)) {
##cooling schedule
- switch(option, mif2={
- cool.sched <- try(mif.cooling2(cooling.fraction, 1 ,
- .ndone +n, ntimes), silent = FALSE)
- },
- { #not mif2
- cool.sched <- try(mif.cooling(cooling.factor, .ndone +
- n), silent = FALSE)
- } )
+ cool.sched <- try(
+ switch(
+ option,
+ mif2=mif2.cooling(cooling.fraction,1,.ndone+n,ntimes),
+ mif.cooling(cooling.factor,.ndone+n)
+ ),
+ silent=FALSE
+ )
if (inherits(cool.sched, "try-error"))
- stop("mif error: cooling schedule error", call. = FALSE)
- sigma.n <- sigma * cool.sched$alpha
+ stop("mif error: cooling schedule error",call.=FALSE)
+ sigma.n <- sigma*cool.sched$alpha
- P <- try(particles(tmp.mif, Np = Np[1], center = theta,
- sd = sigma.n * var.factor), silent = FALSE)
+ P <- try(
+ particles(
+ tmp.mif,
+ Np=Np[1],
+ center=theta,
+ sd=sigma.n*var.factor
+ ),
+ silent = FALSE
+ )
if (inherits(P, "try-error"))
- stop("mif error: error in ", sQuote("particles"),
- call. = FALSE)
- ## Setting up parameter switch
- switch(option, mif2={
- if(!((n==1)&&(missing(paramMatrix)))) #use paramMatrix if it exists
- {
- P[pars, ] <- paramMatrix[[1]][pars,]
- }
+ stop("mif error: error in ",sQuote("particles"),call.=FALSE)
+
+ if (option=="mif2") {
+ if(!((n==1)&&(missing(paramMatrix)))) { #use paramMatrix if it exists
+ P[pars,] <- paramMatrix[[1]][pars,]
+ }
cooling.m=n
- },
- {
-
- }
- )
- pfp <- try(pfilter.internal(object = obj, params = P,
- tol = tol, max.fail = max.fail, pred.mean = (n ==
- Nmif), pred.var = TRUE, filter.mean = TRUE, save.states = FALSE,
- save.params = FALSE, .rw.sd = sigma.n[pars],cooling.m=cooling.m, cooling.fraction=cooling.fraction, paramMatrix=paramMatrix,option=option,
- verbose = verbose, transform = transform, .ndone=.ndone), silent = FALSE)
+ }
+
+ pfp <- try(
+ pfilter.internal(
+ object=obj,
+ params=P,
+ tol=tol,
+ max.fail=max.fail,
+ pred.mean=(n==Nmif),
+ pred.var=TRUE,
+ filter.mean=TRUE,
+ save.states=FALSE,
+ save.params=FALSE,
+ .rw.sd=sigma.n[pars],
+ cooling.m=cooling.m,
+ cooling.fraction=cooling.fraction,
+ paramMatrix=paramMatrix,
+ option=option,
+ verbose=verbose,
+ transform=transform,
+ .ndone=.ndone
+ ),
+ silent=FALSE
+ )
if (inherits(pfp, "try-error"))
- stop("mif error: error in ", sQuote("pfilter"),
- call. = FALSE)
+ stop("mif error: error in ",sQuote("pfilter"),call.=FALSE)
##update switch
- switch(option, mif = {
- v <- pfp$pred.var[pars, , drop = FALSE]
- v1 <- cool.sched$gamma * (1 + var.factor^2) *
- sigma[pars]^2
- theta.hat <- cbind(theta[pars], pfp$filter.mean[pars,
- , drop = FALSE])
- theta[pars] <- theta[pars] + colSums(apply(theta.hat,
- 1, diff)/t(v)) * v1
- }, unweighted = {
- theta.hat <- pfp$filter.mean[pars,,drop=FALSE]
- theta[pars] <- rowMeans(theta.hat)
- }, fp = {
- theta.hat <- pfp$filter.mean[pars, ntimes, drop = FALSE]
- theta[pars] <- theta.hat
- }, mif2 ={
- paramMatrix <- pfp$paramMatrix
- theta.hat <- rowMeans(pfp$paramMatrix[[1]][pars, , drop = FALSE])
- theta[pars] <- theta.hat
-
- },)
- theta[ivps] <- pfp$filter.mean[ivps, ic.lag]
- conv.rec[n + 1, -c(1, 2)] <- theta
- conv.rec[n, c(1, 2)] <- c(pfp$loglik, pfp$nfail)
+ switch(
+ option,
+ mif = {
+ v <- pfp$pred.var[pars, , drop = FALSE]
+ v1 <- cool.sched$gamma * (1 + var.factor^2) * sigma[pars]^2
+ theta.hat <- cbind(theta[pars], pfp$filter.mean[pars,,drop = FALSE])
+ theta[pars] <- theta[pars] + colSums(apply(theta.hat,1,diff)/t(v))*v1
+ },
+ unweighted = {
+ theta.hat <- pfp$filter.mean[pars,,drop=FALSE]
+ theta[pars] <- rowMeans(theta.hat)
+ }, fp = {
+ theta.hat <- pfp$filter.mean[pars, ntimes, drop = FALSE]
+ theta[pars] <- theta.hat
+ }, mif2 ={
+ paramMatrix <- pfp$paramMatrix
+ theta.hat <- rowMeans(pfp$paramMatrix[[1]][pars, , drop = FALSE])
+ theta[pars] <- theta.hat
+ },
+ stop("unrecognized option: ",sQuote(option))
+ )
+ theta[ivps] <- pfp$filter.mean[ivps,ic.lag]
+ conv.rec[n+1,-c(1,2)] <- theta
+ conv.rec[n,c(1,2)] <- c(pfp$loglik,pfp$nfail)
- if (verbose)
- cat("MIF iteration ", n, " of ", Nmif, " completed\n")
+ if (verbose) cat("MIF iteration ",n," of ",Nmif," completed\n")
}
## back transform the parameter estimate if necessary
if (transform)
theta <- partrans(pfp,theta,dir="forward")
- if(option=='mif2' && missing(paramMatrix))
- { paramMatrix <-vector(mode="list", length=1)
- paramMatrix[[1]]<-pfp$paramMatrix[[1]]
-
- }
- if(option!='mif2' && missing(paramMatrix))
- { paramMatrix <-list()
- }
+ if(option=='mif2' && missing(paramMatrix)) {
+ paramMatrix <-vector(mode="list", length=1)
+ paramMatrix[[1]]<-pfp$paramMatrix[[1]]
+ }
+ if(option!='mif2' && missing(paramMatrix)) {
+ paramMatrix <-list()
+ }
new(
"mif",
pfp,
@@ -332,7 +341,7 @@
conv.rec=conv.rec,
option=option,
paramMatrix=paramMatrix,
- cooling.fraction = cooling.fraction
+ cooling.fraction=cooling.fraction
)
}
@@ -430,7 +439,6 @@
var.factor=var.factor,
ic.lag=ic.lag,
option=option,
- method=method,
paramMatrix=paramMatrix,
cooling.fraction = cooling.fraction,
tol=tol,
@@ -526,7 +534,6 @@
var.factor=var.factor,
ic.lag=ic.lag,
option=option,
- method=method,
cooling.fraction=cooling.fraction,
paramMatrix=paramMatrix,
tol=tol,
@@ -685,7 +692,6 @@
ic.lag=ic.lag,
option=option,
cooling.fraction=cooling.fraction,
- method=method,
paramMatrix=paramMatrix,
tol=tol,
max.fail=max.fail,
Modified: branches/mif2/R/pfilter.R
===================================================================
--- branches/mif2/R/pfilter.R 2012-10-17 16:56:16 UTC (rev 788)
+++ branches/mif2/R/pfilter.R 2013-01-07 16:07:24 UTC (rev 789)
@@ -29,8 +29,9 @@
cooling.fraction, cooling.m, option,
.rw.sd, seed, verbose,
save.states, save.params,
- transform,.ndone) {
-
+ transform,
+ .ndone) {
+
if (missing(option)) option<-'mif'
if(missing(paramMatrix)) paramMatrix <- array(dim=c(0,0))
transform <- as.logical(transform)
@@ -183,7 +184,7 @@
if (option=="mif2") {
cool.sched <- try(
- mif.cooling2(cooling.fraction,nt,cooling.m+.ndone,ntimes),
+ mif2.cooling(cooling.fraction,nt,cooling.m+.ndone,ntimes),
silent=FALSE
)
if (inherits(cool.sched,"try-error"))
Modified: branches/mif2/R/pomp-methods.R
===================================================================
--- branches/mif2/R/pomp-methods.R 2012-10-17 16:56:16 UTC (rev 788)
+++ branches/mif2/R/pomp-methods.R 2013-01-07 16:07:24 UTC (rev 789)
@@ -259,6 +259,7 @@
object at params <- params
}
}
+ storage.mode(object at params) <- "double"
object
}
)
Modified: branches/mif2/R/pomp.R
===================================================================
--- branches/mif2/R/pomp.R 2012-10-17 16:56:16 UTC (rev 788)
+++ branches/mif2/R/pomp.R 2013-01-07 16:07:24 UTC (rev 789)
@@ -50,13 +50,20 @@
rmeasure, dmeasure, measurement.model,
skeleton = NULL, skeleton.type = c("map","vectorfield"),
skelmap.delta.t = 1,
- initializer, covar, tcovar,
+ initializer, params, covar, tcovar,
obsnames, statenames, paramnames, covarnames, zeronames,
PACKAGE, parameter.transform, parameter.inv.transform) {
if (missing(data)) stop(sQuote("data")," is a required argument")
if (missing(times)) stop(sQuote("times")," is a required argument")
if (missing(t0)) stop(sQuote("t0")," is a required argument")
+ if (missing(params)) params <- numeric(0)
+
+ if (length(params)>0) {
+ if (is.null(names(params)) || !is.numeric(params))
+ stop("pomp error: ",sQuote("params")," must be a named numeric vector",call.=TRUE)
+ }
+ storage.mode(params) <- 'double'
## check the data
if (is.data.frame(data)) {
@@ -273,6 +280,7 @@
times = times,
t0 = t0,
initializer = initializer,
+ params=params,
covar = covar,
tcovar = tcovar,
obsnames = obsnames,
@@ -374,7 +382,7 @@
rmeasure, dmeasure, measurement.model,
skeleton = NULL, skeleton.type = c("map","vectorfield"),
skelmap.delta.t = 1,
- initializer, covar, tcovar,
+ initializer, params, covar, tcovar,
obsnames, statenames, paramnames, covarnames, zeronames,
PACKAGE, parameter.transform, parameter.inv.transform) {
pomp.constructor(
@@ -390,6 +398,7 @@
skeleton.type=skeleton.type,
skelmap.delta.t=skelmap.delta.t,
initializer=initializer,
+ params=params,
covar=covar,
tcovar=tcovar,
obsnames=obsnames,
@@ -412,7 +421,7 @@
rmeasure, dmeasure, measurement.model,
skeleton = NULL, skeleton.type = c("map","vectorfield"),
skelmap.delta.t = 1,
- initializer, covar, tcovar,
+ initializer, params, covar, tcovar,
obsnames, statenames, paramnames, covarnames, zeronames,
PACKAGE, parameter.transform, parameter.inv.transform) {
pomp.constructor(
@@ -428,6 +437,7 @@
skeleton.type=skeleton.type,
skelmap.delta.t=skelmap.delta.t,
initializer=initializer,
+ params=params,
covar=covar,
tcovar=tcovar,
obsnames=obsnames,
@@ -451,7 +461,7 @@
rmeasure, dmeasure, measurement.model,
skeleton = NULL, skeleton.type = c("map","vectorfield"),
skelmap.delta.t = 1,
- initializer, covar, tcovar,
+ initializer, params, covar, tcovar,
obsnames, statenames, paramnames, covarnames, zeronames,
PACKAGE, parameter.transform, parameter.inv.transform) {
pomp.constructor(
@@ -467,6 +477,7 @@
skeleton.type=skeleton.type,
skelmap.delta.t=skelmap.delta.t,
initializer=initializer,
+ params=params,
covar=covar,
tcovar=tcovar,
obsnames=obsnames,
@@ -488,7 +499,7 @@
function (data, times, t0, ..., rprocess, dprocess,
rmeasure, dmeasure, measurement.model,
skeleton, skeleton.type, skelmap.delta.t,
- initializer, covar, tcovar,
+ initializer, params, covar, tcovar,
obsnames, statenames, paramnames, covarnames, zeronames,
PACKAGE, parameter.transform, parameter.inv.transform) {
mmg <- !missing(measurement.model)
@@ -512,6 +523,7 @@
if (missing(rprocess)) rprocess <- data at rprocess
if (missing(dprocess)) dprocess <- data at dprocess
if (missing(initializer)) initializer <- data at initializer
+ if (missing(params)) params <- data at params
if (missing(covar)) covar <- data at covar
if (missing(tcovar)) tcovar <- data at tcovar
if (missing(obsnames)) obsnames <- data at obsnames
@@ -524,8 +536,6 @@
if (missing(skeleton)) skeleton <- data at skeleton
if (missing(skelmap.delta.t)) skelmap.delta.t <- data at skelmap.delta.t
- pars <- coef(data)
-
if (missing(parameter.transform)) {
if (missing(parameter.inv.transform)) {
par.trans <- data at par.trans
@@ -577,7 +587,7 @@
)
) -> retval
- coef(retval) <- pars
+ coef(retval) <- params
retval
}
Modified: branches/mif2/R/traj-match.R
===================================================================
--- branches/mif2/R/traj-match.R 2012-10-17 16:56:16 UTC (rev 788)
+++ branches/mif2/R/traj-match.R 2013-01-07 16:07:24 UTC (rev 789)
@@ -32,7 +32,7 @@
}
)
-traj.match.objfun <- function (object, params, est, transform = FALSE) {
+traj.match.objfun <- function (object, params, est, transform = FALSE, ...) {
transform <- as.logical(transform)
if (missing(est)) est <- character(0)
@@ -53,7 +53,7 @@
d <- dmeasure(
object,
y=object at data,
- x=trajectory(object,params=tparams),
+ x=trajectory(object,params=tparams,...),
times=time(object),
params=tparams,
log=TRUE
@@ -62,7 +62,7 @@
d <- dmeasure(
object,
y=object at data,
- x=trajectory(object,params=params),
+ x=trajectory(object,params=params,...),
times=time(object),
params=params,
log=TRUE
Modified: branches/mif2/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: branches/mif2/inst/doc/intro_to_pomp.Rnw
===================================================================
--- branches/mif2/inst/doc/intro_to_pomp.Rnw 2012-10-17 16:56:16 UTC (rev 788)
+++ branches/mif2/inst/doc/intro_to_pomp.Rnw 2013-01-07 16:07:24 UTC (rev 789)
@@ -997,22 +997,22 @@
<<>>=
pf.truth <- pfilter(ricker,Np=1000,max.fail=50,seed=1066L)
pf.guess <- pfilter(ricker,params=guess,Np=1000,max.fail=50,seed=1066L)
-#pf.mf <- pfilter(mf,Np=1000,seed=1066L)
-#pf.pm <- pfilter(pm,Np=1000,max.fail=10,seed=1066L)
+pf.mf <- pfilter(mf,Np=1000,seed=1066L)
+pf.pm <- pfilter(pm,Np=1000,max.fail=10,seed=1066L)
pb.mf <- probe(mf,nsim=1000,probes=plist,seed=1066L)
res <- rbind(
cbind(guess=guess,truth=coef(ricker),MLE=coef(mf),PM=coef(pm)),
loglik=c(
pf.guess$loglik,
- pf.truth$loglik
- #pf.mf$loglik
- #pf.pm$loglik
+ pf.truth$loglik,
+ pf.mf$loglik,
+ pf.pm$loglik
),
synth.loglik=c(
summary(pb.guess)$synth.loglik,
summary(pb.truth)$synth.loglik,
- summary(pb.mf)$synth.loglik
- # summary(pm)$synth.loglik
+ summary(pb.mf)$synth.loglik,
+ summary(pm)$synth.loglik
)
)
Modified: branches/mif2/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Modified: branches/mif2/man/pomp.Rd
===================================================================
--- branches/mif2/man/pomp.Rd 2012-10-17 16:56:16 UTC (rev 788)
+++ branches/mif2/man/pomp.Rd 2013-01-07 16:07:24 UTC (rev 789)
@@ -16,25 +16,25 @@
\S4method{pomp}{data.frame}(data, times, t0, \dots, rprocess, dprocess, rmeasure, dmeasure,
measurement.model,
skeleton = NULL, skeleton.type = c("map","vectorfield"), skelmap.delta.t = 1,
- initializer, covar, tcovar,
+ initializer, params, covar, tcovar,
obsnames, statenames, paramnames, covarnames, zeronames,
PACKAGE, parameter.transform, parameter.inv.transform)
\S4method{pomp}{numeric}(data, times, t0, \dots, rprocess, dprocess, rmeasure, dmeasure,
measurement.model,
skeleton = NULL, skeleton.type = c("map","vectorfield"), skelmap.delta.t = 1,
- initializer, covar, tcovar,
+ initializer, params, covar, tcovar,
obsnames, statenames, paramnames, covarnames, zeronames,
PACKAGE, parameter.transform, parameter.inv.transform)
\S4method{pomp}{matrix}(data, times, t0, \dots, rprocess, dprocess, rmeasure, dmeasure,
measurement.model,
skeleton = NULL, skeleton.type = c("map","vectorfield"), skelmap.delta.t = 1,
- initializer, covar, tcovar,
+ initializer, params, covar, tcovar,
obsnames, statenames, paramnames, covarnames, zeronames,
PACKAGE, parameter.transform, parameter.inv.transform)
\S4method{pomp}{pomp}(data, times, t0, \dots, rprocess, dprocess, rmeasure, dmeasure,
measurement.model,
skeleton, skeleton.type, skelmap.delta.t,
- initializer, covar, tcovar,
+ initializer, params, covar, tcovar,
obsnames, statenames, paramnames, covarnames, zeronames,
PACKAGE, parameter.transform, parameter.inv.transform)
}
@@ -108,6 +108,9 @@
These are simply copied over as initial conditions when \code{init.state} is called (see \code{\link{init.state-pomp}}).
The names of the state variables are the same as the corresponding initial value parameters, but with the \dQuote{\code{.0}} dropped.
}
+ \item{params}{
+ optional named numeric vector of parameters.
+ }
\item{covar, tcovar}{
An optional table of covariates: \code{covar} is the table (with one column per variable) and \code{tcovar} the corresponding times (one entry per row of \code{covar}).
\code{covar} can be specified as either a matrix or a data frame.
@@ -158,7 +161,7 @@
}
In general, the specification of process-model codes \code{rprocess} and/or \code{dprocess} can be somewhat nontrivial:
for this reason, \code{\link{plugins}} have been developed to streamline this process for the user.
- Currently, if one's process model evolves in discrete time or one is willing to make such an approximation (e.g., via an Euler approximation), then the \code{\link{euler.sim}} or \code{\link{onestep.sim}} plugin for \code{rprocess} and \code{\link{onestep.dens}} plugin for \code{dprocess} are available.
+ Currently, if one's process model evolves in discrete time or one is willing to make such an approximation (e.g., via an Euler approximation), then the \code{\link{euler.sim}}, \code{\link{discrete.time.sim}}, or \code{\link{onestep.sim}} plugin for \code{rprocess} and \code{\link{onestep.dens}} plugin for \code{dprocess} are available.
For exact simulation of certain continuous-time Markov chains, an implementation of Gillespie's algorithm is available (see \code{\link{gillespie.sim}}).
To use the plugins, consult the help documentation (\code{?\link{plugins}}) and the vignettes.
Modified: branches/mif2/man/traj-match.Rd
===================================================================
--- branches/mif2/man/traj-match.Rd 2012-10-17 16:56:16 UTC (rev 788)
+++ branches/mif2/man/traj-match.Rd 2013-01-07 16:07:24 UTC (rev 789)
@@ -25,7 +25,7 @@
\S4method{traj.match}{traj.matched.pomp}(object, start, est,
method = c("Nelder-Mead","subplex","SANN","BFGS","sannbox"),
gr = NULL, eval.only = FALSE, transform, \dots)
-traj.match.objfun(object, params, est, transform = FALSE)
+traj.match.objfun(object, params, est, transform = FALSE, \dots)
}
\arguments{
\item{object}{
@@ -63,7 +63,10 @@
if \code{TRUE}, optimization is performed on the transformed scale.
}
\item{\dots}{
- Arguments that will be passed to \code{\link{optim}}, \code{\link[subplex]{subplex}} or \code{\link{sannbox}}, via their \code{control} lists.
+ Extra arguments that will be passed either to the optimizer (\code{\link{optim}}, \code{\link[subplex]{subplex}} or \code{\link{sannbox}}, via their \code{control} lists) or to the ODE integrator.
+ In \code{traj.match}, extra arguments will be passed to the optimizer.
+ In \code{traj.match.objfun}, extra arguments are passed to the ODE integrator, but only if an ODE integrator is called, i.e., only when the deterministic skeleton is a vectorfield.
+ If extra arguments are needed by both optimizer and ODE integrator, construct an objective function first using \code{traj.match.objfun}, then give this objective function to the optimizer.
}
}
\details{
@@ -137,6 +140,10 @@
data(ricker)
ofun <- traj.match.objfun(ricker,est=c("r","phi"),transform=TRUE)
optim(fn=ofun,par=c(2,0),method="BFGS")
+
+ data(bbs)
+ ofun <- traj.match.objfun(bbs,est=c("beta","gamma"),transform=TRUE,hmax=0.001)
+ optim(fn=ofun,par=c(0,-1),method="Nelder-Mead")
}
\seealso{
\code{\link{trajectory}},
Modified: branches/mif2/src/dmeasure.c
===================================================================
--- branches/mif2/src/dmeasure.c 2012-10-17 16:56:16 UTC (rev 788)
+++ branches/mif2/src/dmeasure.c 2013-01-07 16:07:24 UTC (rev 789)
@@ -56,7 +56,7 @@
PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params))); nprotect++;
PROTECT(Cnames = GET_COLNAMES(GET_DIMNAMES(GET_SLOT(object,install("covar"))))); nprotect++;
- give_log = *(INTEGER(log));
+ give_log = *(INTEGER(AS_INTEGER(log)));
// set up the covariate table
covariate_table = make_covariate_table(object,&ncovars);
Modified: branches/mif2/src/pomp.h
===================================================================
--- branches/mif2/src/pomp.h 2012-10-17 16:56:16 UTC (rev 788)
+++ branches/mif2/src/pomp.h 2013-01-07 16:07:24 UTC (rev 789)
@@ -6,6 +6,7 @@
#include <R.h>
#include <Rmath.h>
#include <Rdefines.h>
+#include <R_ext/Rdynload.h>
// facility for extracting R objects from the 'userdata' slot
const SEXP get_pomp_userdata (const char *name);
@@ -192,7 +193,14 @@
// a vector of parameters ('coef') against a vector of basis-function values ('basis')
double dot_product (int dim, const double *basis, const double *coef);
+static R_INLINE double logit (double p) {
+ return log(p/(1.0-p));
+}
+static R_INLINE double expit (double x) {
+ return 1.0/(1.0+exp(-x));
+}
+
// prototypes for C-level access to Euler-multinomial distribution functions
// simulate Euler-multinomial transitions
Modified: branches/mif2/src/pomp_fun.c
===================================================================
--- branches/mif2/src/pomp_fun.c 2012-10-17 16:56:16 UTC (rev 788)
+++ branches/mif2/src/pomp_fun.c 2013-01-07 16:07:24 UTC (rev 789)
@@ -33,7 +33,7 @@
PROTECT(f = getListElement(nsi,"address")); nprotect++;
break;
default:
- error("'pomp_fun_handler': invalid 'mode' value");
+ error("operation cannot be completed: some needed function has not been specified");
break;
}
Modified: branches/mif2/src/pomp_internal.h
===================================================================
--- branches/mif2/src/pomp_internal.h 2012-10-17 16:56:16 UTC (rev 788)
+++ branches/mif2/src/pomp_internal.h 2013-01-07 16:07:24 UTC (rev 789)
@@ -240,14 +240,6 @@
return x;
}
-static R_INLINE double expit (double x) {
- return 1.0/(1.0 + exp(-x));
-}
-
-static R_INLINE double logit (double x) {
- return log(x/(1-x));
-}
-
static R_INLINE SEXP getListElement (SEXP list, const char *str)
{
SEXP elmt = R_NilValue;
More information about the pomp-commits
mailing list