[Dream-commits] r26 - in pkg: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed May 12 06:44:34 CEST 2010
Author: felix
Date: 2010-05-12 06:44:34 +0200 (Wed, 12 May 2010)
New Revision: 26
Added:
pkg/R/fitted.dream.R
Modified:
pkg/DESCRIPTION
pkg/LICENSE
pkg/NAMESPACE
pkg/R/CompDensity.R
pkg/R/dream.R
pkg/R/plot.dream.R
pkg/R/predict.dream.R
pkg/R/summary.dream.R
pkg/man/dream.Rd
Log:
rework internal functions a bit; new fitted() method
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-05-11 05:15:27 UTC (rev 25)
+++ pkg/DESCRIPTION 2010-05-12 04:44:34 UTC (rev 26)
@@ -2,7 +2,7 @@
Version: 0.2-1
Date: 2010-05-11
Title: DiffeRential Evolution Adaptive Metropolis
-Author: Jasper Vrugt, CJF ter Braak, et al. (R port by Joseph Guillaume)
+Author: Joseph Guillaume and Felix Andrews, based on MATLAB code by Jasper Vrugt
Maintainer: Joseph Guillaume <joseph.guillaume at anu.edu.au>
Description: Efficient global MCMC even in high-dimensional spaces.
Based on the original MATLAB code written by Jasper Vrugt.
Modified: pkg/LICENSE
===================================================================
--- pkg/LICENSE 2010-05-11 05:15:27 UTC (rev 25)
+++ pkg/LICENSE 2010-05-12 04:44:34 UTC (rev 26)
@@ -1,37 +1,37 @@
-% Copyright (c) 2008, Los Alamos National Security, LLC %
-% All rights reserved. %
-% %
-% Copyright 2008. Los Alamos National Security, LLC. This software was produced under U.S. %
-% Government contract DE-AC52-06NA25396 for Los Alamos National Laboratory (LANL), which is %
-% operated by Los Alamos National Security, LLC for the U.S. Department of Energy. The U.S. %
-% Government has rights to use, reproduce, and distribute this software. %
-% %
-% NEITHER THE GOVERNMENT NOR LOS ALAMOS NATIONAL SECURITY, LLC MAKES A NY WARRANTY, EXPRESS OR %
-% IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. If software is modified to %
-% produce derivative works, such modified software should be clearly marked, so as not to %
-% confuse it with the version available from LANL. %
-% %
-% Additionally, redistribution and use in source and binary forms, with or without %
-% modification, are permitted provided that the following conditions are met: %
-% • Redistributions of source code must retain the above copyright notice, this list of %
-% conditions and the following disclaimer. %
-% • Redistributions in binary form must reproduce the above copyright notice, this list of %
-% conditions and the following disclaimer in the documentation and/or other materials %
-% provided with the distribution. %
-% • Neither the name of Los Alamos National Security, LLC, Los Alamos National Laboratory, LANL %
-% the U.S. Government, nor the names of its contributors may be used to endorse or promote %
-% products derived from this software without specific prior written permission. %
-% %
-% THIS SOFTWARE IS PROVIDED BY LOS ALAMOS NATIONAL SECURITY, LLC AND CONTRIBUTORS "AS IS" AND %
-% ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES %
-% OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL LOS %
-% ALAMOS NATIONAL SECURITY, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, %
-% SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF %
-% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) %
-% 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 %
+Copyright (c) 2008, Los Alamos National Security, LLC
+All rights reserved.
+
+Copyright 2008. Los Alamos National Security, LLC. This software was produced under U.S.
+Government contract DE-AC52-06NA25396 for Los Alamos National Laboratory (LANL), which is
+operated by Los Alamos National Security, LLC for the U.S. Department of Energy. The U.S.
+Government has rights to use, reproduce, and distribute this software.
+
+NEITHER THE GOVERNMENT NOR LOS ALAMOS NATIONAL SECURITY, LLC MAKES A NY WARRANTY, EXPRESS OR
+IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. If software is modified to
+produce derivative works, such modified software should be clearly marked, so as not to
+confuse it with the version available from LANL.
+
+Additionally, redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+* Redistributions of source code must retain the above copyright notice, this list of
+ conditions and the following disclaimer.
+* Redistributions in binary form must reproduce the above copyright notice, this list of
+ conditions and the following disclaimer in the documentation and/or other materials
+ provided with the distribution.
+* Neither the name of Los Alamos National Security, LLC, Los Alamos National Laboratory, LANL
+ the U.S. Government, nor the names of its contributors may be used to endorse or promote
+ products derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY LOS ALAMOS NATIONAL SECURITY, LLC AND CONTRIBUTORS "AS IS" AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL LOS
+ALAMOS NATIONAL SECURITY, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+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/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2010-05-11 05:15:27 UTC (rev 25)
+++ pkg/NAMESPACE 2010-05-12 04:44:34 UTC (rev 26)
@@ -8,6 +8,7 @@
export(possibility.envelope)
S3method(summary,dream)
S3method(coef,dream)
+S3method(fitted,dream)
S3method(plot,dream)
S3method(print,dream)
S3method(predict,dream)
Modified: pkg/R/CompDensity.R
===================================================================
--- pkg/R/CompDensity.R 2010-05-11 05:15:27 UTC (rev 25)
+++ pkg/R/CompDensity.R 2010-05-12 04:44:34 UTC (rev 26)
@@ -112,11 +112,13 @@
p <- sapply(temp,function(x) x[1])
logp <- sapply(temp,function(x) x[2])
- if(class(p)!="numeric") {
+ if (!is.numeric(p)) {
print(p)
stop(sprintf("Expected class numeric, got class %s. Error with multicore? Set control$parallel='none' to not use parallelisation",class(p)))
}
- stopifnot(!any(is.na(p)))
+ if (any(is.na(p))) {
+ stop("likelihood function produced invalid probabilities (NA/NaN)")
+ }
##stopifnot(!any(is.na(logp))) ##Not used anyway
return(list(p=p,logp=logp))
} ##CompDensity
Modified: pkg/R/dream.R
===================================================================
--- pkg/R/dream.R 2010-05-11 05:15:27 UTC (rev 25)
+++ pkg/R/dream.R 2010-05-12 04:44:34 UTC (rev 26)
@@ -1,5 +1,48 @@
+## Copyright (c) 2008, Los Alamos National Security, LLC
+## All rights reserved.
+##
+## Copyright 2008. Los Alamos National Security, LLC. This software was produced under U.S.
+## Government contract DE-AC52-06NA25396 for Los Alamos National Laboratory (LANL), which is
+## operated by Los Alamos National Security, LLC for the U.S. Department of Energy. The U.S.
+## Government has rights to use, reproduce, and distribute this software.
+##
+## NEITHER THE GOVERNMENT NOR LOS ALAMOS NATIONAL SECURITY, LLC MAKES A NY WARRANTY, EXPRESS OR
+## IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. If software is modified to
+## produce derivative works, such modified software should be clearly marked, so as not to
+## confuse it with the version available from LANL.
+##
+## Additionally, redistribution and use in source and binary forms, with or without
+## modification, are permitted provided that the following conditions are met:
+## * Redistributions of source code must retain the above copyright notice, this list of
+## conditions and the following disclaimer.
+## * Redistributions in binary form must reproduce the above copyright notice, this list of
+## conditions and the following disclaimer in the documentation and/or other materials
+## provided with the distribution.
+## * Neither the name of Los Alamos National Security, LLC, Los Alamos National Laboratory, LANL
+## the U.S. Government, nor the names of its contributors may be used to endorse or promote
+## products derived from this software without specific prior written permission.
+##
+## THIS SOFTWARE IS PROVIDED BY LOS ALAMOS NATIONAL SECURITY, LLC AND CONTRIBUTORS "AS IS" AND
+## ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+## OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL LOS
+## ALAMOS NATIONAL SECURITY, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+## SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+## SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+## 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
+################################################################################################
+## This R code has been converted and modified from the original MATLAB code by ##
+## Joseph Guillaume and Felix Andrews, 2010. ##
+################################################################################################
+
+
dreamDefaults <- function()
list(
nCR = 3, ## Crossover values used to generate proposals (geometric series)
@@ -270,17 +313,19 @@
test.pars <- FUN.pars
test.pars[[names(formals(FUN))[1]]] <- x[1,]
modpred <- do.call(FUN,test.pars)
- if (!inherits(modpred,"numeric")) stop(sprintf("Result of FUN should be of class numeric, not %s",class(modpred)))
+ if (!is.numeric(modpred))
+ stop(sprintf("Result of FUN should be of class numeric, not %s", class(modpred)))
## make each element of pars a list and extract lower / upper
lower <- sapply(pars, function(x) min(x[[1]]))
upper <- sapply(pars, function(x) max(x[[1]]))
##Step 2: Calculate posterior density associated with each value in x
- tmp<-do.call(CompDensity,list(pars=x,control=control,FUN=FUN,func.type=func.type,measurement=measurement,FUN.pars=FUN.pars))
+ tmp <- CompDensity(pars = x, control = control, FUN = FUN, func.type = func.type,
+ measurement = measurement, FUN.pars = FUN.pars)
##Save the initial population, density and log density in one list X
- X<-cbind(x=x,p=tmp$p,logp=tmp$logp)
+ X <- cbind(x = x, p = tmp$p, logp = tmp$logp)
colnames(X) <- c(names(pars), "p", "logp")
##Initialise the sequences
@@ -323,7 +368,8 @@
CR[,gen.number] <- tmp$CR
## Now compute the likelihood of the new points
- tmp <- CompDensity(pars=x.new,control=control,FUN=FUN,func.type=func.type,measurement=measurement,FUN.pars=FUN.pars)
+ tmp <- CompDensity(pars = x.new, control = control, FUN = FUN, func.type = func.type,
+ measurement = measurement, FUN.pars = FUN.pars)
p.new <- tmp$p
logp.new <- tmp$logp
@@ -480,12 +526,10 @@
obj$Sequences <- as.mcmc.list(lapply(1:NSEQ,function(i) as.mcmc(Sequences[,1:NDIM,i])))
if (!is.na(control$thin.t)){
Reduced.Seq <- Reduced.Seq[1:counter.redseq,,]
- obj$Reduced.Seq <- as.mcmc.list(lapply(1:NSEQ,function(i) mcmc(
- Reduced.Seq[,1:NDIM,i],
- start=1,
- end=counter-1,
- thin=control$thin.t)
- ))
+ obj$Reduced.Seq <- as.mcmc.list(lapply(1:NSEQ, function(i) {
+ mcmc(Reduced.Seq[,1:NDIM,i], start = 1,
+ end = counter-1, thin = control$thin.t)
+ }))
}
## TODO: make these 'ts' objects and sync with Reduced.Seq by thinning
Added: pkg/R/fitted.dream.R
===================================================================
--- pkg/R/fitted.dream.R (rev 0)
+++ pkg/R/fitted.dream.R 2010-05-12 04:44:34 UTC (rev 26)
@@ -0,0 +1,10 @@
+
+
+fitted.dream <-
+ function(object,
+ start = 1+(end(object$Sequences)-1)*(1-fraction),
+ fraction = 0.5,
+ ...)
+{
+ window(object$Sequences, start = start, ...)
+}
Modified: pkg/R/plot.dream.R
===================================================================
--- pkg/R/plot.dream.R 2010-05-11 05:15:27 UTC (rev 25)
+++ pkg/R/plot.dream.R 2010-05-12 04:44:34 UTC (rev 26)
@@ -4,39 +4,31 @@
##' Uses second half of sequences
-plot.dream <- function(x,interactive=TRUE,...){
- devAskNewPage(interactive)
+plot.dream <- function(x, interactive = TRUE, ...){
+ opar <- devAskNewPage(interactive)
+ on.exit(devAskNewPage(opar))
- ss <- window(x$Sequences, start = end(x$Sequences)/2 + 1)
+ ss <- fitted(x, ...)
- ## Trace and parameter density
+ ## Convergence (Gelman plot)
- plot(ss)
+ tmp <- try(gelman.plot(ss))
- print(xyplot(ss))
- print(densityplot(ss))
+ if (inherits(tmp, "try-error")) {
+ 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")
+ }
- ## Acceptance rate
- plot(table(x$AR[,2]),main="Distribution of % acceptance rate")
+ ## Trace and parameter density and auto-correlation
- ##Convergence
-
- try(gelman.plot(ss))
+ print(densityplot(ss))
- 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")
+ print(xyplot(ss))
- ## Multi-variate density for first chain
+ print(acfplot(ss))
+
+ ## Acceptance rate
+ print(barchart(table(x$AR[,2]), main="Distribution of % acceptance rate"))
- print(splom(as.data.frame(x$Sequences[[1]]),
- upper.panel = panel.smoothScatter, nrpoints = 0,
- lower.panel = function(x, y, ...) {
- panel.grid(-1, -1)
- panel.loess(x, y, span = 1/3, lwd = 1)
- panel.loess(x, y, span = 2/3, lwd = 2)
- grid::grid.text(paste("cor =", round(cor(x, y),2)),
- y = 0.1)
- }))
-
}##plot.dream
Modified: pkg/R/predict.dream.R
===================================================================
--- pkg/R/predict.dream.R 2010-05-11 05:15:27 UTC (rev 25)
+++ pkg/R/predict.dream.R 2010-05-12 04:44:34 UTC (rev 26)
@@ -37,7 +37,7 @@
##' @param use.thinned Whether to use existing thinned chains
##' @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,
+ method="uni.mode",level=0.99,
last.prop=0.5,use.thinned=TRUE,...
){
Modified: pkg/R/summary.dream.R
===================================================================
--- pkg/R/summary.dream.R 2010-05-11 05:15:27 UTC (rev 25)
+++ pkg/R/summary.dream.R 2010-05-12 04:44:34 UTC (rev 26)
@@ -1,7 +1,7 @@
-summary.dream <- function(object,...){
+summary.dream <- function(object, fraction = 0.5, ...){
- coda.sum <- summary(window(object$Sequences, start = end(object$Sequences)/2 + 1))
+ coda.sum <- summary(fitted(object, fraction = fraction), ...)
cat(sprintf("
Exit message: %s
@@ -18,7 +18,7 @@
if (all(R.stat.last<0)) {
R.stat.last <- gelman.diag(object$Sequences)$psrf[,1]
}
- names(R.stat.last) <-colnames(object$R.stat)[-1]
+ names(R.stat.last) <- colnames(object$R.stat)[-1]
for (i in 1:length(R.stat.last)){
cat(sprintf("\t%s:\t%f\n",
Modified: pkg/man/dream.Rd
===================================================================
--- pkg/man/dream.Rd 2010-05-11 05:15:27 UTC (rev 25)
+++ pkg/man/dream.Rd 2010-05-12 04:44:34 UTC (rev 26)
@@ -179,7 +179,7 @@
\references{
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
+ Hyman, J. M., Higdon, D., 2009. Accelerating Markov chain Monte Carlo
simulation by differential evolution with self-adaptive randomized
subspace sampling. \emph{International Journal of Nonlinear Sciences
and Numerical Simulation} 10 (3), 273-290.
More information about the Dream-commits
mailing list