[Dream-commits] r18 - in pkg: . R demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Mar 4 00:12:01 CET 2010
Author: josephguillaume
Date: 2010-03-04 00:12:00 +0100 (Thu, 04 Mar 2010)
New Revision: 18
Modified:
pkg/DESCRIPTION
pkg/R/CompDensity.R
pkg/R/dream.R
pkg/R/offde.R
pkg/demo/FME.nonlinear.model.R
pkg/man/dream.Rd
Log:
Multicore support + associated help changes
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-03-03 07:48:11 UTC (rev 17)
+++ pkg/DESCRIPTION 2010-03-03 23:12:00 UTC (rev 18)
@@ -7,6 +7,7 @@
Description: Efficient global MCMC even in high-dimensional spaces.
Based on the original MATLAB code written by Jasper Vrugt.
Depends: coda
+Suggests: multicore
Imports: stats, utils
License: file LICENSE
URL: http://dream.r-forge.r-project.org/
Modified: pkg/R/CompDensity.R
===================================================================
--- pkg/R/CompDensity.R 2010-03-03 07:48:11 UTC (rev 17)
+++ pkg/R/CompDensity.R 2010-03-03 23:12:00 UTC (rev 18)
@@ -35,8 +35,7 @@
## SSR scalar
## temp. list of length nseq with elements of length 2: p and logp
- ## Ready for multicore
- temp <- lapply(1:nrow(pars),function (ii){
+ do.calc <- function (ii){
## Call model to generate simulated data
## TODO: correct use of optional pars?
@@ -87,8 +86,11 @@
logp <- -0.5*SSR
}) ##switch
c(p,logp)
- }) ## lapply
-
+ }
+
+ if (control$use.multicore) temp <- mclapply(1:nrow(pars),do.calc,mc.preschedule=FALSE)
+ else temp <- lapply(1:nrow(pars),do.calc)
+
p <- sapply(temp,function(x) x[1])
logp <- sapply(temp,function(x) x[2])
Modified: pkg/R/dream.R
===================================================================
--- pkg/R/dream.R 2010-03-03 07:48:11 UTC (rev 17)
+++ pkg/R/dream.R 2010-03-03 23:12:00 UTC (rev 18)
@@ -16,8 +16,9 @@
Rthres=1.01, ## R value at which to stop. Vrugt suggests 1.2
### Thinning
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
+### Efficiency improvements
+ REPORT = 1000, ## approximate number of function evaluations between reports. >0. 0=none
+ use.multicore=TRUE,
### Parameters with auto-set values
ndim=NA, ## number of parameters (automatically set from length of pars)
DEpairs = NA, ## Number of DEpairs. defaults to max val floor((nseq-1)/2)
@@ -133,6 +134,8 @@
control$REPORT <- (control$REPORT%/%control$nseq) * control$nseq
if (control$burnin.length<1) control$burnin.length <- control$burnin.length*control$ndraw
+
+ if (control$use.multicore) control$use.multicore <- require(multicore)
## Check validity of settings
if (control$DEpairs==0) stop("control$DEpairs set to 0. Increase nseq?")
@@ -140,6 +143,7 @@
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)
+
############################
## Initialize variables
Modified: pkg/R/offde.R
===================================================================
--- pkg/R/offde.R 2010-03-03 07:48:11 UTC (rev 17)
+++ pkg/R/offde.R 2010-03-03 23:12:00 UTC (rev 18)
@@ -19,7 +19,6 @@
nseq <- control$nseq
ndim <- control$ndim
-
## dimensions
## eps. matrix nseq x ndim
## DEversion. vector of length nseq. range [0,DEpairs]
@@ -57,6 +56,9 @@
delta.x<-matrix(0,nseq,ndim)
## Each chain evolves using information from other chains to create offspring
+ ## TODO: can this loop be parallelised?
+ ## is an expensive part of dream. 36.6% of time
+ ## but depends on state of other loops
for (qq in 1:nseq){
## Define ii and remove current member as an option
Modified: pkg/demo/FME.nonlinear.model.R
===================================================================
--- pkg/demo/FME.nonlinear.model.R 2010-03-03 07:48:11 UTC (rev 17)
+++ pkg/demo/FME.nonlinear.model.R 2010-03-03 23:12:00 UTC (rev 18)
@@ -51,8 +51,8 @@
control <- list(
nseq=4
-## REPORT=0
-## ndraw=1000
+ ## REPORT=0
+ ## ndraw=1000
## Rthres=1+1e-3
)
Modified: pkg/man/dream.Rd
===================================================================
--- pkg/man/dream.Rd 2010-03-03 07:48:11 UTC (rev 17)
+++ pkg/man/dream.Rd 2010-03-03 23:12:00 UTC (rev 18)
@@ -32,34 +32,48 @@
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
+ \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
+ \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
+ \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
+ \cr
+
+ \code{burnin.length} \tab 0.1 \tab If \code{<1} proportion of
+ function evaluations to use as burn-in period. If \code{>1} number
+ of function evaluations to use.
+ In the burnin-period, adaptive tuning of crossover values is
+ performed if \code{pCR.Update=TRUE} and the method specified by
+ \code{outlierTest} is used to remove outliers.
+ If outliers are detected outside the burn-in period, dream emits a
+ warning and returns to the burn-in period for another \code{burnin.length} iterations.
+ If \code{=0} there is no burn-in period, and all outlier warnings
+ are suppressed. \cr
+
+ \code{thin.t} \tab NA \tab MCMC chain 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
+ Frequent calculation is likely to slow performance. The summary
+ function can use \code{\link{gelman.diag}} to calculate convergence statistic after completion
+ \cr
+ \code{use.multicore} \tab \code{TRUE} \tab Whether to use multicore
+ package if available. Parallelises function evaluations. \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
+ \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
+ \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
+ \cr
}
Execution terminates when Gelman-Rubin statistic < \code{control$Rthres}, or
\code{control$ndraw} or \code{control$maxtime} are reached
@@ -69,10 +83,16 @@
Default settings may not be sufficient for model to run. For few
parameters, nseq may need to be increased.
+
+ If dream frequently warns that it is reentering the burn-period,
+ \code{control$burnin.length} may need to be increased.
- There are methods for print, summary, plot, coef
+ There are S3 methods for print, summary, plot,
+ \code{\link[=coef.dream]{coef}},
+ \code{\link[=predict.dream]{predict}},
+ simulate
- coef uses \code{\link{maxLikPars}}, a naive approach to extract the most likely parameter value.
+ coef uses \code{\link{maxLikPars}}, a naive approach to extracting the most likely parameter value.
}
\value{
@@ -84,7 +104,10 @@
\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)}
+ \item{CR}{Probability of crossover at each step. matrix n.elem/steps x
+ 1+length(pCR)}
+ \item{in.burnin}{Boolean showing whether dream was in the burn-in
+ period at termination (see description above)}
}
\seealso{
More information about the Dream-commits
mailing list