[Dream-commits] r33 - in pkg: . R demo inst inst/extdata man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 16 06:53:11 CET 2010
Author: josephguillaume
Date: 2010-12-16 06:53:07 +0100 (Thu, 16 Dec 2010)
New Revision: 33
Added:
pkg/R/compareToMatlab.R
pkg/demo/example2.R
pkg/inst/extdata/
pkg/inst/extdata/example1a.mat
pkg/inst/extdata/example1b.mat
pkg/inst/extdata/example2a.mat
Removed:
pkg/demo/run_example1.R
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/RemOutlierChains.R
pkg/R/dream.R
pkg/demo/00Index
pkg/demo/example1.R
pkg/man/simulate.dream.Rd
Log:
- Allowed outlierTest='None' (and burnin.length=Inf) (don't run removeOutliers if burnin.length=0)
- Padding of numbers in default parameter names to guarantee alphabetic order
- Help file for simulate
- Added a example adapted from Vrugt's code
- Added undocumented functions to compare R output to the original matlab output
- Added mat6 files containing matlab output for example1 and example2
- Added ks.test and qq plot comparisons to matlab results for both example1 and example2. The R version appears to behave identically to the matlab version
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-06-11 01:42:19 UTC (rev 32)
+++ pkg/DESCRIPTION 2010-12-16 05:53:07 UTC (rev 33)
@@ -1,6 +1,6 @@
Package: dream
-Version: 0.3-2
-Date: 2010-06-11
+Version: 0.3-3
+Date: 2010-12-14
Title: DiffeRential Evolution Adaptive Metropolis
Author: Joseph Guillaume and Felix Andrews, based on MATLAB code by Jasper Vrugt
Maintainer: Joseph Guillaume <joseph.guillaume at anu.edu.au>
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2010-06-11 01:42:19 UTC (rev 32)
+++ pkg/NAMESPACE 2010-12-16 05:53:07 UTC (rev 33)
@@ -20,3 +20,10 @@
export(calc.weighted.rmse)
export(calc.loglik)
S3method(predict,dream_model)
+
+export(compareToMatlab,
+ plotMCMCQQ,
+ getMatlabSeq,
+ getMatlabControl,
+ writeMatlabDREAMSettings
+ )
Modified: pkg/R/RemOutlierChains.R
===================================================================
--- pkg/R/RemOutlierChains.R 2010-06-11 01:42:19 UTC (rev 32)
+++ pkg/R/RemOutlierChains.R 2010-12-16 05:53:07 UTC (rev 33)
@@ -79,6 +79,9 @@
chain.id <- idx
}
},
+ 'None'={
+ stop("Outlier detection reached when it should have been turned off")
+ },
stop("Unknown outlierTest specified")
) ##switch
Added: pkg/R/compareToMatlab.R
===================================================================
--- pkg/R/compareToMatlab.R (rev 0)
+++ pkg/R/compareToMatlab.R 2010-12-16 05:53:07 UTC (rev 33)
@@ -0,0 +1,161 @@
+compareToMatlab <- function(mat.file,
+ dream.obj){
+ library(R.matlab)
+ mat <- readMat(mat.file)
+ ## CHECK INPUTS
+ mat.control <- getMatlabControl(mat)
+ cat("\nDid matlab and R use same inputs?")
+ cat("\n names in Matlab not in R:")
+ cat(setdiff(names(mat.control),names(dream.obj$control)))
+ cat("\n names in R not in Matlab: ")
+ nn <- setdiff(names(dream.obj$control),names(mat.control))
+ nn <- nn[!nn %in% c("Rthres","parallel","burnin.length","maxtime","REPORT")]
+ cat(nn)
+ common.names <- intersect(names(dream.obj$control),names(mat.control))
+ mat.control <- mat.control[common.names]
+ r.control <- dream.obj$control[common.names]
+ cat("\n identical:",identical(mat.control,r.control))
+ cat("\n all.equal:",all.equal(mat.control,r.control))
+ eq <- sapply(1:length(r.control),
+ function(i) r.control[[i]]==mat.control[[i]])
+ names(eq) <- names(r.control)
+ cat("\n element-wise:\n")
+ print(eq)
+ cat("\n")
+ ## CHECK OUTPUTS
+ ## Compare sequences
+ mat.seq <- getMatlabSeq(mat)
+ R.seq <- dream.obj$Sequences
+ cat("\nDo outputs have same dimensions?\n")
+ print(sapply(list(matlab=mat.seq,R=R.seq),function(x) c(
+ variables=nvar(x),
+ chains=nchain(x),
+ iterations=niter(x),
+ thin=thin(x)
+ )))
+ cut.start=1 + (end(mat.seq) - 1) * (1 - 0.5)
+ mat.m <- as.matrix(window(mat.seq,start=cut.start,thin=100))
+ R.m <- as.matrix(window(R.seq,start=cut.start,thin=100))
+ cat("\n ks.test that samples are from same distribution for each variable\n")
+ pvals <- sapply(1:ncol(mat.m),function(i) ks.test(R.m[,i],mat.m[,i])$p.value)
+ names(pvals) <- colnames(R.m)
+ print(round(pvals,2))
+ ## Graphically
+ plotMCMCQQ(mat.m,R.m)
+}##compareToMatlab
+
+plotMCMCQQ <- function(mat.m,R.m){
+ library(reshape)
+ library(lattice)
+ colnames(mat.m) <- colnames(R.m)
+ mm <- rbind(
+ cbind(melt(mat.m),source="mat"),
+ cbind(melt(R.m),source="R")
+ )
+ mm2 <- cast(mm,X1+X2~source)
+ stopifnot(!any(is.na(mm2)))
+ print(xyplot(mat~R|X2,data=mm2,as.table=T,scales='free',
+ main="QQ plot of distribution of R and matlab code",
+ panel=function(x,y,...){
+ e <- sort(x) ##quantile(x,1:length(x)/length(x))
+ o <- sort(y) ##quantile(y,1:length(y)/length(y))
+ panel.points(e,o)
+ panel.lines(e,e)
+ }
+ ))
+}##plotMCMCQQ
+
+getMatlabSeq <-
+ function(mat) {
+ NSEQ <- dim(mat$Reduced.Seq)[3]
+ NDIM <- dim(mat$Reduced.Seq)[2]-2
+ as.mcmc.list(lapply(1:NSEQ,function(i)
+ mcmc(mat$Reduced.Seq[,1:NDIM,i],
+ thin=as.numeric(mat$Extra[dimnames(mat$Extra)[[1]]=="T"])
+ )))
+ }
+
+getMatlabControl <- function(mat){
+ MCMCPar <- lapply(mat$MCMCPar[,,1],as.vector)
+ Extra <- lapply(mat$Extra[,,1],function(x) ifelse(all(dim(x)==c(1,1)),as.vector(x),x))
+ list(
+ thin.t=Extra$T,
+ pCR.Update=Extra$pCR=="Update",
+ ndim=MCMCPar$n,
+ nseq=MCMCPar$seq,
+ ndraw=MCMCPar$ndraw,
+ nCR=MCMCPar$nCR,
+ gamma=MCMCPar$Gamma,
+ DEpairs=MCMCPar$DEpairs,
+ steps=MCMCPar$steps,
+ eps=MCMCPar$eps,
+ outlierTest=MCMCPar$outlierTest,
+ boundHandling=tolower(Extra$BoundHandling)
+ )
+}
+
+
+
+writeMatlabDREAMSettings <- function(dream.obj,ModelName,InitPopulation,
+ matlab.dream.dir,
+ in.mat.file="in.mat",
+ out.mat.file="out.mat",
+ run.m.file="run_once.m"
+ ){
+ library(R.matlab)
+ control <- dream.obj$control
+ Extra <- list(
+ pCR=as.matrix(ifelse(control$pCR.Update,"Update","Fixed")),
+ reduced_sample_collection="Yes",
+ T=control$thin.t,
+ InitPopulation=InitPopulation,
+ save_in_memory="No",
+ BoundHandling=paste(toupper(substring(control$boundHandling,1,1)),
+ substring(control$boundHandling,2),sep=""),
+ DR="No",
+ DRscale=10
+ )
+ MCMCPar <-
+ with(control,list(
+ n=as.numeric(ndim),
+ seq=nseq,
+ DEpairs=DEpairs,
+ Gamma=gamma,
+ nCR=nCR,
+ ndraw=ndraw,
+ steps=steps,
+ eps=eps,
+ outlierTest=outlierTest
+ ))
+ pars <- eval(dream.obj$call$pars)
+ ParRange <- list(
+ minn=matrix(sapply(pars, min),nrow=1),
+ maxn=matrix(sapply(pars, max),nrow=1)
+ )
+ Measurement <- list(
+ MeasData=NA,
+ Sigma=NA,
+ N=0
+ )
+ func.types <- c("posterior.density","calc.loglik","calc.rmse","logposterior.density","calc.weighted.rmse")
+ ## Write and run
+ oldwd <- setwd(matlab.dream.dir)
+ writeMat(in.mat.file,
+ Extra=Extra,
+ MCMCPar=MCMCPar,
+ ParRange=ParRange,
+ Measurement=Measurement,
+ ModelName=ModelName,
+ option=as.matrix(which(func.types==dream.obj$call$func.type))
+ )
+ cat(sprintf('
+load %s;
+[Sequences,Reduced_Seq,X,output,hist_logp] = dream(MCMCPar,ParRange,Measurement,ModelName,Extra,option);
+save -6 %s
+',in.mat.file,out.mat.file),file=run.m.file)
+ setwd(oldwd)
+ cat(sprintf("
+Please ensure that there is a file named '%s.m' in the directory '%s'.
+Run '%s' in Matlab and then run readMat('%s') in R to obtain Matlab's output.
+",ModelName,matlab.dream.dir,run.m.file,out.mat.file))
+} ##writeMatlabDREAMSettings
Property changes on: pkg/R/compareToMatlab.R
___________________________________________________________________
Added: svn:executable
+ *
Modified: pkg/R/dream.R
===================================================================
--- pkg/R/dream.R 2010-06-11 01:42:19 UTC (rev 32)
+++ pkg/R/dream.R 2010-12-16 05:53:07 UTC (rev 33)
@@ -135,8 +135,10 @@
stopifnot(is.list(pars))
stopifnot(length(pars) > 0)
pars <- lapply(pars, function(x) if (is.list(x)) x else list(x))
- if (is.null(names(pars)))
- names(pars) <- paste("p", 1:length(pars), sep = "")
+ if (is.null(names(pars))){
+ pad.length <- nchar(as.character(length(pars)))
+ names(pars) <- sprintf(paste("p%0",pad.length,"d",sep=""),1:length(pars))
+ }
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"))
@@ -180,6 +182,7 @@
control$REPORT <- (control$REPORT%/%control$nseq) * control$nseq
if (control$burnin.length<1) control$burnin.length <- control$burnin.length*control$ndraw
+ if (identical(tolower(control$outlierTest),'none')) control$burnin.length <- 0
##Choice of parallel backend
if (control$parallel!="none"){
@@ -252,6 +255,7 @@
end.burnin <- control$burnin.length
+
############################
## Initialise output object
@@ -430,7 +434,7 @@
## Do this to get rounded iteration numbers
if (counter.outloop == 2) control$steps <- control$steps + 1
- outliers <- RemOutlierChains(X,hist.logp[1:(counter-1),],control)
+ if (control$burnin.length!=0) outliers <- RemOutlierChains(X,hist.logp[1:(counter-1),],control)
if (counter.fun.evals <= end.burnin) {
## Check whether to update individual pCR values
Modified: pkg/demo/00Index
===================================================================
--- pkg/demo/00Index 2010-06-11 01:42:19 UTC (rev 32)
+++ pkg/demo/00Index 2010-12-16 05:53:07 UTC (rev 33)
@@ -1,4 +1,4 @@
FME.nonlinear.model Nonlinear model calibration as with FME vignette.
FME.nonlinear.model_parallelisation Example of parallelisation options
-run_example1 Function call for n-dimensional banana shaped gaussian distribution. From Vrugt's matlab code.
-example1 Setup for n-dimensional banana shaped gaussian distribution. From Vrugt's matlab code.
+example1 n-dimensional banana shaped gaussian distribution. From Vrugt's matlab code.
+example2 n-dimensional Gaussian distribution from Vrugt's matlab code
Modified: pkg/demo/example1.R
===================================================================
--- pkg/demo/example1.R 2010-06-11 01:42:19 UTC (rev 32)
+++ pkg/demo/example1.R 2010-12-16 05:53:07 UTC (rev 33)
@@ -1,4 +1,3 @@
-unloadNamespace("dream")
library(dream)
## n-dimensional banana shaped gaussian distribution
@@ -7,12 +6,14 @@
## @param bpar banana-ness. scalar.
## @param imat matrix ndim x ndim
## @return SSR. scalar
-## Cursorily verified to match matlab version
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,
@@ -25,22 +26,75 @@
outlierTest='IQR_test',
pCR.Update=TRUE,
thin.t=10,
- boundHandling='none'
+ boundHandling='none',
+ burnin.length=Inf, ##for compatibility with matlab code
+ REPORT=1e4 ##reduce frequency of progress reports
)
-pars <- lapply(1:control$ndim,function(x) c(-Inf,Inf))
-names(pars) <- paste("b",1:control$ndim,sep="")
+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
-FUN=Banshp
-func.type="logposterior.density"
-FUN.pars=list(
- imat=solve(cmat),
- bpar=0.1
- )
+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 <- 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))
+}
+
+## 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)
Added: pkg/demo/example2.R
===================================================================
--- pkg/demo/example2.R (rev 0)
+++ pkg/demo/example2.R 2010-12-16 05:53:07 UTC (rev 33)
@@ -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(system.file("extdata/example2a.mat",package="dream"),ddd)
+
Property changes on: pkg/demo/example2.R
___________________________________________________________________
Added: svn:executable
+ *
Deleted: pkg/demo/run_example1.R
===================================================================
--- pkg/demo/run_example1.R 2010-06-11 01:42:19 UTC (rev 32)
+++ pkg/demo/run_example1.R 2010-12-16 05:53:07 UTC (rev 33)
@@ -1,21 +0,0 @@
-
-source("example1.R")
-
-set.seed(11)
-
-dd <- dream(
- FUN=Banshp, func.type="logposterior.density",
- pars = pars,
- FUN.pars=list(
- imat=solve(cmat),
- bpar=0.1),
- INIT = CovInit,
- INIT.pars=list(
- muX=muX,
- qcov=qcov,
- bound.handling=control$boundHandling
- ),
- control = control
- )
-
-plot(dd)
Added: pkg/inst/extdata/example1a.mat
===================================================================
(Binary files differ)
Property changes on: pkg/inst/extdata/example1a.mat
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ application/octet-stream
Added: pkg/inst/extdata/example1b.mat
===================================================================
(Binary files differ)
Property changes on: pkg/inst/extdata/example1b.mat
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: pkg/inst/extdata/example2a.mat
===================================================================
(Binary files differ)
Property changes on: pkg/inst/extdata/example2a.mat
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Modified: pkg/man/simulate.dream.Rd
===================================================================
--- pkg/man/simulate.dream.Rd 2010-06-11 01:42:19 UTC (rev 32)
+++ pkg/man/simulate.dream.Rd 2010-12-16 05:53:07 UTC (rev 33)
@@ -2,3 +2,23 @@
\alias{simulate.dream}
\title{Simulate values from a distribution using converged MCMC chains
from a DREAM object}
+\usage{simulate(object,nsim=1000,seed=NULL,...)}
+\arguments{
+ \item{object}{\code{\link{dream}} object}
+ \item{nsim}{Approximate number of function evaluations (length of
+ sequence will be lower due to thinning)}
+ \item{seed}{passed to \code{\link{set.seed}} before continuing}
+ \item{...}{ignored, for compatibility with S3 generic}
+}
+\description{
+ Continues a converged MCMC chain to provide an additional
+ sample.
+
+ Outlier detection, reporting and convergence diagnostics are
+ turned off.
+
+ If no thinning was previously used, the sample is thinned with \code{t=10}. This value can be overridden by specifying a value for \code{thin.t} in the \code{control} parameter to \code{\link{dream}}.
+}
+\value{
+ a dream object with approximately the requested number of function evaluations
+}
More information about the Dream-commits
mailing list