[Dream-commits] r20 - in pkg: R demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Mar 12 07:40:09 CET 2010
Author: josephguillaume
Date: 2010-03-12 07:40:09 +0100 (Fri, 12 Mar 2010)
New Revision: 20
Modified:
pkg/R/CompDensity.R
pkg/R/coef.dream.R
pkg/R/predict.dream.R
pkg/demo/FME.nonlinear.model.R
pkg/man/predict.dream.Rd
Log:
Added check that posterior density is strictly positive.
Allow newFUN with same parameters for predict, for when FUN returns posterior density to dream
Allow newFUN to return list of numerics. predict calculates CIs for each element of the list.
Modified: pkg/R/CompDensity.R
===================================================================
--- pkg/R/CompDensity.R 2010-03-05 03:49:04 UTC (rev 19)
+++ pkg/R/CompDensity.R 2010-03-12 06:40:09 UTC (rev 20)
@@ -45,6 +45,7 @@
## Model directly computes posterior density
posterior.density={
p <- modpred
+ if (any(modpred<=0)) stop("Posterior density returned by FUN should be strictly positive. Otherwise use logposterior.density?")
logp <- log(modpred)
},
## Model computes output simulation
Modified: pkg/R/coef.dream.R
===================================================================
--- pkg/R/coef.dream.R 2010-03-05 03:49:04 UTC (rev 19)
+++ pkg/R/coef.dream.R 2010-03-12 06:40:09 UTC (rev 20)
@@ -6,10 +6,19 @@
##' @return named vector of parameter values
coef.dream <- function(object,last.prop=.5,use.thinned=FALSE,
method=c("maxLik","mean","median"),...){
+
stopifnot(last.prop>0)
+
+ if (use.thinned & is.null(object$Reduced.Seq)) {
+ warning("Attempted to use.thinned when no thinned chains available: setting use.thinned=FALSE")
+ use.thinned <- FALSE
+ }
+
if (use.thinned) sss <- object$Reduced.Seq
else sss <- object$Sequences
+ stopifnot(!is.null(sss))
+
if (class(method)!="function") {
method <- switch(
match.arg(method),
@@ -21,7 +30,7 @@
if (last.prop==1) return(method(sss))
else {
- ss <- window(object$Sequences, start = end(sss)*(1-last.prop) + 1)
+ ss <- window(sss, start = end(sss)*(1-last.prop) + 1)
return(method(ss))
}
}##coef.dream
Modified: pkg/R/predict.dream.R
===================================================================
--- pkg/R/predict.dream.R 2010-03-05 03:49:04 UTC (rev 19)
+++ pkg/R/predict.dream.R 2010-03-12 06:40:09 UTC (rev 20)
@@ -30,12 +30,13 @@
##' using various methods of summarising the posterior parameter and output distributions
##' @param object dream object
##' @param newdata. new FUN.pars list. If NULL, use object's.
+##' @param newFUN. a new function to run with same arguments as the original FUN
##' @param method CI or a \code{\link{method}} of coef
##' @param level. Requested two-sided level of confidence. For CI method.
##' @param last.prop Proportion of MCMC chains to keep
##' @param use.thinned Whether to use existing thinned chains
-##' @return data frame with each column corresponding to a returned vector
-predict.dream <- function(object,newdata=NULL,
+##' @return whatever FUN returns (either numeric, ts or list). For CI, either a matrix with upper and lower bound or list of matrices.
+predict.dream <- function(object,newdata=NULL,newFUN=NULL,
method="maxLik",level=0.99,
last.prop=0.5,use.thinned=TRUE,...
){
@@ -44,16 +45,26 @@
stopifnot(is.null(newdata) || is.list(newdata))
stopifnot(!"CI" %in% method || !is.null(level))
stopifnot(last.prop>0)
+
+ if (use.thinned & is.null(object$Reduced.Seq)) {
+ warning("Attempted to use.thinned when no thinned chains available: setting use.thinned=FALSE")
+ use.thinned <- FALSE
+ }
+
###
## Fetch function and parameters from dream object
+
+ if (!is.null(newFUN) & !identical(formals(eval(object$call$FUN)),formals(newFUN))) stop("FUN used in dream and newFUN must have same parameters")
+ if (is.null(newFUN)) newFUN <- eval(object$call$FUN)
if (is.null(newdata)) newdata <- eval(object$call$FUN.pars)
- par.name <- names(formals(eval(object$call$FUN)))[1]
+ par.name <- names(formals(newFUN))[1]
+
wrap <- function(p) {
newdata[[par.name]] <- p
- as.numeric(do.call(eval(object$call$FUN),newdata))
+ do.call(newFUN,newdata)
}
###
@@ -67,7 +78,17 @@
ff <- apply(as.matrix(sss),1,wrap)
- return(t(apply(ff,1,quantile,c((1-level)/2,1-(1-level)/2))))
+ if (inherits(ff,"matrix")) return(t(apply(ff,1,quantile,c((1-level)/2,1-(1-level)/2))))
+ else if (inherits(ff,"list")) {
+ ## Calculate CI for each series separately
+ ## list is of format list[[run.number]][[series.number]]=numeric
+ return(
+ lapply(1:length(ff[[1]]),function(s){
+ ff.s <- sapply(ff,function(x) x[[s]])
+ t(apply(ff.s,1,quantile,c((1-level)/2,1-(1-level)/2)))
+ })
+ )
+ } else stop("Unexpected output from application of newFUN to matrix of parameters. newFUN should return a numeric or list of numerics")
} else {
return(wrap(coef(object,method=method,
Modified: pkg/demo/FME.nonlinear.model.R
===================================================================
--- pkg/demo/FME.nonlinear.model.R 2010-03-05 03:49:04 UTC (rev 19)
+++ pkg/demo/FME.nonlinear.model.R 2010-03-12 06:40:09 UTC (rev 20)
@@ -120,7 +120,6 @@
lines(Obs$x,predict(dd,method="median"),col="orange")
plotCIs(Obs$x,predict(dd,method="CI"),col="black")
-
##Obs2
plotFME()
lines(Obs2$x,predict(dd,list(x=Obs2$x)),col="blue")
Modified: pkg/man/predict.dream.Rd
===================================================================
--- pkg/man/predict.dream.Rd 2010-03-05 03:49:04 UTC (rev 19)
+++ pkg/man/predict.dream.Rd 2010-03-12 06:40:09 UTC (rev 20)
@@ -2,15 +2,20 @@
\alias{predict.dream}
\title{Predict values using dream object}
\usage{
-\S3method{predict}{dream}(object, newdata=NULL,method="maxLik", level=0.99, last.prop=0.5, use.thinned=TRUE,\dots)
+\S3method{predict}{dream}(object, newdata=NULL,newFUN=NULL,method="maxLik", level=0.99, last.prop=0.5, use.thinned=TRUE,\dots)
}
\description{
Predict values using function calibrated by dream, optionally with new data,
using various methods of summarising the posterior parameter and output distributions}
-\value{data frame with each column corresponding to a returned vector}
+\value{whatever newFUN returns (either numeric, ts or list). For CI,
+ either a matrix with upper and lower bound or list of matrices.
+
+ N.B. newFUN can return a list, though FUN for use in dream must return
+a numeric}
\arguments{
\item{object}{dream object}
\item{newdata}{new FUN.pars list. If NULL, use object's.}
+ \item{newFUN}{A new function run, for the case where FUN returned a posterior density, but the model result is now required. newFUN must have exactly the same parameters as FUN. If NULL, don't use a different function.}
\item{method}{CI or a \code{method} of \code{\link{coef.dream}}}
\item{level}{Requested two-sided level of confidence. For CI method.}
\item{last.prop}{Proportion of MCMC chains to keep}
More information about the Dream-commits
mailing list