[Pomp-commits] r760 - branches/mif2/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Aug 7 16:40:18 CEST 2012
Author: nxdao2000
Date: 2012-08-07 16:40:18 +0200 (Tue, 07 Aug 2012)
New Revision: 760
Modified:
branches/mif2/R/pfilter.R
Log:
new update
Modified: branches/mif2/R/pfilter.R
===================================================================
--- branches/mif2/R/pfilter.R 2012-08-07 14:38:18 UTC (rev 759)
+++ branches/mif2/R/pfilter.R 2012-08-07 14:40:18 UTC (rev 760)
@@ -1,400 +1,401 @@
## particle filtering codes
setClass(
- "pfilterd.pomp",
- contains="pomp",
- representation=representation(
- pred.mean="array",
- pred.var="array",
- filter.mean="array",
- paramMatrix = "array",
- eff.sample.size="numeric",
- cond.loglik="numeric",
- saved.states="list",
- saved.params="list",
- seed="integer",
- Np="integer",
- tol="numeric",
- nfail="integer",
- loglik="numeric"
- )
- )
+ "pfilterd.pomp",
+ contains="pomp",
+ representation=representation(
+ pred.mean="array",
+ pred.var="array",
+ filter.mean="array",
+ paramMatrix = "array",
+ eff.sample.size="numeric",
+ cond.loglik="numeric",
+ saved.states="list",
+ saved.params="list",
+ seed="integer",
+ Np="integer",
+ tol="numeric",
+ nfail="integer",
+ loglik="numeric"
+ )
+)
## question: when pfilter.internal is called by mif, do we need to compute the prediction means and variances of the state variables each time, or only at the end?
pfilter.internal <- function (object, params, Np,
- tol, max.fail,
- pred.mean, pred.var, filter.mean, paramMatrix, cooling.scalar,cooling.m,option,
- .rw.sd, seed, verbose,
- 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
- runif(1)
- }
- save.seed <- get(".Random.seed",pos=.GlobalEnv)
- set.seed(seed)
- }
-
- if (length(params)==0)
- stop(sQuote("pfilter")," error: ",sQuote("params")," must be specified",call.=FALSE)
-
- if (missing(tol))
- stop(sQuote("pfilter")," error: ",sQuote("tol")," must be specified",call.=FALSE)
- one.par <- FALSE
- times <- time(object,t0=TRUE)
- ntimes <- length(times)-1
-
- if (missing(Np))
- Np <- NCOL(params)
- if (is.function(Np)) {
- Np <- try(
- vapply(seq.int(from=0,to=ntimes,by=1),Np,numeric(1)),
- silent=FALSE
- )
- if (inherits(Np,"try-error"))
- stop("if ",sQuote("Np")," is a function, it must return a single positive integer")
- }
- if (length(Np)==1)
- Np <- rep(Np,times=ntimes+1)
- else if (length(Np)!=(ntimes+1))
- stop(sQuote("Np")," must have length 1 or length ",ntimes+1)
- if (any(Np<=0))
- stop("number of particles, ",sQuote("Np"),", must always be positive")
- if (!is.numeric(Np))
- stop(sQuote("Np")," must be a number, a vector of numbers, or a function")
- Np <- as.integer(Np)
-
- if (is.null(dim(params))) {
- one.par <- TRUE # there is only one parameter vector
- coef(object) <- params # set params slot to the parameters
- params <- matrix(
- params,
- nrow=length(params),
- ncol=Np[1],
- dimnames=list(
- names(params),
- NULL
- )
- )
- }
- paramnames <- rownames(params)
- if (is.null(paramnames))
- stop(sQuote("pfilter")," error: ",sQuote("params")," must have rownames",call.=FALSE)
-
- x <- init.state(
- object,
- params=if (transform) {
- partrans(object,params,dir="forward")
- } else {
- params
- }
- )
- statenames <- rownames(x)
- nvars <- nrow(x)
-
- ## set up storage for saving samples from filtering distributions
- if (save.states)
- xparticles <- vector(mode="list",length=ntimes)
- else
- xparticles <- list()
- if (save.params)
- pparticles <- vector(mode="list",length=ntimes)
- else
- pparticles <- list()
-
- random.walk <- !missing(.rw.sd)
- if (random.walk) {
- rw.names <- names(.rw.sd)
- if (is.null(rw.names)||!is.numeric(.rw.sd))
- stop(sQuote("pfilter")," error: ",sQuote(".rw.sd")," must be a named vector",call.=FALSE)
- if (any(!(rw.names%in%paramnames)))
- stop(
- sQuote("pfilter")," error: the rownames of ",
- sQuote("params")," must include all of the names of ",
- sQuote(".rw.sd"),"",call.=FALSE
- )
- sigma <- .rw.sd
- } else {
- rw.names <- character(0)
- sigma <- NULL
- }
-
- loglik <- rep(NA,ntimes)
- eff.sample.size <- numeric(ntimes)
- nfail <- 0
- npars <- length(rw.names)
-
- ## set up storage for prediction means, variances, etc.
- if (pred.mean)
- pred.m <- matrix(
- data=0,
- nrow=nvars+npars,
- ncol=ntimes,
- dimnames=list(c(statenames,rw.names),NULL)
- )
- else
- pred.m <- array(dim=c(0,0))
-
- if (pred.var)
- pred.v <- matrix(
- data=0,
- nrow=nvars+npars,
- ncol=ntimes,
- dimnames=list(c(statenames,rw.names),NULL)
- )
- else
- pred.v <- array(dim=c(0,0))
-
- if (filter.mean)
- if (random.walk)
- filt.m <- matrix(
- data=0,
- nrow=nvars+length(paramnames),
- ncol=ntimes,
- dimnames=list(c(statenames,paramnames),NULL)
- )
- else
- filt.m <- matrix(
- data=0,
- nrow=nvars,
- ncol=ntimes,
- dimnames=list(statenames,NULL)
- )
- else
- filt.m <- array(dim=c(0,0))
- if (paramMatrix)
- { paramMatrix<-matrix(
- data=0,
- nrow=length(paramnames),
- ncol=ntimes,
- dimnames=list(paramnames,NULL)
- )
-
- }
- else
- paramMatrix <- array(dim=c(0,0))
+ tol, max.fail,
+ pred.mean, pred.var, filter.mean, paramMatrix, cooling.scalar,cooling.m, option,
+ .rw.sd, seed, verbose,
+ save.states, save.params,
+ transform) {
- for (nt in seq_len(ntimes)) {
- if (option=="mif2")
- { cool.sched <- try(mif.cooling2(cooling.scalar, nt , cooling.m, ntimes), silent = FALSE)
- if (inherits(cool.sched, "try-error"))
- stop("pfilter error: cooling schedule error", call. = FALSE)
- sigma1=sigma*cool.sched$alpha
- }
- else
- {
- sigma1=sigma
-
- }
- ## 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=if (transform) tparams else params,
- offset=1
- ),
- silent=FALSE
- )
- if (inherits(X,'try-error'))
- stop(sQuote("pfilter")," error: process simulation error",call.=FALSE)
-
- if (pred.var) { ## check for nonfinite state variables and parameters
- problem.indices <- unique(which(!is.finite(X),arr.ind=TRUE)[,1])
- if (length(problem.indices)>0) { # state variables
- stop(
- sQuote("pfilter")," error: non-finite state variable(s): ",
- paste(rownames(X)[problem.indices],collapse=', '),
- call.=FALSE
- )
- }
- if (random.walk) { # parameters (need to be checked only if 'random.walk=TRUE')
- problem.indices <- unique(which(!is.finite(params[rw.names,,drop=FALSE]),arr.ind=TRUE)[,1])
- if (length(problem.indices)>0) {
- stop(
- sQuote("pfilter")," error: non-finite parameter(s): ",
- paste(rw.names[problem.indices],collapse=', '),
- call.=FALSE
- )
- }
- }
- }
-
- ## determine the weights
- weights <- try(
- dmeasure(
- object,
- y=object at data[,nt,drop=FALSE],
- x=X,
- times=times[nt+1],
- params=if (transform) tparams else params,
- log=FALSE
- ),
- silent=FALSE
- )
- if (inherits(weights,'try-error'))
- stop(sQuote("pfilter")," error: error in calculation of weights",call.=FALSE)
- if (any(!is.finite(weights))) {
- stop(sQuote("pfilter")," error: ",sQuote("dmeasure")," returns non-finite value",call.=FALSE)
- }
-
- ## compute prediction mean, prediction variance, filtering mean,
- ## effective sample size, log-likelihood
- ## also do resampling if filtering has not failed
- xx <- try(
- .Call(
- pfilter_computations,
- X,params,Np[nt+1],
- random.walk,sigma1,
- pred.mean,pred.var,
- filter.mean,one.par,
- weights,tol
- ),
- silent=FALSE
- )
- if (inherits(xx,'try-error')) {
- stop(sQuote("pfilter")," error",call.=FALSE)
- }
- all.fail <- xx$fail
- loglik[nt] <- xx$loglik
- eff.sample.size[nt] <- xx$ess
-
- x <- xx$states
- params <- xx$params
-
- if (pred.mean)
- pred.m[,nt] <- xx$pm
- if (pred.var)
- pred.v[,nt] <- xx$pv
- if (filter.mean)
- filt.m[,nt] <- xx$fm
-
- if (all.fail) { ## all particles are lost
- nfail <- nfail+1
- if (verbose)
- message("filtering failure at time t = ",times[nt+1])
- if (nfail>max.fail)
- stop(sQuote("pfilter")," error: too many filtering failures",call.=FALSE)
- }
-
- if (save.states) {
- xparticles[[nt]] <- x
- }
-
- if (save.params) {
- pparticles[[nt]] <- params
- }
-
- if (verbose && (nt%%5==0))
- cat("pfilter timestep",nt,"of",ntimes,"finished\n")
-
- }
-
- if (!is.null(seed)) {
- assign(".Random.seed",save.seed,pos=.GlobalEnv)
- seed <- save.seed
- }
- if(length(paramMatrix)!=0)
- paramMatrix <- params
- new(
- "pfilterd.pomp",
- object,
- pred.mean=pred.m,
- pred.var=pred.v,
- filter.mean=filt.m,
- paramMatrix = paramMatrix,
- eff.sample.size=eff.sample.size,
- cond.loglik=loglik,
- saved.states=xparticles,
- saved.params=pparticles,
- seed=as.integer(seed),
- Np=as.integer(Np),
- tol=tol,
- nfail=as.integer(nfail),
- loglik=sum(loglik)
- )
+ transform <- as.logical(transform)
+
+ if (missing(seed)) seed <- NULL
+ if (!is.null(seed)) {
+ if (!exists(".Random.seed",where=.GlobalEnv)) { # need to initialize the RNG
+ runif(1)
+ }
+ save.seed <- get(".Random.seed",pos=.GlobalEnv)
+ set.seed(seed)
+ }
+
+ if (length(params)==0)
+ stop(sQuote("pfilter")," error: ",sQuote("params")," must be specified",call.=FALSE)
+
+ if (missing(tol))
+ stop(sQuote("pfilter")," error: ",sQuote("tol")," must be specified",call.=FALSE)
+
+ one.par <- FALSE
+ times <- time(object,t0=TRUE)
+ ntimes <- length(times)-1
+
+ if (missing(Np))
+ Np <- NCOL(params)
+ if (is.function(Np)) {
+ Np <- try(
+ vapply(seq.int(from=0,to=ntimes,by=1),Np,numeric(1)),
+ silent=FALSE
+ )
+ if (inherits(Np,"try-error"))
+ stop("if ",sQuote("Np")," is a function, it must return a single positive integer")
+ }
+ if (length(Np)==1)
+ Np <- rep(Np,times=ntimes+1)
+ else if (length(Np)!=(ntimes+1))
+ stop(sQuote("Np")," must have length 1 or length ",ntimes+1)
+ if (any(Np<=0))
+ stop("number of particles, ",sQuote("Np"),", must always be positive")
+ if (!is.numeric(Np))
+ stop(sQuote("Np")," must be a number, a vector of numbers, or a function")
+ Np <- as.integer(Np)
+
+ if (is.null(dim(params))) {
+ one.par <- TRUE # there is only one parameter vector
+ coef(object) <- params # set params slot to the parameters
+ params <- matrix(
+ params,
+ nrow=length(params),
+ ncol=Np[1],
+ dimnames=list(
+ names(params),
+ NULL
+ )
+ )
+ }
+ paramnames <- rownames(params)
+ if (is.null(paramnames))
+ stop(sQuote("pfilter")," error: ",sQuote("params")," must have rownames",call.=FALSE)
+
+ x <- init.state(
+ object,
+ params=if (transform) {
+ partrans(object,params,dir="forward")
+ } else {
+ params
+ }
+ )
+ statenames <- rownames(x)
+ nvars <- nrow(x)
+
+ ## set up storage for saving samples from filtering distributions
+ if (save.states)
+ xparticles <- vector(mode="list",length=ntimes)
+ else
+ xparticles <- list()
+ if (save.params)
+ pparticles <- vector(mode="list",length=ntimes)
+ else
+ pparticles <- list()
+
+ random.walk <- !missing(.rw.sd)
+ if (random.walk) {
+ rw.names <- names(.rw.sd)
+ if (is.null(rw.names)||!is.numeric(.rw.sd))
+ stop(sQuote("pfilter")," error: ",sQuote(".rw.sd")," must be a named vector",call.=FALSE)
+ if (any(!(rw.names%in%paramnames)))
+ stop(
+ sQuote("pfilter")," error: the rownames of ",
+ sQuote("params")," must include all of the names of ",
+ sQuote(".rw.sd"),"",call.=FALSE
+ )
+ sigma <- .rw.sd
+ } else {
+ rw.names <- character(0)
+ sigma <- NULL
+ }
+
+ loglik <- rep(NA,ntimes)
+ eff.sample.size <- numeric(ntimes)
+ nfail <- 0
+ npars <- length(rw.names)
+
+ ## set up storage for prediction means, variances, etc.
+ if (pred.mean)
+ pred.m <- matrix(
+ data=0,
+ nrow=nvars+npars,
+ ncol=ntimes,
+ dimnames=list(c(statenames,rw.names),NULL)
+ )
+ else
+ pred.m <- array(dim=c(0,0))
+
+ if (pred.var)
+ pred.v <- matrix(
+ data=0,
+ nrow=nvars+npars,
+ ncol=ntimes,
+ dimnames=list(c(statenames,rw.names),NULL)
+ )
+ else
+ pred.v <- array(dim=c(0,0))
+
+ if (filter.mean)
+ if (random.walk)
+ filt.m <- matrix(
+ data=0,
+ nrow=nvars+length(paramnames),
+ ncol=ntimes,
+ dimnames=list(c(statenames,paramnames),NULL)
+ )
+ else
+ filt.m <- matrix(
+ data=0,
+ nrow=nvars,
+ ncol=ntimes,
+ dimnames=list(statenames,NULL)
+ )
+ else
+ filt.m <- array(dim=c(0,0))
+ if (paramMatrix)
+ { paramMatrix<-matrix(
+ data=0,
+ nrow=length(paramnames),
+ ncol=ntimes,
+ dimnames=list(paramnames,NULL)
+ )
+
+ }
+ else
+ paramMatrix <- array(dim=c(0,0))
+
+
+
+ for (nt in seq_len(ntimes)) {
+ if (option=="mif2")
+ { cool.sched <- try(mif.cooling2(cooling.scalar, nt , cooling.m, ntimes), silent = FALSE)
+ if (inherits(cool.sched, "try-error"))
+ stop("pfilter error: cooling schedule error", call. = FALSE)
+ sigma1=sigma*cool.sched$alpha
+
+ }
+ else
+ sigma1=sigma
+ ## 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=if (transform) tparams else params,
+ offset=1
+ ),
+ silent=FALSE
+ )
+ if (inherits(X,'try-error'))
+ stop(sQuote("pfilter")," error: process simulation error",call.=FALSE)
+
+ if (pred.var) { ## check for nonfinite state variables and parameters
+ problem.indices <- unique(which(!is.finite(X),arr.ind=TRUE)[,1])
+ if (length(problem.indices)>0) { # state variables
+ stop(
+ sQuote("pfilter")," error: non-finite state variable(s): ",
+ paste(rownames(X)[problem.indices],collapse=', '),
+ call.=FALSE
+ )
+ }
+ if (random.walk) { # parameters (need to be checked only if 'random.walk=TRUE')
+ problem.indices <- unique(which(!is.finite(params[rw.names,,drop=FALSE]),arr.ind=TRUE)[,1])
+ if (length(problem.indices)>0) {
+ stop(
+ sQuote("pfilter")," error: non-finite parameter(s): ",
+ paste(rw.names[problem.indices],collapse=', '),
+ call.=FALSE
+ )
+ }
+ }
+ }
+
+ ## determine the weights
+ weights <- try(
+ dmeasure(
+ object,
+ y=object at data[,nt,drop=FALSE],
+ x=X,
+ times=times[nt+1],
+ params=if (transform) tparams else params,
+ log=FALSE
+ ),
+ silent=FALSE
+ )
+ if (inherits(weights,'try-error'))
+ stop(sQuote("pfilter")," error: error in calculation of weights",call.=FALSE)
+ if (any(!is.finite(weights))) {
+ stop(sQuote("pfilter")," error: ",sQuote("dmeasure")," returns non-finite value",call.=FALSE)
+ }
+
+ ## compute prediction mean, prediction variance, filtering mean,
+ ## effective sample size, log-likelihood
+ ## also do resampling if filtering has not failed
+ xx <- try(
+ .Call(
+ pfilter_computations,
+ X,params,Np[nt+1],
+ random.walk,sigma1,
+ pred.mean,pred.var,
+ filter.mean,one.par,
+ weights,tol
+ ),
+ silent=FALSE
+ )
+ if (inherits(xx,'try-error')) {
+ stop(sQuote("pfilter")," error",call.=FALSE)
+ }
+ all.fail <- xx$fail
+ loglik[nt] <- xx$loglik
+ eff.sample.size[nt] <- xx$ess
+
+ x <- xx$states
+ params <- xx$params
+
+ if (pred.mean)
+ pred.m[,nt] <- xx$pm
+ if (pred.var)
+ pred.v[,nt] <- xx$pv
+ if (filter.mean)
+ filt.m[,nt] <- xx$fm
+
+ if (all.fail) { ## all particles are lost
+ nfail <- nfail+1
+ if (verbose)
+ message("filtering failure at time t = ",times[nt+1])
+ if (nfail>max.fail)
+ stop(sQuote("pfilter")," error: too many filtering failures",call.=FALSE)
+ }
+
+ if (save.states) {
+ xparticles[[nt]] <- x
+ }
+
+ if (save.params) {
+ pparticles[[nt]] <- params
+ }
+
+ if (verbose && (nt%%5==0))
+ cat("pfilter timestep",nt,"of",ntimes,"finished\n")
+
+ }
+
+ if (!is.null(seed)) {
+ assign(".Random.seed",save.seed,pos=.GlobalEnv)
+ seed <- save.seed
+ }
+ if(length(paramMatrix)!=0)
+ paramMatrix <- params
+ new(
+ "pfilterd.pomp",
+ object,
+ pred.mean=pred.m,
+ pred.var=pred.v,
+ filter.mean=filt.m,
+ paramMatrix = paramMatrix,
+ eff.sample.size=eff.sample.size,
+ cond.loglik=loglik,
+ saved.states=xparticles,
+ saved.params=pparticles,
+ seed=as.integer(seed),
+ Np=as.integer(Np),
+ tol=tol,
+ nfail=as.integer(nfail),
+ loglik=sum(loglik)
+ )
}
## generic particle filter
setGeneric("pfilter",function(object,...)standardGeneric("pfilter"))
setMethod(
- "pfilter",
- signature=signature(object="pomp"),
- function (object, params, Np,
- tol = 1e-17,
- max.fail = 0,
- pred.mean = FALSE,
- pred.var = FALSE,
- filter.mean = FALSE,
- paramMatrix = FALSE,
- save.states = FALSE,
- save.params = FALSE,
- seed = NULL,
- verbose = getOption("verbose"),
- ...) {
- if (missing(params)) params <- coef(object)
- pfilter.internal(
- object=object,
- params=params,
- Np=Np,
- tol=tol,
- max.fail=max.fail,
- pred.mean=pred.mean,
- pred.var=pred.var,
- filter.mean=filter.mean,
- paramMatrix = paramMatrix,
- save.states=save.states,
- save.params=save.params,
- seed=seed,
- verbose=verbose,
- transform=FALSE
- )
- }
- )
+ "pfilter",
+ signature=signature(object="pomp"),
+ function (object, params, Np,
+ tol = 1e-17,
+ max.fail = 0,
+ pred.mean = FALSE,
+ pred.var = FALSE,
+ filter.mean = FALSE,
+ paramMatrix = FALSE,
+ save.states = FALSE,
+ save.params = FALSE,
+ seed = NULL,
+ verbose = getOption("verbose"),
+ ...) {
+ if (missing(params)) params <- coef(object)
+ pfilter.internal(
+ object=object,
+ params=params,
+ Np=Np,
+ tol=tol,
+ max.fail=max.fail,
+ pred.mean=pred.mean,
+ pred.var=pred.var,
+ filter.mean=filter.mean,
+ paramMatrix = paramMatrix,
+ save.states=save.states,
+ save.params=save.params,
+ seed=seed,
+ verbose=verbose,
+ transform=FALSE
+ )
+ }
+)
setMethod(
- "pfilter",
- signature=signature(object="pfilterd.pomp"),
- function (object, params, Np,
- tol,
- max.fail = 0,
- pred.mean = FALSE,
- pred.var = FALSE,
- filter.mean = FALSE,
- paramMatrix = FALSE,
- save.states = FALSE,
- save.params = FALSE,
- seed = NULL,
- verbose = getOption("verbose"),
- ...) {
- if (missing(params)) params <- coef(object)
- if (missing(Np)) Np <- object at Np
- if (missing(tol)) tol <- object at tol
- pfilter.internal(
- object=as(object,"pomp"),
- params=params,
- Np=Np,
- tol=tol,
- max.fail=max.fail,
- pred.mean=pred.mean,
- pred.var=pred.var,
- filter.mean=filter.mean,
- paramMatrix = paramMatrix,
- save.states=save.states,
- save.params=save.params,
- seed=seed,
- verbose=verbose,
- transform=FALSE
- )
- }
- )
+ "pfilter",
+ signature=signature(object="pfilterd.pomp"),
+ function (object, params, Np,
+ tol,
+ max.fail = 0,
+ pred.mean = FALSE,
+ pred.var = FALSE,
+ filter.mean = FALSE,
+ paramMatrix = FALSE,
+ save.states = FALSE,
+ save.params = FALSE,
+ seed = NULL,
+ verbose = getOption("verbose"),
+ ...) {
+ if (missing(params)) params <- coef(object)
+ if (missing(Np)) Np <- object at Np
+ if (missing(tol)) tol <- object at tol
+ pfilter.internal(
+ object=as(object,"pomp"),
+ params=params,
+ Np=Np,
+ tol=tol,
+ max.fail=max.fail,
+ pred.mean=pred.mean,
+ pred.var=pred.var,
+ filter.mean=filter.mean,
+ paramMatrix = paramMatrix,
+ save.states=save.states,
+ save.params=save.params,
+ seed=seed,
+ verbose=verbose,
+ transform=FALSE
+ )
+ }
+)
More information about the pomp-commits
mailing list