[Dream-commits] r13 - in pkg: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Feb 26 07:59:11 CET 2010
Author: josephguillaume
Date: 2010-02-26 07:59:10 +0100 (Fri, 26 Feb 2010)
New Revision: 13
Added:
pkg/man/
pkg/man/CovInit.Rd
pkg/man/LHSInit.Rd
pkg/man/dream.Rd
pkg/man/maxLikPars.Rd
pkg/man/possibility.envelope.Rd
Modified:
pkg/NAMESPACE
pkg/R/coef.dream.R
pkg/R/dream.R
pkg/R/plot.dream.R
pkg/R/print.dream.R
pkg/R/summary.dream.R
Log:
Added documentation. Passing R CMD check.
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2010-02-26 04:31:57 UTC (rev 12)
+++ pkg/NAMESPACE 2010-02-26 06:59:10 UTC (rev 13)
@@ -4,7 +4,8 @@
export(CovInit)
export(LHSInit)
export(maxLikPars)
+export(possibility.envelope)
S3method(summary,dream)
S3method(coef,dream)
S3method(plot,dream)
-S3method(print,dream)
\ No newline at end of file
+S3method(print,dream)
Modified: pkg/R/coef.dream.R
===================================================================
--- pkg/R/coef.dream.R 2010-02-26 04:31:57 UTC (rev 12)
+++ pkg/R/coef.dream.R 2010-02-26 06:59:10 UTC (rev 13)
@@ -3,14 +3,14 @@
##' @last.prop proportion of total sequence to use (0,1]
##' if 1, use whole sequence
##' @return named vector of parameter values
-coef.dream <- function(dream.obj,last.prop=.5,use.thinned=FALSE){
+coef.dream <- function(object,last.prop=.5,use.thinned=FALSE,...){
stopifnot(last.prop>0)
- if (use.thinned) sss <- dream.obj$Reduced.Seq
- else sss <- dream.obj$Sequences
+ if (use.thinned) sss <- object$Reduced.Seq
+ else sss <- object$Sequences
if (last.prop==1) return(maxLikPars(sss))
else {
- ss <- window(dream.obj$Sequences, start = end(sss)*(1-last.prop) + 1)
+ ss <- window(object$Sequences, start = end(sss)*(1-last.prop) + 1)
return(maxLikPars(ss))
}
}##coef.dream
Modified: pkg/R/dream.R
===================================================================
--- pkg/R/dream.R 2010-02-26 04:31:57 UTC (rev 12)
+++ pkg/R/dream.R 2010-02-26 06:59:10 UTC (rev 13)
@@ -14,7 +14,6 @@
maxtime = Inf, ## maximum duration of optimization in seconds
Rthres=1.01, ## R value at which to stop. Vrugt suggests 1.2
### Thinning
- thin=FALSE, ## do reduced sample collection
thin.t=NA, ## parameter for reduced sample collection
### Reporting
REPORT = 1000, ## approximate number of function evaluations between reports. >0. 0=none TODO: when trace >= 1
@@ -56,8 +55,7 @@
##' MATLAB function:
##' function [Sequences,Reduced_Seq,X,output,hist_logp] = dream(MCMCPar,ParRange,Measurement,ModelName,Extra,option)
-dream <- function(FUN, func.type,
- pars = list(x = range(0, 1e6)),
+dream <- function(FUN, func.type,pars,
FUN.pars=list(),
INIT = LHSInit,
INIT.pars=list(),
@@ -87,25 +85,26 @@
############################
## Process parameters
+ ## Check validity of parameters
if (is.character(FUN)) FUN <- get(FUN, mode = "function")
stopifnot(is.function(FUN))
stopifnot(is.list(pars))
stopifnot(length(pars) > 0)
+ pars <- lapply(pars, function(x) if (is.list(x)) x else list(x))
+ stopifnot(is.list(control))
+ stopifnot(func.type %in% c("calc.rmse","calc.loglik","calc.weighted.rmse","posterior.density","logposterior.density"))
stopifnot(!is.null(measurement) || func.type %in% c("posterior.density","logposterior.density"))
- stopifnot(func.type!="calc.rmse" || "data" %in% names(measurement))
+ stopifnot(!func.type %in% c("calc.rmse","calc.loglik","calc.weighted.rmse") || "data" %in% names(measurement))
- stopifnot(control$boundHandling %in% c("reflect", "bound", "fold", "none"))
-
+ ## Check INIT and FUN have required extra parameters in INIT.pars & FUN.pars
req.args.init <- names(formals(INIT))
req.args.FUN <- names(formals(FUN))
-
+
if(!all(req.args.init %in% c("pars","nseq",names(INIT.pars)))) stop(paste(c("INIT Missing extra arguments:",req.args.init[!req.args.init %in% c("pars","nseq",names(INIT.pars))]),sep=" "))
if(!all(req.args.FUN[2:length(req.args.FUN)] %in% c(names(FUN.pars)))) stop(paste(c("FUN Missing extra arguments:",req.args.FUN[!req.args.FUN %in% c("x",names(FUN.pars))]),sep=" "))
-
- pars <- lapply(pars, function(x) if (is.list(x)) x else list(x))
+
+ ## Update default settings with supplied settings
- ## update default options with supplied options
- stopifnot(is.list(control))
control <- modifyList(dreamDefaults(), control)
isValid <- names(control) %in% names(dreamDefaults())
if (any(!isValid))
@@ -116,18 +115,22 @@
if (! "sigma" %in% names(measurement)) measurement$sigma <- sd(measurement$data)
if (! "N" %in% names(measurement)) measurement$N <- length(measurement$data)
}
-
- ## determine number of variables to be optimized
+
+ ## Set automatically determined values
control$ndim<-length(pars)
if (is.na(control$nseq)) control$nseq <- control$ndim
if (is.na(control$DEpairs)) control$DEpairs <- floor((control$nseq-1)/2)
+ ## Correct to match nseq
+ control$REPORT <- (control$REPORT%/%control$nseq) * control$nseq
+
+
+ ## Check validity of settings
+ if (control$DEpairs==0) stop("control$DEpairs set to 0. Increase nseq?")
stopifnot(control$DEpairs<=(control$nseq-1)/2) ## Requirement of offde
-
+ stopifnot(control$boundHandling %in% c("reflect", "bound", "fold", "none"))
if (control$boundHandling == 'none') warning("No bound handling in use, parameters may cause errors elsewhere")
-
stopifnot(control$REPORT>=0)
- control$REPORT <- (control$REPORT%/%control$nseq) * control$nseq
############################
## Initialize variables
@@ -136,8 +139,8 @@
NCR <- control$nCR
NSEQ <- control$nseq
- ## for each iteration...
- counter.fun.evals <- NSEQ #? 1
+ ## Counters
+ counter.fun.evals <- NSEQ
counter <- 2
iloc <- 1
counter.outloop <- 2
@@ -209,7 +212,7 @@
## Sequences[1,] <- sapply(pars, mean) ## TODO: include?
## Check whether will save a reduced sample
- if (control$thin){
+ if (!is.na(control$thin.t)){
iloc.2 <- 0
Reduced.Seq <- array(NA,c(floor(n.elem/control$thin.t),NDIM+2,NSEQ))
} else Reduced.Seq <- NULL
@@ -302,7 +305,7 @@
Sequences[iloc,,] <- t(newgen)
## Check reduced sample collection
- if (control$thin && counter.thin == control$thin.t){
+ if (!is.na(control$thin.t) && counter.thin == control$thin.t){
## Update iloc_2 and counter.thin
iloc.2 <- iloc.2+1
counter.thin <- 0
@@ -431,7 +434,7 @@
## Convert sequences to mcmc objects
Sequences <- Sequences[1:iloc,,]
obj$Sequences <- as.mcmc.list(lapply(1:NSEQ,function(i) as.mcmc(Sequences[,1:NDIM,i])))
- if (control$thin){
+ if (!is.na(control$thin)){
Reduced.Seq <- Reduced.Seq[1:iloc.2,,]
obj$Reduced.Seq <- as.mcmc.list(lapply(1:NSEQ,function(i) mcmc(
Reduced.Seq[,1:NDIM,i],
Modified: pkg/R/plot.dream.R
===================================================================
--- pkg/R/plot.dream.R 2010-02-26 04:31:57 UTC (rev 12)
+++ pkg/R/plot.dream.R 2010-02-26 06:59:10 UTC (rev 13)
@@ -4,10 +4,10 @@
##' Uses second half of sequences
-plot.dream <- function(dream.obj,interactive=TRUE){
+plot.dream <- function(x,interactive=TRUE,...){
devAskNewPage(interactive)
- ss <- window(dream.obj$Sequences, start = end(dream.obj$Sequences)/2 + 1)
+ ss <- window(x$Sequences, start = end(x$Sequences)/2 + 1)
## Trace and parameter density
@@ -17,19 +17,19 @@
densityplot(ss)
## Acceptance rate
- plot(table(dd$AR[,2]),main="Distribution of % acceptance rate")
+ plot(table(x$AR[,2]),main="Distribution of % acceptance rate")
##Convergence
try(gelman.plot(ss))
- plot(dd$R.stat[,1],dd$R.stat[,2],type="l",ylim=c(0,2))
- for (i in 2:dd$control$ndim) lines(dd$R.stat[,1],dd$R.stat[,i+1],ylim=c(0,2))
+ plot(x$R.stat[,1],x$R.stat[,2],type="l",ylim=c(0,2))
+ for (i in 2:x$control$ndim) lines(x$R.stat[,1],x$R.stat[,i+1],ylim=c(0,2))
title(main="Evolution of R.stat",sub="Equivalent to gelman.plot")
## Multi-variate density for first chain
- splom(as.data.frame(dd$Sequences[[1]]),
+ splom(as.data.frame(x$Sequences[[1]]),
upper.panel = panel.smoothScatter, nrpoints = 0,
lower.panel = function(x, y, ...) {
panel.grid(-1, -1)
Modified: pkg/R/print.dream.R
===================================================================
--- pkg/R/print.dream.R 2010-02-26 04:31:57 UTC (rev 12)
+++ pkg/R/print.dream.R 2010-02-26 06:59:10 UTC (rev 13)
@@ -1,18 +1,19 @@
-print.dream <- function(dream.obj){
+print.dream <- function(x,...){
cat("
Call:
")
- print(dream.obj$call)
+ print(x$call)
cat("
Control:
")
- for (i in names(dream.obj$control)){
- v <- dream.obj$control[[i]]
+ for (i in names(x$control)){
+ v <- x$control[[i]]
cat(sprintf("%15s: ",i))
cat(v,"\n")
}
- cat("\nExit condition:",dream.obj$EXITMSG,"\n")
+ cat("\nExit condition:",x$EXITMSG,"\n")
+
}
Modified: pkg/R/summary.dream.R
===================================================================
--- pkg/R/summary.dream.R 2010-02-26 04:31:57 UTC (rev 12)
+++ pkg/R/summary.dream.R 2010-02-26 06:59:10 UTC (rev 13)
@@ -1,5 +1,5 @@
-summary.dream <- function(dream.obj){
+summary.dream <- function(object,...){
cat(sprintf("
Exit message: %s
@@ -7,15 +7,15 @@
Time (secs): %.1f
Final R.stats:
",
- dream.obj$EXITMSG,
- dream.obj$fun.evals,
- dream.obj$time
+ object$EXITMSG,
+ object$fun.evals,
+ object$time
))
- R.stat.last <- tail(dream.obj$R.stat,1)
- for (i in 2:ncol(dream.obj$R.stat)){
+ R.stat.last <- tail(object$R.stat,1)
+ for (i in 2:ncol(object$R.stat)){
cat(sprintf("\t%s:\t%f\n",
- colnames(dream.obj$R.stat)[i],
+ colnames(object$R.stat)[i],
R.stat.last[i]
))
} ##for
@@ -23,11 +23,11 @@
cat("
CODA summary for last 50% of MCMC chains:
")
- print(summary(window(dream.obj$Sequences, start = end(dream.obj$Sequences)/2 + 1)))
+ print(summary(window(object$Sequences, start = end(object$Sequences)/2 + 1)))
cat("
Acceptance Rate
")
- summary(dream.obj$AR[,2])
+ summary(object$AR[,2])
} ##summary.dream
Added: pkg/man/CovInit.Rd
===================================================================
--- pkg/man/CovInit.Rd (rev 0)
+++ pkg/man/CovInit.Rd 2010-02-26 06:59:10 UTC (rev 13)
@@ -0,0 +1,16 @@
+\name{CovInit}
+\alias{CovInit}
+\title{Alternative sampling method...}
+\usage{CovInit(pars, nseq, muX, qcov, bound.handling)}
+\description{Alternative initial sampling method
+used in examples, may not be useful in practice}
+\value{
+matrix dim nseq x ndim. parameter wise range [xmin,xmax] assured by handleBounds
+When used as input to dream, muX, qcov and bound.handling must be passed as extra parameters}
+\arguments{
+ \item{pars}{list of parameter vectors}
+ \item{nseq}{Number of chains to initialise. scalar}
+ \item{muX}{vector of length ndim}
+ \item{qcov}{Covariance matrix}
+ \item{bound.handling}{one of: "none","reflect","bound","fold","rand"}
+}
Property changes on: pkg/man/CovInit.Rd
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/man/LHSInit.Rd
===================================================================
--- pkg/man/LHSInit.Rd (rev 0)
+++ pkg/man/LHSInit.Rd 2010-02-26 06:59:10 UTC (rev 13)
@@ -0,0 +1,8 @@
+\name{LHSInit}
+\alias{LHSInit}
+\title{Latin Hypercube sampling...}
+\usage{LHSInit(pars, nseq)}
+\description{Latin Hypercube sampling}
+\value{x matrix nseq x ndim. parameter-wise range [xmin, xmax]}
+\arguments{\item{pars}{list. length ndim - one element per parameter}
+\item{nseq}{scalar = number of sequences}}
Property changes on: pkg/man/LHSInit.Rd
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/man/dream.Rd
===================================================================
--- pkg/man/dream.Rd (rev 0)
+++ pkg/man/dream.Rd 2010-02-26 06:59:10 UTC (rev 13)
@@ -0,0 +1,113 @@
+\name{dream}
+\alias{dream}
+\title{DiffeRential Evolution Adaptive Metropolis (DREAM)}
+\usage{dream(FUN, func.type, pars, FUN.pars=list(), INIT=LHSInit,
+INIT.pars=list(), control=list(), measurement)}
+\description{
+ Efficient global MCMC even in high-dimensional spaces.
+ From J.A. Vrugt, C.J.F. ter Braak et al.
+}
+\arguments{
+ \item{FUN}{model function with first argument a vector of parameter values of length ndim}
+ \item{func.type}{type of value FUN returns.
+ one of: posterior.density, logposterior.density,calc.loglik,calc.rmse,calc.weighted.rmse}
+ \item{pars}{a list of variable ranges. Any names will be propagated to output}
+ \item{FUN.pars}{A named list of extra arguments to FUN}
+ \item{INIT}{A function f(pars,nseq,...) returning an nseq x ndim
+ matrix of initial parameter values}
+ \item{INIT.pars}{A named list of extra arguments to INIT}
+ \item{control}{list of settings for dream. see below}
+ \item{measurement}{Required parameter for func.types starting with calc. list with
+ element data, a vector containing observed dependent variable values.
+ }
+}
+\details{
+ Elements of control are:
+ \tabular{lll}{
+ \emph{Element} \tab \emph{Default} \tab \emph{Description} \cr
+ \code{ndim} \tab \tab Number of parameters. Calculated from
+ parameters \code{pars} \cr
+ \code{nseq} \tab \code{ndim} \tab Number of parallel chains to evolve \cr
+ \code{DEpairs} \tab \code{(nseq-1)/2} \tab Number of pairs to evolve at each
+ generation \cr
+ \code{nCR} \tab 3 \tab Crossover values used to generate proposals (geometric
+ series). scalar \cr
+ \code{gamma} \tab 0 \tab Kurtosis parameter Bayesian Inference Scheme \cr
+ \code{steps} \tab 10 \tab Number of steps in sem \cr
+ \code{eps} \tab 5e-2 \tab Random error for ergodicity \cr
+ \code{outlierTest} \tab 'IQR_test' \tab Test used to detect outlier
+ chains.
+ One off: IQR_test, Grubbs_test,Mahal_test \cr
+ \code{pCR.Update} \tab \code{TRUE} \tab Whether to use adaptive tuning of crossover
+ values \cr
+ \code{boundHandling} \tab 'reflect' \tab Method used to handle parameter
+ values outside of parameter bounds.
+ One of: "reflect", "bound", "fold", "none","rand"
+ \cr
+ \code{thin.t} \tab NA \tab Thinning interval. NA if thinning is not
+ desired \cr
+ \code{REPORT} \tab 1000 \tab Approximate number of function evaluations
+ between calculation and reporting of convergence diagnostics. 0 if
+ no reporting or calculation is desired.
+ Frequent calculation is likely to slow performance.
+ \code{\link{gelman.diag}} can calculate convergence statistic after completion
+ \cr
+ \code{ndraw} \tab 1e5 \tab Maximum number of function evaluations. May
+ terminate before convergence. \cr
+ \code{maxtime} \tab Inf \tab Maximum duration of optimization in seconds. May
+ terminate before convergence.\cr
+ \code{Rthres} \tab 1.01 \tab Value of Gelman & Rubin's convergence
+ diagnostic R value below which the sequences are considered to have
+ converged, and execution is terminated. Vrugt suggests 1.2
+ \cr
+ }
+ Execution terminates when Gelman-Rubin statistic < \code{control$Rthres}, or
+ \code{control$ndraw} or \code{control$maxtime} are reached
+
+ Two options for INIT are provided: latin hypercube sampling
+ \code{\link{LHSInit}} and covariance-based sampling \code{\link{CovInit}}.
+
+ Default settings may not be sufficient for model to run. For few
+ parameters, nseq may need to be increased.
+
+ There are methods for print, summary, plot, coef
+
+ coef uses \code{\link{maxLikPars}}, a naive approach to extract the most likely parameter value.
+
+}
+\value{
+ A list with elements:
+ \item{X}{converged nseq points in parameter space. matrix nseq x ndim}
+ \item{Sequences}{ \code{\link{mcmc.list}}. nseq mcmc elements of ndim variables}
+ \item{Reduced.Seq}{ \code{\link{mcmc.list}}. nseq mcmc elements of ndim
+ variables, if \code{control$thin.t!=NA}}
+ \item{AR}{Acceptance rate for each draw. matrix n.elem x 2}
+ \item{outlier}{Iterations at which chains were replaced. vector of variable length}
+ \item{R.stat}{gelman.diag statistic for each variable at each step. matrix n.elem/steps x 1+ndim}
+ \item{CR}{Probability of crossover at each step. matrix n.elem/steps x 1+length(pCR)}
+}
+
+\seealso{
+ \code{\link{possibility.envelope}}
+
+ Examples in demo folder:
+ \itemize{
+ \item{FME non linear model: }{Calibrating the non-linear model shown in the FME vignette}
+ \item{example1: }{Fitting a banana shaped distribution - the first example
+ in Vrugt matlab code}
+ }
+ }
+
+\references{
+
+ Accelerating Markov Chain Monte Carlo Simulation by Differential Evolution with Self-Adaptive Randomized Subspace Sampling
+Jasper A. Vrugt, C.J.F. ter Braak, C.G.H. Diks, Bruce A. Robinson, James M. Hyman and Dave Higdon
+ International Journal of Nonlinear Sciences and Numerical Simulation
+ (in press 2009)
+
+ Equifinality of formal (DREAM) and informal (GLUE)
+Bayesian approaches in hydrologic modeling?
+Jasper A. Vrugt, Cajo J. F. ter Braak, Hoshin V. Gupta, Bruce A. Robinson
+Stoch Environ Res Risk Assess (2009) 23:1011\2261026
+DOI 10.1007/s00477-008-0274-y
+ }
\ No newline at end of file
Property changes on: pkg/man/dream.Rd
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/man/maxLikPars.Rd
===================================================================
--- pkg/man/maxLikPars.Rd (rev 0)
+++ pkg/man/maxLikPars.Rd 2010-02-26 06:59:10 UTC (rev 13)
@@ -0,0 +1,10 @@
+\name{maxLikPars}
+\alias{maxLikPars}
+\title{Naive maximum density parameter selection...}
+\usage{maxLikPars(ss, ...)}
+\description{Naive maximum density parameter selection}
+\value{named character vector of parameter values
+
+Uses which.max and density function}
+\arguments{\item{ss}{MCMC chains. mcmc or mcmclist object}
+\item{...}{extra parameters to pass to density}}
Property changes on: pkg/man/maxLikPars.Rd
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/man/possibility.envelope.Rd
===================================================================
--- pkg/man/possibility.envelope.Rd (rev 0)
+++ pkg/man/possibility.envelope.Rd 2010-02-26 06:59:10 UTC (rev 13)
@@ -0,0 +1,9 @@
+\name{possibility.envelope}
+\alias{possibility.envelope}
+\title{Obtain prediction confidence intervals for model, around new input values...}
+\usage{possibility.envelope(dd, FUN.pars, ndraw=1000, conf=99)}
+\description{Obtain prediction confidence intervals for model, around new input values}
+\arguments{\item{dd}{dream object}
+\item{FUN.pars}{new par values. if missing, use those from dream object}
+\item{ndraw}{number of new iterations}
+\item{conf}{Percentage two-sided confidence interval}}
Property changes on: pkg/man/possibility.envelope.Rd
___________________________________________________________________
Name: svn:executable
+ *
More information about the Dream-commits
mailing list