[Pomp-commits] r629 - 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 Mar 27 13:57:42 CEST 2012
Author: kingaa
Date: 2012-03-27 13:57:41 +0200 (Tue, 27 Mar 2012)
New Revision: 629
Modified:
pkg/DESCRIPTION
pkg/R/bsmc.R
pkg/R/mif-class.R
pkg/R/mif-methods.R
pkg/R/mif.R
pkg/R/nlf-objfun.R
pkg/R/nlf.R
pkg/R/pfilter.R
pkg/R/pmcmc.R
pkg/R/pomp.R
pkg/R/probe-match.R
pkg/R/traj-match.R
pkg/data/blowflies.rda
pkg/data/dacca.rda
pkg/data/euler.sir.rda
pkg/data/gillespie.sir.rda
pkg/data/gompertz.rda
pkg/data/ou2.rda
pkg/data/ricker.rda
pkg/data/rw2.rda
pkg/data/verhulst.rda
pkg/inst/TODO
pkg/inst/data-R/dacca.R
pkg/inst/data-R/gompertz.R
pkg/inst/data-R/ricker.R
pkg/inst/data-R/sir.R
pkg/inst/doc/advanced_topics_in_pomp.Rnw
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/gompertz-multi-mif.rda
pkg/inst/doc/gompertz-trajmatch.rda
pkg/inst/doc/intro_to_pomp.Rnw
pkg/inst/doc/intro_to_pomp.pdf
pkg/inst/doc/manual.pdf
pkg/inst/doc/nlf-block-boot.rda
pkg/inst/doc/nlf-boot.rda
pkg/inst/doc/nlf-fit-from-truth.rda
pkg/inst/doc/nlf-fits.rda
pkg/inst/doc/nlf-lag-tests.rda
pkg/inst/doc/nlf-multi-short.rda
pkg/inst/doc/ricker-mif.rda
pkg/inst/doc/ricker-probe-match.rda
pkg/man/bsmc.Rd
pkg/man/dacca.Rd
pkg/man/mif-class.Rd
pkg/man/mif.Rd
pkg/man/nlf.Rd
pkg/man/pmcmc.Rd
pkg/man/probe.Rd
pkg/man/sir.Rd
pkg/man/traj-match.Rd
pkg/src/cholmodel.c
pkg/src/gompertz.c
pkg/src/ricker.c
pkg/src/sir.c
pkg/tests/dacca.R
pkg/tests/dacca.Rout.save
pkg/tests/dimchecks.R
pkg/tests/dimchecks.Rout.save
pkg/tests/gillespie.R
pkg/tests/gillespie.Rout.save
pkg/tests/gompertz.R
pkg/tests/gompertz.Rout.save
pkg/tests/ou2-mif.R
pkg/tests/ou2-mif.Rout.save
pkg/tests/ou2-pmcmc.R
pkg/tests/ou2-trajmatch.Rout.save
pkg/tests/partrans.R
pkg/tests/partrans.Rout.save
pkg/tests/pfilter.R
pkg/tests/pfilter.Rout.save
pkg/tests/pomppomp.R
pkg/tests/pomppomp.Rout.save
pkg/tests/ricker-probe.R
pkg/tests/ricker-probe.Rout.save
pkg/tests/ricker-spect.R
pkg/tests/ricker-spect.Rout.save
pkg/tests/ricker.R
pkg/tests/ricker.Rout.save
pkg/tests/skeleton.R
pkg/tests/skeleton.Rout.save
pkg/tests/steps.R
pkg/tests/steps.Rout.save
Log:
- put parameter transformation options into 'mif', 'bsmc', 'pmcmc', 'probe-match', 'traj-match', and 'nlf'
- reverse direction of parameter transformations in data()-loadable examples
- update documentation
- make 'bsmc' into an S4 method
- plot method for 'bsmcd.pomp' objects
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2012-03-21 19:52:13 UTC (rev 628)
+++ pkg/DESCRIPTION 2012-03-27 11:57:41 UTC (rev 629)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.40-10
-Date: 2012-03-21
+Version: 0.41-1
+Date: 2012-03-27
Revision: $Rev$
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/R/bsmc.R
===================================================================
--- pkg/R/bsmc.R 2012-03-21 19:52:13 UTC (rev 628)
+++ pkg/R/bsmc.R 2012-03-27 11:57:41 UTC (rev 629)
@@ -11,6 +11,24 @@
## lower = lower bounds on prior
## upper = upper bounds on prior
+setClass(
+ "bsmcd.pomp",
+ contains="pomp",
+ representation=representation(
+ transform="logical",
+ post="array",
+ prior="array",
+ est="character",
+ eff.sample.size="numeric",
+ smooth="numeric",
+ seed="integer",
+ nfail="integer",
+ cond.log.evidence="numeric",
+ log.evidence="numeric",
+ weights="numeric"
+ )
+ )
+
setGeneric("bsmc",function(object,...)standardGeneric("bsmc"))
setMethod(
@@ -24,8 +42,11 @@
seed = NULL,
verbose = getOption("verbose"),
max.fail = 0,
+ transform = FALSE,
...) {
+ transform <- as.logical(transform)
+
if (missing(seed)) seed <- NULL
if (!is.null(seed)) {
if (!exists(".Random.seed",where=.GlobalEnv))
@@ -37,8 +58,8 @@
error.prefix <- paste(sQuote("bsmc"),"error: ")
if (missing(params)) {
- if (length(object at params)>0) {
- params <- object at params
+ if (length(coef(object))>0) {
+ params <- coef(object)
} else {
stop(error.prefix,sQuote("params")," must be supplied",call.=FALSE)
}
@@ -46,8 +67,11 @@
if (missing(Np)) Np <- NCOL(params)
else if (is.matrix(params)&&(Np!=ncol(params)))
- stop(error.prefix,sQuote("Np")," cannot be other than ",sQuote("ncol(params)"),call.=FALSE)
+ warning(sQuote("Np")," is ignored when ",sQuote("params")," is a matrix")
+ if (transform)
+ params <- partrans(object,params,dir="inverse")
+
ntimes <- length(time(object))
if (is.null(dim(params))) {
params <- matrix(
@@ -60,10 +84,10 @@
)
)
}
- prior <- params
npars <- nrow(params)
paramnames <- rownames(params)
+ prior <- params
if (missing(est))
est <- paramnames[apply(params,1,function(x)diff(range(x))>0)]
@@ -108,7 +132,14 @@
}
}
- xstart <- init.state(object,params=params)
+ xstart <- init.state(
+ object,
+ params=if (transform) {
+ partrans(object,params,dir="forward")
+ } else {
+ params
+ }
+ )
statenames <- rownames(xstart)
nvars <- nrow(xstart)
@@ -141,17 +172,18 @@
)
}
- X <- matrix(data=x,nrow=nvars,ncol=Np*ntries)
- rownames(X) <- statenames
- P <- matrix(data=params,nrow=npars,ncol=Np*ntries)
- rownames(P) <- paramnames
## update mean of states at time nt as per L&W AGM (1)
tries <- rprocess(
object,
- xstart=X,
+ xstart=parmat(x,nrep=ntries),
times=times[c(nt,nt+1)],
- params=P
- )[,,2,drop=FALSE]
+ params=if (transform) {
+ partrans(object,params,dir="forward")
+ } else {
+ params
+ },
+ offset=1
+ )
dim(tries) <- c(nvars,Np,ntries,1)
mu <- apply(tries,c(1,2,4),mean)
rownames(mu) <- statenames
@@ -164,7 +196,11 @@
y=object at data[,nt,drop=FALSE],
x=mu,
times=times[nt+1],
- params=m
+ params=if (transform) {
+ partrans(object,m,dir="forward")
+ } else {
+ m
+ }
)
storeForEvidence1 <- log(sum(g))
## sample indices -- From L&W AGM (2)
@@ -190,13 +226,21 @@
stop(error.prefix,"extreme particle depletion",call.=FALSE)
params[estind,] <- m[estind,]+t(pvec)
+ if (transform)
+ tparams <- partrans(object,params,dir="forward")
+
## sample current state vector x^(g)_(t+1) as per L&W AGM (4)
X <- rprocess(
object,
xstart=x[,k,drop=FALSE],
times=times[c(nt,nt+1)],
- params=params
- )[,,2,drop=FALSE]
+ params=if (transform) {
+ tparams
+ } else {
+ params
+ },
+ offset=1
+ )
## evaluate likelihood of observation given X (from L&W AGM (4))
numer <- dmeasure(
@@ -204,7 +248,11 @@
y=object at data[,nt,drop=FALSE],
x=X,
times=times[nt+1],
- params=params
+ params=if (transform) {
+ tparams
+ } else {
+ params
+ }
)
## evaluate weights as per L&W AGM (5)
@@ -281,16 +329,62 @@
seed <- save.seed
}
- list(
- post=params,
- prior=prior,
- eff.sample.size=eff.sample.size,
- smooth=smooth,
- seed=seed,
- nfail=nfail,
- cond.log.evidence=evidence,
- log.evidence=sum(evidence),
- weights=weights
- )
+ ## if (transform) {
+ ## params <- partrans(object,params,dir="forward")
+ ## prior <- partrans(object,prior,dir="forward")
+ ## }
+
+ new(
+ "bsmcd.pomp",
+ object,
+ transform=transform,
+ post=params,
+ prior=prior,
+ est=as.character(est),
+ eff.sample.size=eff.sample.size,
+ smooth=smooth,
+ seed=as.integer(seed),
+ nfail=as.integer(nfail),
+ cond.log.evidence=evidence,
+ log.evidence=sum(evidence),
+ weights=weights
+ )
}
)
+
+setMethod("$",signature(x="bsmcd.pomp"),function (x,name) slot(x,name))
+
+bsmc.plot <- function (prior, post, est, nbreaks, thin, ...) {
+ if (missing(thin)) thin <- Inf
+ prior <- t(prior[est,sample.int(n=nrow(prior),size=min(thin,nrow(prior)))])
+ post <- t(post[est,sample.int(n=nrow(post),size=min(thin,nrow(post)))])
+ all <- rbind(prior,post)
+ pairs(
+ all,
+ labels=est,
+ panel=function (x, y, ...) { ## prior, posterior pairwise scatterplot
+ op <- par(new=TRUE)
+ on.exit(par(op))
+ i <- which(x[1]==all[1,])
+ j <- which(y[1]==all[1,])
+ points(x,y,pch=20,col=rgb(0.85,0.85,0.85,0.1),xlim=range(all[,i]),ylim=range(all[,j]))
+ points(post[,i],post[,j],pch=20,col=rgb(0,0,1,0.01))
+ },
+ diag.panel=function (x, ...) { ## marginal posterior histogram
+ i <- which(x[1]==all[1,])
+ breaks <- hist(c(post[,i],prior[,i]),breaks=nbreaks,plot=FALSE)$breaks
+ y1 <- hist(post[,i],breaks=breaks,plot=FALSE)$counts
+ usr <- par('usr')
+ op <- par(usr=c(usr[1:2],0,1.5*max(y1)))
+ on.exit(par(op))
+ rect(head(breaks,-1),0,tail(breaks,-1),y1,col=rgb(0,0,1,0.5),border=NA,...)
+ }
+ )
+}
+
+setMethod(
+ "plot",
+ signature(x="bsmcd.pomp"),
+ function (x, ..., thin, breaks) bsmc.plot(prior=x at prior,post=x at post,est=x at est,
+ nbreaks=breaks,thin=thin,...)
+ )
Modified: pkg/R/mif-class.R
===================================================================
--- pkg/R/mif-class.R 2012-03-21 19:52:13 UTC (rev 628)
+++ pkg/R/mif-class.R 2012-03-27 11:57:41 UTC (rev 629)
@@ -3,6 +3,7 @@
'mif',
contains='pfilterd.pomp',
representation=representation(
+ transform = "logical",
ivps = 'character',
pars = 'character',
Nmif = 'integer',
Modified: pkg/R/mif-methods.R
===================================================================
--- pkg/R/mif-methods.R 2012-03-21 19:52:13 UTC (rev 628)
+++ pkg/R/mif-methods.R 2012-03-27 11:57:41 UTC (rev 629)
@@ -20,7 +20,7 @@
partrans(
object,
params=t(object at conv.rec[,pars.proper]),
- dir="inverse"
+ dir="forward"
)
),
object at conv.rec[,pars.improper]
Modified: pkg/R/mif.R
===================================================================
--- pkg/R/mif.R 2012-03-21 19:52:13 UTC (rev 628)
+++ pkg/R/mif.R 2012-03-27 11:57:41 UTC (rev 629)
@@ -39,14 +39,20 @@
rw.sd,
Np, cooling.factor, var.factor, ic.lag,
method, tol, max.fail,
- verbose, .ndone) {
+ verbose, transform, .ndone) {
+ transform <- as.logical(transform)
+
if (length(start)==0)
stop(
"mif error: ",sQuote("start")," must be specified if ",
sQuote("coef(object)")," is NULL",
call.=FALSE
)
+
+ if (transform)
+ start <- partrans(object,start,dir="inverse")
+
start.names <- names(start)
if (missing(start.names))
stop("mif error: ",sQuote("start")," must be a named vector",call.=FALSE)
@@ -237,7 +243,8 @@
save.states=FALSE,
save.params=FALSE,
.rw.sd=sigma.n[pars],
- verbose=verbose
+ verbose=verbose,
+ transform=transform
),
silent=FALSE
)
@@ -274,9 +281,14 @@
}
+ ## back transform the parameter estimate if necessary
+ if (transform)
+ theta <- partrans(pfp,theta,dir="forward")
+
new(
"mif",
pfp,
+ transform=transform,
params=theta,
ivps=ivps,
pars=pars,
@@ -303,8 +315,11 @@
Np, ic.lag, var.factor, cooling.factor,
weighted, method = c("mif","unweighted","fp"),
tol = 1e-17, max.fail = 0,
- verbose = getOption("verbose"), ...) {
+ verbose = getOption("verbose"),
+ transform = FALSE, ...) {
+ transform <- as.logical(transform)
+
if (missing(start)) start <- coef(object)
if (missing(rw.sd))
stop("mif error: ",sQuote("rw.sd")," must be specified",call.=FALSE)
@@ -338,21 +353,7 @@
}
if (missing(particles)) { # use default: normal distribution
- particles <- function (Np, center, sd, ...) {
- matrix(
- data=rnorm(
- n=Np*length(center),
- mean=center,
- sd=sd
- ),
- nrow=length(center),
- ncol=Np,
- dimnames=list(
- names(center),
- NULL
- )
- )
- }
+ particles <- default.pomp.particles.fun
} else {
particles <- match.fun(particles)
if (!all(c('Np','center','sd','...')%in%names(formals(particles))))
@@ -381,6 +382,7 @@
tol=tol,
max.fail=max.fail,
verbose=verbose,
+ transform=transform,
.ndone=0
)
@@ -398,8 +400,11 @@
Np, ic.lag, var.factor, cooling.factor,
weighted, method = c("mif","unweighted","fp"),
tol = 1e-17, max.fail = 0,
- verbose = getOption("verbose"), ...) {
+ verbose = getOption("verbose"),
+ transform = FALSE, ...) {
+ transform <- as.logical(transform)
+
if (missing(start)) start <- coef(object)
if (missing(rw.sd))
stop("mif error: ",sQuote("rw.sd")," must be specified",call.=FALSE)
@@ -462,6 +467,7 @@
tol=tol,
max.fail=max.fail,
verbose=verbose,
+ transform=transform,
.ndone=0
)
}
@@ -477,7 +483,8 @@
Np, ic.lag, var.factor, cooling.factor,
weighted, method = c("mif","unweighted","fp"),
tol = 1e-17, max.fail = 0,
- verbose = getOption("verbose"), ...) {
+ verbose = getOption("verbose"),
+ transform, ...) {
if (missing(Nmif)) Nmif <- object at Nmif
if (missing(start)) start <- coef(object)
@@ -490,6 +497,8 @@
if (missing(var.factor)) var.factor <- object at var.factor
if (missing(cooling.factor)) cooling.factor <- object at cooling.factor
if (missing(tol)) tol <- object at tol
+ if (missing(transform)) transform <- object at transform
+ transform <- as.logical(transform)
method <- match.arg(method)
if (!missing(weighted)) {
@@ -523,6 +532,7 @@
tol=tol,
max.fail=max.fail,
verbose=verbose,
+ transform=transform,
.ndone=0
)
}
@@ -538,7 +548,8 @@
Np, ic.lag, var.factor, cooling.factor,
weighted, method = c("mif","unweighted","fp"),
tol = 1e-17, max.fail = 0,
- verbose = getOption("verbose"), ...) {
+ verbose = getOption("verbose"),
+ transform, ...) {
ndone <- object at Nmif
if (missing(start)) start <- coef(object)
@@ -551,6 +562,8 @@
if (missing(var.factor)) var.factor <- object at var.factor
if (missing(cooling.factor)) cooling.factor <- object at cooling.factor
if (missing(tol)) tol <- object at tol
+ if (missing(transform)) transform <- object at transform
+ transform <- as.logical(transform)
method <- match.arg(method)
if (!missing(weighted)) {
@@ -584,6 +597,7 @@
tol=tol,
max.fail=max.fail,
verbose=verbose,
+ transform=transform,
.ndone=ndone
)
Modified: pkg/R/nlf-objfun.R
===================================================================
--- pkg/R/nlf-objfun.R 2012-03-21 19:52:13 UTC (rev 628)
+++ pkg/R/nlf-objfun.R 2012-03-27 11:57:41 UTC (rev 629)
@@ -1,4 +1,4 @@
-NLF.LQL <- function (params.fitted, object, params, par.index,
+NLF.LQL <- function (params.fitted, object, params, par.index, transform.params = FALSE,
times, t0, lags, period, tensor, seed = NULL, transform = identity,
nrbf = 4, verbose = FALSE,
bootstrap = FALSE, bootsamp = NULL) {
@@ -9,10 +9,13 @@
### so a large NEGATIVE value is used to flag bad parameters
###>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ transform.params <- as.logical(transform.params)
+
FAILED = -99999999999
params[par.index] <- params.fitted
- params <- as.matrix(params)
+ if (transform.params)
+ params <- partrans(object,params,dir="forward")
## Need to extract number of state variables (nvar) from pomp object
## Need to include simulation times in problem specification
@@ -20,7 +23,7 @@
## Version 0.1, 3 Dec. 2007, Bruce E. Kendall & Stephen P. Ellner
## Version 0.2, May 2008, Stephen P. Ellner
- data.ts <- object at data
+ data.ts <- obs(object)
y <- try(
simulate(object,times=times,t0=t0,params=params,seed=seed,obs=TRUE,states=FALSE),
@@ -60,30 +63,5 @@
LQL
}
-
-
-nlf.objfun <- function (params.fitted, object, params, par.index,
- times, t0, lags, period, tensor, seed, transform = function(x)x,
- nrbf = 4, verbose = FALSE, bootstrap = FALSE, bootsamp = NULL)
-{
- -sum(
- NLF.LQL(
- params.fitted=params.fitted,
- object=object,
- params=params,
- par.index=par.index,
- times=times,
- t0=t0,
- lags=lags,
- period=period,
- tensor=tensor,
- seed=seed,
- transform=transform,
- nrbf=nrbf,
- verbose=verbose,
- bootstrap=bootstrap,
- bootsamp=bootsamp
- ),
- na.rm=T
- )
-}
+nlf.objfun <- function (...)
+ -sum(NLF.LQL(...),na.rm=TRUE)
Modified: pkg/R/nlf.R
===================================================================
--- pkg/R/nlf.R 2012-03-21 19:52:13 UTC (rev 628)
+++ pkg/R/nlf.R 2012-03-27 11:57:41 UTC (rev 629)
@@ -1,12 +1,12 @@
nlf <- function (object, start, est, lags,
period = NA, tensor = FALSE,
nconverge = 1000, nasymp = 1000,
- seed = 1066, transform = function(x)x,
+ seed = 1066, transform = identity,
nrbf = 4, method = "subplex",
skip.se = FALSE, verbose = FALSE, gr = NULL,
bootstrap = FALSE, bootsamp = NULL,
lql.frac = 0.1, se.par.frac = 0.1,
- eval.only = FALSE, ...) {
+ eval.only = FALSE, transform.params = FALSE, ...) {
## Fit a POMP object using NLF
## v. 0.1, 3 Dec. 2007
@@ -25,33 +25,39 @@
if (!is(object,'pomp'))
stop("'object' must be a 'pomp' object")
- if (!is.function(transform))
- stop(sQuote("transform")," must be a function!")
+ transform <- match.fun(transform)
- if (eval.only) est <- 1
+ if (eval.only) est <- 1L
- if (is.character(est)) {
- if (!all(est%in%names(start)))
- stop(sQuote("nlf")," error: parameters named in ",sQuote("est")," must exist in ",sQuote("start"))
+ if (missing(start)) start <- coef(object)
- par.index <- which(names(start)%in%est)
+ transform.params <- as.logical(transform.params)
+ if (transform.params)
+ params <- partrans(object,start,dir="inverse")
+ else
+ params <- start
+ if (is.character(est)) {
+ if (!all(est%in%names(params)))
+ stop("parameters named in ",sQuote("est")," must exist in ",sQuote("start"))
+ par.index <- which(names(params)%in%est)
} else if (is.numeric(est)) {
est <- as.integer(est)
- if (any((est<1)|(est>length(start))))
- stop(sQuote("nlf")," error: trying to estimate parameters that don't exist!")
- par.index <- as.integer(est)
+ if (any((est<1)|(est>length(params))))
+ stop("indices in ",sQuote("est")," are not appropriate")
+ par.index <- est
}
- if (!(is.numeric(lql.frac)&&(lql.frac>0)&&(lql.frac<1)))
- stop(sQuote("nlf")," error: ",sQuote("lql.frac")," must be in (0,1)")
+ guess <- params[par.index]
- if (!(is.numeric(se.par.frac)&&(se.par.frac>0)&&(se.par.frac<1)))
- stop(sQuote("nlf")," error: ",sQuote("se.par.frac")," must be in (0,1)")
+ lql.frac <- as.numeric(lql.frac)
+ if ((lql.frac<=0)||(lql.frac>=1))
+ stop(sQuote("lql.frac")," must be in (0,1)")
+
+ se.par.frac <- as.numeric(se.par.frac)
+ if ((se.par.frac<=0)||(se.par.frac>=1))
+ stop(sQuote("se.par.frac")," must be in (0,1)")
- params <- start
- guess <- params[par.index]
-
dt.tol <- 1e-3
times <- time(object,t0=FALSE)
t0 <- timezero(object)
@@ -73,6 +79,7 @@
object=object,
params=params,
par.index=par.index,
+ transform.params=transform.params,
times=times,
t0=t0,
lags=lags,
@@ -95,6 +102,7 @@
object=object,
params=params,
par.index=par.index,
+ transform.params=transform.params,
times=times,
t0=t0,
lags=lags,
@@ -117,6 +125,7 @@
object=object,
params=params,
par.index=par.index,
+ transform.params=transform.params,
times=times,
t0=t0,
lags=lags,
@@ -135,15 +144,19 @@
opt$est <- est
opt$value <- -opt$value
params[par.index] <- opt$par
- opt$params <- params
+ opt$params <- if (transform.params) partrans(object,params,dir="forward") else params
opt$par <- NULL
- if (!skip.se) { ## compute estimated Variance-Covariance matrix of fitted parameters
+ if (!skip.se) { ## compute estimated Variance-Covariance matrix of fitted parameters
fitted <- params[par.index]
nfitted <- length(fitted)
Jhat <- matrix(0,nfitted,nfitted)
Ihat <- Jhat
- f0 <- NLF.LQL(fitted,object=object, params=params, par.index=par.index,
+ f0 <- NLF.LQL(fitted,
+ object=object,
+ params=params,
+ par.index=par.index,
+ transform.params=transform.params,
times=times, t0=t0,
lags=lags, period=period, tensor=tensor, seed=seed,
transform=transform, nrbf=4,
@@ -165,24 +178,28 @@
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,
+ transform.params=transform.params,
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,
+ transform.params=transform.params,
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,
+ transform.params=transform.params,
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,
+ transform.params=transform.params,
times=times, t0=t0, lags=lags, period=period, tensor=tensor,
seed=seed, transform=transform, nrbf=4,
verbose=FALSE),na.rm=T)
@@ -201,12 +218,14 @@
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,
+ transform.params=transform.params,
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,
+ transform.params=transform.params,
times=times, t0=t0, lags=lags, period=period, tensor=tensor,
seed=seed, transform=transform, nrbf=4,
verbose=FALSE)
@@ -217,6 +236,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,
+ transform.params=transform.params,
times=times, t0=t0, lags=lags, period=period, tensor=tensor,
seed=seed, transform=transform, nrbf=4,
verbose=FALSE)
@@ -236,6 +256,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,
+ transform.params=transform.params,
times=times, t0=t0, lags=lags, period=period, tensor=tensor,
seed=seed, transform=transform, nrbf=4,
verbose=FALSE),na.rm=T)
@@ -244,6 +265,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,
+ transform.params=transform.params,
times=times, t0=t0, lags=lags, period=period, tensor=tensor,
seed=seed, transform=transform, nrbf=4,
verbose=FALSE),na.rm=T)
@@ -252,6 +274,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,
+ transform.params=transform.params,
times=times, t0=t0, lags=lags, period=period, tensor=tensor,
seed=seed, transform=transform, nrbf=4,
verbose=FALSE),na.rm=T)
@@ -260,6 +283,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,
+ transform.params=transform.params,
times=times, t0=t0, lags=lags, period=period, tensor=tensor,
seed=seed, transform=transform, nrbf=4,
verbose=FALSE),na.rm=T)
@@ -272,13 +296,14 @@
Ihat[j,i] <- Ihat[i,j]
}
}
+ opt$transform.params <- transform.params
opt$Jhat <- Jhat
opt$Ihat <- Ihat
negJinv <- -solve(Jhat)
Qhat <- negJinv%*%Ihat%*%negJinv
opt$Qhat <- Qhat
opt$se <- sqrt(diag(Qhat))/sqrt(npts)
- names(opt$se) <- names(start)[par.index]
+ names(opt$se) <- names(params)[par.index]
opt$npts <- npts
}
Modified: pkg/R/pfilter.R
===================================================================
--- pkg/R/pfilter.R 2012-03-21 19:52:13 UTC (rev 628)
+++ pkg/R/pfilter.R 2012-03-27 11:57:41 UTC (rev 629)
@@ -25,8 +25,11 @@
tol, max.fail,
pred.mean, pred.var, filter.mean,
.rw.sd, seed, verbose,
- save.states, save.params) {
+ save.states, save.params,
+ transform) {
+ transform <- as.logical(transform)
+
if (missing(seed)) seed <- NULL
if (!is.null(seed)) {
if (!exists(".Random.seed",where=.GlobalEnv)) { # need to initialize the RNG
@@ -83,7 +86,14 @@
if (is.null(paramnames))
stop(sQuote("pfilter")," error: ",sQuote("params")," must have rownames",call.=FALSE)
- x <- init.state(object,params=params)
+ x <- init.state(
+ object,
+ params=if (transform) {
+ partrans(object,params,dir="forward")
+ } else {
+ params
+ }
+ )
statenames <- rownames(x)
nvars <- nrow(x)
@@ -160,13 +170,16 @@
for (nt in seq_len(ntimes)) {
+ ## transform the parameters if necessary
+ if (transform) tparams <- partrans(object,params,dir="forward")
+
## advance the state variables according to the process model
X <- try(
rprocess(
object,
xstart=x,
times=times[c(nt,nt+1)],
- params=params,
+ params=if (transform) tparams else params,
offset=1
),
silent=FALSE
@@ -202,7 +215,7 @@
y=object at data[,nt,drop=FALSE],
x=X,
times=times[nt+1],
- params=params,
+ params=if (transform) tparams else params,
log=FALSE
),
silent=FALSE
@@ -213,7 +226,8 @@
stop(sQuote("pfilter")," error: ",sQuote("dmeasure")," returns non-finite value",call.=FALSE)
}
- ## prediction mean, prediction variance, filtering mean, effective sample size, log-likelihood
+ ## compute prediction mean, prediction variance, filtering mean,
+ ## effective sample size, log-likelihood
## also do resampling if filtering has not failed
xx <- try(
.Call(
@@ -317,7 +331,8 @@
save.states=save.states,
save.params=save.params,
seed=seed,
- verbose=verbose
+ verbose=verbose,
+ transform=FALSE
)
}
)
@@ -351,7 +366,8 @@
save.states=save.states,
save.params=save.params,
seed=seed,
- verbose=verbose
+ verbose=verbose,
+ transform=FALSE
)
}
)
Modified: pkg/R/pmcmc.R
===================================================================
--- pkg/R/pmcmc.R 2012-03-21 19:52:13 UTC (rev 628)
+++ pkg/R/pmcmc.R 2012-03-27 11:57:41 UTC (rev 629)
@@ -4,6 +4,7 @@
contains='pfilterd.pomp',
representation(
pars = 'character',
+ transform = 'logical',
Nmcmc = 'integer',
dprior = 'function',
hyperparams = 'list',
@@ -31,14 +32,16 @@
rw.sd, Np,
hyperparams,
tol, max.fail,
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 629
More information about the pomp-commits
mailing list