[Dream-commits] r36 - pkg pkg/R pkg/demo pkg/inst www www/matlab_test
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri May 25 02:01:57 CEST 2012
Author: josephguillaume
Date: 2012-05-25 02:01:56 +0200 (Fri, 25 May 2012)
New Revision: 36
Added:
www/matlab_test/
www/matlab_test/example1.R
www/matlab_test/example1a.mat
www/matlab_test/example1b.mat
www/matlab_test/example2.R
www/matlab_test/example2a.mat
Removed:
pkg/inst/extdata/
Modified:
pkg/DESCRIPTION
pkg/LICENSE
pkg/NEWS
pkg/R/GenCR.R
pkg/demo/FME.nonlinear.model.R
pkg/demo/example1.R
pkg/demo/example2.R
pkg/inst/CITATION
www/index.php
Log:
- Moved matlab files from extdata to website
- Added bare-bones website
- Fixed GenCR R CMD check warning
- Added authors at R field to DESCRIPTION
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2012-04-16 04:04:07 UTC (rev 35)
+++ pkg/DESCRIPTION 2012-05-25 00:01:56 UTC (rev 36)
@@ -1,9 +1,10 @@
Package: dream
-Version: 0.4-1
+Version: 0.4-2
Date: 2012-04-16
Title: DiffeRential Evolution Adaptive Metropolis
-Author: Joseph Guillaume <joseph.guillaume at anu.edu.au>, Felix Andrews <felix at nfrac.org>
-Maintainer: Joseph Guillaume <joseph.guillaume at anu.edu.au>
+Authors at R: c(person("Jasper","Vrugt",role="aut",comment="Matlab original"),
+ person("Joseph","Guillaume",role=c("aut","cre"),comment="R implementation and added functionality",email="joseph.guillaume at anu.edu.au"),
+ person("Felix","Andrews",role="aut",comment="R implementation and added functionality",email="felix at nfrac.org"))
Description: Efficient global MCMC even in high-dimensional spaces.
Based on the original MATLAB code written by Jasper Vrugt.
Methods for calibration and prediction using the DREAM algorithm
Modified: pkg/LICENSE
===================================================================
--- pkg/LICENSE 2012-04-16 04:04:07 UTC (rev 35)
+++ pkg/LICENSE 2012-05-25 00:01:56 UTC (rev 36)
@@ -31,7 +31,3 @@
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-
-MATLAB code written by Jasper A. Vrugt, Center for NonLinear Studies (CNLS)
-
-Written by Jasper A. Vrugt: vrugt at lanl.gov
Modified: pkg/NEWS
===================================================================
--- pkg/NEWS 2012-04-16 04:04:07 UTC (rev 35)
+++ pkg/NEWS 2012-05-25 00:01:56 UTC (rev 36)
@@ -1,5 +1,6 @@
-Changes in Version 0.4-1
+Changes in Version 0.4-2
+ o Remove matlab files to reduce size of package
o turn the calculation of Cb and Wb off if a likelihood function is provided.
o cleaned up demo parallelisation_chain_id, added example of extracting parameter values and likelihood function
o expanded 'see also' section in dream man page
Modified: pkg/R/GenCR.R
===================================================================
--- pkg/R/GenCR.R 2012-04-16 04:04:07 UTC (rev 35)
+++ pkg/R/GenCR.R 2012-05-25 00:01:56 UTC (rev 36)
@@ -15,7 +15,7 @@
## How many candidate points for each crossover value?
## TODO: verify result matches matlab
## MATLAB: [L] = multrnd(MCMCPar.seq * MCMCPar.steps,pCR);
- L <- as.numeric(rmultinom(1,size=control$nseq*control$steps,p=pCR))
+ L <- as.numeric(rmultinom(1,size=control$nseq*control$steps,prob=pCR))
L2 <- c(0,cumsum(L))
## Then select which candidate points are selected with what CR
Modified: pkg/demo/FME.nonlinear.model.R
===================================================================
--- pkg/demo/FME.nonlinear.model.R 2012-04-16 04:04:07 UTC (rev 35)
+++ pkg/demo/FME.nonlinear.model.R 2012-05-25 00:01:56 UTC (rev 36)
@@ -1,7 +1,5 @@
-
## Example from FME vignette, p21, section 7. Soetaert & Laine
## Nonlinear model parameter estimation
-## TODO: document more clearly, better outputs
##
## http://r-forge.r-project.org/plugins/scmsvn/viewcvs.php/pkg/FME/inst/doc/FMEmcmc.Rnw?rev=96&root=fme&view=markup
Modified: pkg/demo/example1.R
===================================================================
--- pkg/demo/example1.R 2012-04-16 04:04:07 UTC (rev 35)
+++ pkg/demo/example1.R 2012-05-25 00:01:56 UTC (rev 36)
@@ -11,9 +11,6 @@
S <- -0.5 * (x %*% imat) %*% as.matrix(x)
return(S)
}
-## Output manually verified against matlab version
-##Banshp(1:10,bpar,imat) ## Banshp(1:10,Extra)
-##Banshp(rep(1,10),bpar,imat) ## Banshp(ones(1,10),Extra)
control <- list(
ndim=10,
@@ -60,41 +57,12 @@
summary(dd)
-## Show bananity
+## Show banana
library(lattice)
ddm <- as.matrix(window(dd))
plot(ddm[,1],ddm[,2])
splom(ddm)
-## Compare to two matlab runs
-fn.example1a <- system.file("extdata/example1a.mat",package="dream")
-fn.example1b <- system.file("extdata/example1b.mat",package="dream")
-
-for (fn.example in c(fn.example1a,fn.example1b)){
- compareToMatlab(fn.example,dd)
- mat <- readMat(fn.example)
- all.equal(muX,as.numeric(mat$Extra[,,1]$muX))
- all.equal(qcov,mat$Extra[,,1]$qcov)
- all.equal(imat,mat$Extra[,,1]$imat)
- all.equal(bpar,as.numeric(mat$Extra[,,1]$bpar))
-}
-
+## Results were compared to two matlab runs
## While banana doesn't appear to match, it appears this is because it is a difficult problem
## Separate matlab results do not match either
-
-## Compare matlab runs
-matb <- getMatlabSeq(readMat(fn.example1b))
-matb <- as.matrix(window(matb,start=1+(end(matb)-1)*(1-0.5)))
-
-mata <- getMatlabSeq(readMat(fn.example1a))
-mata <- as.matrix(window(mata,start=1+(end(mata)-1)*(1-0.5)))
-
-pvals <- sapply(1:ncol(mata), function(i) ks.test(mata[,i], matb[, i])$p.value)
-round(pvals,2)
-
-plotMCMCQQ(matb,mata)
-
-## Results from matlab version
-plot(mata[,1],mata[,2])
-
-splom(mata)
Modified: pkg/demo/example2.R
===================================================================
--- pkg/demo/example2.R 2012-04-16 04:04:07 UTC (rev 35)
+++ pkg/demo/example2.R 2012-05-25 00:01:56 UTC (rev 36)
@@ -1,5 +1,5 @@
## Example 2
-##n-dimensional Gaussian distribution
+## n-dimensional Gaussian distribution
## From Vrugt's matlab code
library(dream)
@@ -33,22 +33,13 @@
}
}
invC <- solve(C) ## checked against matlab for ndim=3 and 7
-##The following lines are defined in the matlab code but not used
-##qcov <- C
-##muX <- rep(0,ndim)
-##Extra.mu = zeros(1,MCMCPar.n); ## Define center of target
normalfunc <- function(x){
- ##p = mvnpdf(x,Extra.mu,Extra.qcov); ##commented out in matlab code
as.numeric(-0.5* x %*% invC %*% matrix(x))
}
-## matches output from matlab
-##normalfunc(replicate(ndim,-3)) ##normalfunc(-3*ones(1,MCMCPar.n),Extra)
-##normalfunc(1:ndim) ## normalfunc(1:MCMCPar.n,Extra)
ddd <- dream(normalfunc,
func.type="logposterior.density",
- ##INIT=LHSInit,
pars=pars,
control=control
)
@@ -74,8 +65,5 @@
checkNormality(window(simulate(ddd,nsim=1e4)))
-##############################################
-
## Output matches that from matlab
-compareToMatlab(system.file("extdata/example2a.mat",package="dream"),ddd)
Modified: pkg/inst/CITATION
===================================================================
--- pkg/inst/CITATION 2012-04-16 04:04:07 UTC (rev 35)
+++ pkg/inst/CITATION 2012-05-25 00:01:56 UTC (rev 36)
@@ -4,7 +4,7 @@
title = paste("Accelerating Markov chain Monte Carlo",
"simulation by differential evolution with",
"self-adaptive randomized subspace sampling."),
- author = personList(as.person("J A Vrugt"),
+ author = c(as.person("J A Vrugt"),
as.person("C J F ter Braak"),
as.person("C G H Diks"),
as.person("B A Robinson"),
@@ -32,7 +32,7 @@
citEntry(header = "To cite the 'dream' package use:",
entry = "Manual",
title = "dream: DiffeRential Evolution Adaptive Metropolis",
- author = personList(as.person("Joseph Guillaume"),
+ author = c(as.person("Joseph Guillaume"),
as.person("Felix Andrews")),
year = year,
note = vers,
Modified: www/index.php
===================================================================
--- www/index.php 2012-04-16 04:04:07 UTC (rev 35)
+++ www/index.php 2012-05-25 00:01:56 UTC (rev 36)
@@ -43,9 +43,62 @@
<!-- end of project description -->
-<p> No content added. </p>
-
<p> The <strong>project summary page</strong> you can find <a href="http://<?php echo $domain; ?>/projects/<?php echo $group_name; ?>/"><strong>here</strong></a>. </p>
+For information on how to use dream, please run in R:<ul>
+<li><i>help("dreamCalibrate")</i> - to calibrate a function using dream</li>
+<li><i>help("dream")</i> - for low-level interface </li>
+<li><i>demo(example1)</i> - Fitting a banana shaped distribution</li>
+<li><i>demo(example2)</i> - Fitting an n-dimensional Gaussian distribution</li>
+<li><i>demo(FME.nonlinear.model)</i> - Calibrating the non-linear model shown
+ in the <a href="http://fme.r-forge.r-project.org/">FME package</a> vignette</li>
+<li><i>demo(FME.nonlinear.model_parallelisation)</i> - Example of parallelisation using the <a href="http://cran.r-project.org/web/packages/snow/index.html">SNOW package</a></li>
+<li><i>demo(parallelisation_chain_id)</i> - Example of parallelisation when DREAM calls an external model using batch files in separate folders</li>
+</ul>
+<p>
+To cite the DREAM algorithm please use:
+<ul>
+<li>
+ Vrugt, J. A., ter Braak, C. J. F., Diks, C. G. H., Robinson, B. A.,
+ Hyman, J. M., Higdon, D., 2009. Accelerating Markov chain Monte Carlo
+ simulation by differential evolution with self-adaptive randomized
+ subspace sampling. <b>International Journal of Nonlinear Sciences
+ and Numerical Simulation</b> 10 (3), 273-290. DOI: <a href="http://dx.doi.org/10.1515/IJNSNS.2009.10.3.273">10.1515/IJNSNS.2009.10.3.273</a>
+</li></ul>
+<p>
+To cite the dream package, please use:<br><ul><li>
+ Joseph Guillaume and Felix Andrews (2012). dream: DiffeRential
+ Evolution Adaptive Metropolis. R package version 0.4-2. URL
+ <a href="http://CRAN.R-project.org/package=dream">http://CRAN.R-project.org/package=dream</a>
+</li></ul>
+
+For additional information on the algorithm also see:<ul>
+<li>
+ Vrugt, J. A., ter Braak, C. J. F., Gupta, H. V., Robinson, B. A.,
+ 2009. Equifinality of formal (DREAM) and informal (GLUE) Bayesian
+ approaches in hydrologic modeling?
+ <b>Stochastic Environmental Research and Risk Assessment</b> 23 (7), 1011--1026.
+ DOI: <a href="http://dx.doi.org/10.1007/s00477-008-0274-y">10.1007/s00477-008-0274-y</a>
+</li>
+</ul>
+
+<p>
+This implementation of DREAM has been tested against the original Matlab implementation. See <a href="./matlab_test/example1.R">example1.R</a> and <a href="./matlab_test/example2.R">example2.R</a>
+<p>
+Please note that the dream_zs and dream_d algorithms may be superior in your circumstances. These are not implemented in this package. Please read the following references for details:
+<ul>
+<li>
+Vrugt, J. A. and Ter Braak, C. J. F. (2011) DREAM(D): an adaptive Markov Chain Monte Carlo simulation algorithm to solve discrete, noncontinuous, and combinatorial posterior parameter estimation problems, <b>Hydrol. Earth Syst. Sci.</b>, 15, 3701-3713, DOI: <a href="http://dx.doi.org/10.5194/hess-15-3701-2011">10.5194/hess-15-3701-2011</a>
+</li>
+<li>
+ter Braak, C. and J. Vrugt (2008). Differential Evolution Markov Chain
+with snooker updater and fewer chains. <b>Statistics and Computing</b> 18(4): 435-446 DOI: <a href="http://dx.doi.org/10.1007/s11222-008-9104-9">10.1007/s11222-008-9104-9</a>
+</li>
+<li>
+Laloy,E., and J.A. Vrugt. 2012. High-dimensional posterior exploration
+of hydrologic models using multiple-try DREAM(ZS) and high-performance
+computing. <b>Water Resources Research</b>, 48, W0156. DOI <a href="http://dx.doi.org/10.1029/2011WR010608">10.1029/2011WR010608</a>
+</li></ul>
+
</body>
</html>
Added: www/matlab_test/example1.R
===================================================================
--- www/matlab_test/example1.R (rev 0)
+++ www/matlab_test/example1.R 2012-05-25 00:01:56 UTC (rev 36)
@@ -0,0 +1,100 @@
+library(dream)
+
+## n-dimensional banana shaped gaussian distribution
+## Hyperbolic shaped posterior probability distribution
+## @param x vector of length ndim
+## @param bpar banana-ness. scalar.
+## @param imat matrix ndim x ndim
+## @return SSR. scalar
+Banshp <- function(x,bpar,imat){
+ x[2] <- x[2]+ bpar * x[1]^2 - 100*bpar
+ S <- -0.5 * (x %*% imat) %*% as.matrix(x)
+ return(S)
+}
+## Output manually verified against matlab version
+##Banshp(1:10,bpar,imat) ## Banshp(1:10,Extra)
+##Banshp(rep(1,10),bpar,imat) ## Banshp(ones(1,10),Extra)
+
+control <- list(
+ ndim=10,
+ DEpairs=3,
+ gamma=0,
+ nCR=3,
+ ndraw=1e5,
+ steps=10,
+ eps=5e-2,
+ outlierTest='IQR_test',
+ pCR.Update=TRUE,
+ thin.t=10,
+ boundHandling='none',
+ burnin.length=Inf, ##for compatibility with matlab code
+ REPORT=1e4 ##reduce frequency of progress reports
+ )
+
+
+pars <- replicate(control$ndim,c(-Inf,Inf),simplify=F)
+cmat <- diag(control$ndim)
+cmat[1,1] <- 100
+bpar <- 0.1
+imat <- solve(cmat)
+muX=rep(0,control$ndim)
+qcov=diag(control$ndim)*5
+
+set.seed(11)
+
+dd <- dream(
+ FUN=Banshp, func.type="logposterior.density",
+ pars = pars,
+ FUN.pars=list(
+ imat=imat,
+ bpar=bpar),
+ INIT = CovInit,
+ INIT.pars=list(
+ muX=muX,
+ qcov=qcov,
+ bound.handling=control$boundHandling
+ ),
+ control = control
+ )
+
+summary(dd)
+
+
+## Show bananity
+library(lattice)
+ddm <- as.matrix(window(dd))
+plot(ddm[,1],ddm[,2])
+splom(ddm)
+
+## Compare to two matlab runs
+fn.example1a <- "http://dream.r-forge.r-project.org/matlab_test/example1a.mat"
+fn.example1b <- "http://dream.r-forge.r-project.org/matlab_test/example1b.mat"
+
+for (fn.example in c(fn.example1a,fn.example1b)){
+ compareToMatlab(fn.example,dd)
+ mat <- readMat(fn.example)
+ all.equal(muX,as.numeric(mat$Extra[,,1]$muX))
+ all.equal(qcov,mat$Extra[,,1]$qcov)
+ all.equal(imat,mat$Extra[,,1]$imat)
+ all.equal(bpar,as.numeric(mat$Extra[,,1]$bpar))
+}
+
+## While banana doesn't appear to match, it appears this is because it is a difficult problem
+## Separate matlab results do not match either
+
+## Compare matlab runs
+matb <- getMatlabSeq(readMat(fn.example1b))
+matb <- as.matrix(window(matb,start=1+(end(matb)-1)*(1-0.5)))
+
+mata <- getMatlabSeq(readMat(fn.example1a))
+mata <- as.matrix(window(mata,start=1+(end(mata)-1)*(1-0.5)))
+
+pvals <- sapply(1:ncol(mata), function(i) ks.test(mata[,i], matb[, i])$p.value)
+round(pvals,2)
+
+plotMCMCQQ(matb,mata)
+
+## Results from matlab version
+plot(mata[,1],mata[,2])
+
+splom(mata)
Property changes on: www/matlab_test/example1.R
___________________________________________________________________
Added: svn:executable
+ *
Added: www/matlab_test/example1a.mat
===================================================================
(Binary files differ)
Property changes on: www/matlab_test/example1a.mat
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ application/octet-stream
Added: www/matlab_test/example1b.mat
===================================================================
(Binary files differ)
Property changes on: www/matlab_test/example1b.mat
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ application/octet-stream
Added: www/matlab_test/example2.R
===================================================================
--- www/matlab_test/example2.R (rev 0)
+++ www/matlab_test/example2.R 2012-05-25 00:01:56 UTC (rev 36)
@@ -0,0 +1,81 @@
+## Example 2
+##n-dimensional Gaussian distribution
+## From Vrugt's matlab code
+
+library(dream)
+## Original version used ndim=100. The current R version has memory problems with such high dimensions
+ndim <- 9
+## Commented lines are default settings
+control <- list(
+ thin.t=10,
+ ##pCR.Update=T,
+ ##ndim=100,
+ nseq=ndim,
+ ndraw=1e5, ##The number of runs was reduced from 1e6 for speed
+ ##nCR=3,
+ ##gamma=0,
+ DEpairs=3,
+ ##steps=10,
+ ##eps=5e-2,
+ ##outlierTest='IQR_test',
+ boundHandling='none',
+ burnin.length=Inf, ##compatibility with matlab version, remove outliers all the time
+ REPORT=1e4 ## Reduce frequency of progress reporting
+ )
+
+pars <- replicate(ndim,c(-5,15),simplify=F)
+
+A <- 0.5*diag(ndim)+0.5
+C <- matrix(NA,nrow=ndim,ncol=ndim)
+for (i in 1:ndim){
+ for (j in 1:ndim){
+ C[i,j] <- A[i,j]*sqrt(i*j)
+ }
+}
+invC <- solve(C) ## checked against matlab for ndim=3 and 7
+##The following lines are defined in the matlab code but not used
+##qcov <- C
+##muX <- rep(0,ndim)
+##Extra.mu = zeros(1,MCMCPar.n); ## Define center of target
+
+normalfunc <- function(x){
+ ##p = mvnpdf(x,Extra.mu,Extra.qcov); ##commented out in matlab code
+ as.numeric(-0.5* x %*% invC %*% matrix(x))
+}
+## matches output from matlab
+##normalfunc(replicate(ndim,-3)) ##normalfunc(-3*ones(1,MCMCPar.n),Extra)
+##normalfunc(1:ndim) ## normalfunc(1:MCMCPar.n,Extra)
+
+ddd <- dream(normalfunc,
+ func.type="logposterior.density",
+ ##INIT=LHSInit,
+ pars=pars,
+ control=control
+ )
+
+################################################
+
+library(reshape)
+library(latticeExtra)
+checkNormality <- function(seqs){
+ print(round(sapply(apply(as.matrix(seqs),2,shapiro.test),
+ function(x) x$p.value),2))
+ xxw.long <- melt(as.matrix(seqs))
+ qqmath(~value|X2,xxw.long,as.table=T)+layer(panel.qqmathline(x))
+}
+
+
+summary(ddd)
+
+## Using fraction of sequences
+checkNormality(window(ddd,frac=0.45))
+
+## Using extension of sequences
+checkNormality(window(simulate(ddd,nsim=1e4)))
+
+
+##############################################
+
+## Output matches that from matlab
+compareToMatlab("http://dream.r-forge.r-project.org/matlab_test/example2a",ddd)
+
Property changes on: www/matlab_test/example2.R
___________________________________________________________________
Added: svn:executable
+ *
Added: www/matlab_test/example2a.mat
===================================================================
(Binary files differ)
Property changes on: www/matlab_test/example2a.mat
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ application/octet-stream
More information about the Dream-commits
mailing list