[Dream-commits] r6 - in pkg: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Feb 17 08:24:18 CET 2010
Author: josephguillaume
Date: 2010-02-17 08:24:18 +0100 (Wed, 17 Feb 2010)
New Revision: 6
Added:
pkg/NAMESPACE
Modified:
pkg/R/AdaptpCR.R
pkg/R/CalcCbWb.R
pkg/R/CalcDelta.R
pkg/R/CompDensity.R
pkg/R/CovInit.R
pkg/R/dream.R
pkg/R/metrop.R
pkg/R/offde.R
Log:
Added namespace. Fixed numerous runtime errors using vrugt's example1. Still non-functional due to all chains being rejected from first iteration.
Added: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE (rev 0)
+++ pkg/NAMESPACE 2010-02-17 07:24:18 UTC (rev 6)
@@ -0,0 +1,6 @@
+import(coda)
+
+export(dream)
+export(CovInit)
+export(LHSInit)
+
Property changes on: pkg/NAMESPACE
___________________________________________________________________
Name: svn:executable
+ *
Modified: pkg/R/AdaptpCR.R
===================================================================
--- pkg/R/AdaptpCR.R 2010-02-15 06:48:13 UTC (rev 5)
+++ pkg/R/AdaptpCR.R 2010-02-17 07:24:18 UTC (rev 6)
@@ -10,6 +10,8 @@
##' lCR vector of length nCR
AdaptpCR <- function(CR,delta.tot,lCR.old,control){
+ stopifnot(sum(delta.tot)>0)
+
## dimensions:
## CR vector of length nseq*steps
## zz iter. 1:nCR
Modified: pkg/R/CalcCbWb.R
===================================================================
--- pkg/R/CalcCbWb.R 2010-02-15 06:48:13 UTC (rev 5)
+++ pkg/R/CalcCbWb.R 2010-02-17 07:24:18 UTC (rev 6)
@@ -11,6 +11,6 @@
## And use these to derive Cb and Wb
Cb <- (A1/A2)^(1/(1+beta))
Wb <- sqrt(A1)/((1+beta)*(A2^1.5))
- return(c(Cb=Cb,Wb=Wb))
+ return(list(Cb=Cb,Wb=Wb))
}
Modified: pkg/R/CalcDelta.R
===================================================================
--- pkg/R/CalcDelta.R 2010-02-15 06:48:13 UTC (rev 5)
+++ pkg/R/CalcDelta.R 2010-02-17 07:24:18 UTC (rev 6)
@@ -8,6 +8,8 @@
##' @return delta.tot vector of length nCR
CalcDelta <- function(nCR,delta.tot,delta.normX,CR){
+ stopifnot(sum(delta.tot)>0 || sum(delta.normX)>0)
+
## Dimensions:
## zz. iter 1:nCR
## idx. vector. length [0,nseq]. range [1,nseq]
@@ -22,5 +24,7 @@
delta.tot[zz] <- delta.tot[zz]+sum(delta.normX[idx])
} ## for CRs
+
+ stopifnot(!any(is.na(delta.tot)) && sum(delta.tot)>0)
return(delta.tot)
} ##CalcDelta
Modified: pkg/R/CompDensity.R
===================================================================
--- pkg/R/CompDensity.R 2010-02-15 06:48:13 UTC (rev 5)
+++ pkg/R/CompDensity.R 2010-02-17 07:24:18 UTC (rev 6)
@@ -5,7 +5,7 @@
##' @param FUN - the model to run
##' R function with first argument a vector of length ndim.
##' returns a scalar or vector corresponding to one of the options below.
-##' @param option Type of function output. One of:
+##' @param func.type Type of function output. One of:
##' posterior.density, logposterior.density,
##' calc.loglik. requires optional parameter measurement with elements data & sigma
##' calc.rmse, calc.weighted.rmse. requires measurement$data
@@ -24,7 +24,8 @@
measurement=NULL,...){
stopifnot(!is.null(measurement) || func.type%in% c("posterior.density","logposterior.density"))
-
+ stopifnot(!any(is.na(x)))
+
## dimensions:
## i. iter 1:nseq
## modpred. scalar or vector commensurate to measurement$data
@@ -86,5 +87,8 @@
logp[ii] <- -0.5*SSR
}) ##switch
} ## for rows
+
+ stopifnot(!any(is.na(p)))
+ stopifnot(!any(is.na(logp)))
return(list(p=p,logp=logp))
} ##CompDensity
Modified: pkg/R/CovInit.R
===================================================================
--- pkg/R/CovInit.R 2010-02-15 06:48:13 UTC (rev 5)
+++ pkg/R/CovInit.R 2010-02-17 07:24:18 UTC (rev 6)
@@ -5,8 +5,8 @@
##' @param muX vector of length ndim
##' @param qcov
##' @param bound.handling. character one of: none,reflect,bound,fold,rand
-##' @return ...
-##' x. dim nseq x ndim. parameter wise range [xmin,xmax] assured by handleBounds
+##' @return x
+##' matrix dim nseq x ndim. parameter wise range [xmin,xmax] assured by handleBounds
##
## depends on:
## handleBounds.R
Modified: pkg/R/dream.R
===================================================================
--- pkg/R/dream.R 2010-02-15 06:48:13 UTC (rev 5)
+++ pkg/R/dream.R 2010-02-17 07:24:18 UTC (rev 6)
@@ -17,7 +17,7 @@
REPORT = 10, ## number of iterations between reports when trace >= 1
thin=FALSE, ## do reduced sample collection
thin.t=NA, ## parameter for reduced sample collection
- metrop.opt=NA ## ?? which option to use for metropolis acceptance/rejection. Some also need measurement$sigma
+ ndim=NA ## number of parameters (automatically set from length of pars)
)
library(coda)
@@ -25,10 +25,28 @@
## MATLAB function:
## function [Sequences,Reduced_Seq,X,output,hist_logp] = dream(MCMCPar,ParRange,Measurement,ModelName,Extra,option)
+##' @param FUN model function with first argument a vector of length ndim
+##' @param func.type. one of posterior.density, logposterior.density
+##' @param pars a list of variable ranges
+##' @param INIT f(pars,nseq,...) returns nseq x ndim matrix of initial parameter values
+##' @param control
+##' @param measurement list N,sigma,data. must be included unless func.type=posterior.density or logposterior.density is selected
+
+##' @return ...
+##' TODO
+##' Sequences array n.elem*1.125 x ndim+2 x nseq
+##' Reduced.Seq array n.elem*1.125 x ndim+2 x nseq
+##' AR matrix n.elem x 2
+##' outlier vector of variable length
+##' R.stat matrix n.elem/steps x 1+ndim
+##' CR matrix n.elem/steps x 1+length(pCR)
+
+
dream <- function(FUN, func.type,
pars = list(x = range(0, 1e6)),
- ...,
+ FUN.pars=list(),
INIT = LHSInit,
+ INIT.pars=list(),
control = list(),
measurement=NULL
)
@@ -49,11 +67,18 @@
if (is.character(FUN))
FUN <- get(FUN, mode = "function")
stopifnot(is.function(FUN))
- stopifnot(is.list(par))
- stopifnot(length(par) > 0)
+ stopifnot(is.list(pars))
+ stopifnot(length(pars) > 0)
stopifnot(!is.null(measurement) || func.type %in% c("posterior.density","logposterior.density"))
+ stopifnot(control$boundHandling %in% c("reflect", "bound", "fold", "none"))
+ 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 %in% c("x",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))
############################
@@ -68,12 +93,13 @@
toString(names(control)[!isValid]))
## determine number of variables to be optimized
- control$ndim<-length(par)
+ control$ndim<-length(pars)
if (is.na(control$nseq)) control$nseq <- control$ndim
NDIM <- control$ndim
- NCR <- control$ncr
+ NCR <- control$nCR
NSEQ <- control$nseq
+
## for each iteration...
Iter <- NSEQ #? 1
@@ -112,14 +138,6 @@
############################
## Initialise output object
- ## dimensions:
- ## AR matrix n.elem x 2
- ## outlier vector of variable length
- ## R.stat matrix n.elem/steps x 1+ndim
- ## CR matrix n.elem/steps x 1+length(pCR)
- ## Sequences array n.elem*1.125 x ndim+2 x nseq
- ## Reduced.Seq array n.elem*1.125 x ndim+2 x nseq
-
obj <- list()
class(obj) <- c("dream", class(obj))
obj$call <- match.call()
@@ -149,6 +167,7 @@
## Check whether will save a reduced sample
if (control$thin){
+ iloc.2 <- 0
Reduced.Seq <- array(NA,c(floor(n.elem/control$thin.t),NDIM+2,NSEQ))
} else Reduced.Seq <- NULL
@@ -164,14 +183,16 @@
################################
## Step 1: Sample s points in the parameter space
- x <- INIT(pars, NSEQ,...)
+ x <- do.call(INIT,modifyList(INIT.pars,list(pars=pars,nseq=NSEQ)))
+
## 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<-CompDensity(x, control = control, FUN = FUN, func.type=func.type,measurement=measurement...)
+ tmp<-do.call(CompDensity,modifyList(FUN.pars,list(x=x,control=control,FUN=FUN,func.type=func.type,measurement=measurement)))
+
##Save the initial population, density and log density in one list X
X<-cbind(x=x,p=tmp$p,logp=tmp$logp)
if (!is.null(names(pars))) colnames(X) <- c(names(pars),"p","logp")
@@ -183,7 +204,7 @@
##Save N_CR in memory and initialize delta.tot
obj$CR[1,] <- c(Iter,pCR)
- delta.tot <- rep(NA,NCR)
+ delta.tot <- rep(0,NCR)
##Save history log density of individual chains
hist.logp[1,] <- X[,"logp"]
@@ -191,7 +212,8 @@
##Compute R-statistic. Using coda package
## TODO: more elegant way of using coda. And check for correctness
## TODO: alternatively, convert matlab implementation
- obj$R.stat[1,] <- c(Iter,gelman.diag(as.mcmc.list(lapply(1:NSEQ,function(i) as.mcmc(Sequences[1:iloc,1:NDIM,i])))))
+ ## n<10 matlab: -2 * ones(1,MCMCPar.n);
+ obj$R.stat[1,] <- c(Iter,rep(-2,NDIM))
################################
##Start iteration
@@ -199,36 +221,41 @@
for (gen.number in 1:control$steps) {
+ ## TODO: A logic error in CovInit, offde or CompDensity is causing everything to be rejected in metrop, even on the first iteration
+
## Initialize DR properties
new_teller <- new_teller + 1 ## counter for thinning
## Define the current locations and associated posterior densities
x.old <- X[,1:NDIM]
- p.old <- X[,1:(NDIM+1)]
- logp.old <- X[,1:(NDIM+2)]
+ p.old <- X[,NDIM+1]
+ logp.old <- X[,NDIM+2]
## Now generate candidate in each sequence using current point and members of X
+ ## Table.JumpRate appears to match matlab version
tmp <- offde(x.old, control = control,
CR=CR[,gen.number],
lower = lower, upper = upper,
Table.JumpRate=Table.JumpRate)
x.new <- tmp$x.new
+ stopifnot(!identical(x.new,x.old))
CR[,gen.number] <- tmp$CR
## Now compute the likelihood of the new points
- tmp <- CompDensity(x.new, control = control, FUN = FUN, ...)
+ tmp<-do.call(CompDensity,modifyList(FUN.pars,list(x=x.new,control=control,FUN=FUN,func.type=func.type,measurement=measurement)))
p.new <- tmp$p
logp.new <- tmp$logp
## Now apply the acceptance/rejectance rule for the chain itself
tmp <- metrop(x.new,p.new,logp.new,
x.old,p.old,logp.old,
- measurement,control
+ func.type,control,
+ measurement
)
newgen <- tmp$newgen
alpha12 <- tmp$alpha
accept <- tmp$accept
-
+ ## stopifnot(any(accept))
## NOTE: original MATLAB code had option for DR Delayed Rejection here)
## accept2,ItExtra not required
@@ -249,14 +276,20 @@
## And update X using current members of Sequences
X <- newgen; rm(newgen)
+
if (control$pCR.Update) {
## Calculate the standard deviation of each dimension of X
- ## TODO: matlab syntax is unclear - seems to be elementwise?
- r <- sd(c(X[,1:NDIM]))
+ ## TODO: matlab syntax is unclear - seems to be columnwise
+ ## element-wise: sd(c(X[,1:NDIM]))
+ r <- apply(X[,1:NDIM],2,sd)
## Compute the Euclidean distance between new X and old X
delta.normX <- rowSums(((x.old-X[,1:NDIM])/r)^2)
## Use this information to update sum_p2 to update N_CR
delta.tot <- CalcDelta(NCR,delta.tot,delta.normX,CR[,gen.number])
+
+ ##0s in delta.tot, delta.normX -> pCR has NaN in pCR.Update
+ stopifnot(any(delta.tot!=0) | any(delta.normX!=0))
+
}
## Update hist.logp
@@ -273,10 +306,10 @@
## ---------------------------------------------------------------------
## Store Important Diagnostic information -- Probability of individual crossover values
- obj$CR[teller, ] <- pCR
+ obj$CR[teller, ] <- c(Iter,pCR)
## Do this to get rounded iteration numbers
- if (teller == 2) steps <- steps + 1
+ if (teller == 2) control$steps <- control$steps + 1
## Check whether to update individual pCR values
if (Iter <= 0.1 * control$ndraw) {
@@ -313,7 +346,7 @@
## Calculate Gelman and Rubin convergence diagnostic
## Compute the R-statistic using 50% burn-in from Sequences
- obj$R.stat <- gelman.diag(as.mcmc.list(lapply(1:NSEQ,function(i) as.mcmc(Sequences[,1:NDIM,i]))),autoburnin=TRUE)
+ obj$R.stat <- gelman.diag(as.mcmc.list(lapply(1:NSEQ,function(i) as.mcmc(Sequences[i,1:NDIM,]))),autoburnin=TRUE)
## break if maximum time exceeded
toc <- as.numeric(Sys.time()) - tic
Modified: pkg/R/metrop.R
===================================================================
--- pkg/R/metrop.R 2010-02-15 06:48:13 UTC (rev 5)
+++ pkg/R/metrop.R 2010-02-17 07:24:18 UTC (rev 6)
@@ -6,8 +6,12 @@
##' @param x.old matrix nseq x ndim
##' @param p.old vector of length nseq
##' @param logp.old vector of length nseq
+##' @param func.type Type of function output. One of:
+##' posterior.density, logposterior.density,
+##' calc.loglik. requires optional parameter measurement with elements data & sigma
+##' calc.rmse, calc.weighted.rmse. requires measurement$data
+##' @param control list with elements:
##' @param measurement. needs N and sigma
-##' @param control list with elements:
##' gamma
##' metrop.opt range [1,5]
##' @return ... list with elements
@@ -17,9 +21,13 @@
## TODO: names for control$metrop.opt.
metrop<-function(x,p.x,logp.x,
x.old,p.old,logp.old,
- measurement,control
+ func.type,control,
+ measurement=NULL
){
+ stopifnot(!is.null(measurement) || func.type%in% c("posterior.density","logposterior.density"))
+ stopifnot(!any(is.na(p.x)))
+
## dimensions:
## nr.chains scalar. should be = control$nseq
## Z vector of length nseq range [0,1]
@@ -35,28 +43,33 @@
## And initialize accept with false
accept <- rep(FALSE,nr.chains)
- switch(control$metrop.opt,
- "1" = {
+ switch(func.type,
+ posterior.density = {
alpha <- min(p.x/p.old,1)
},
- "2" = { ## Lnp probability evaluation
+ calc.loglik = { ## Lnp probability evaluation
alpha <- min(exp(p.x-p.old),1)
},
- "3" = { ## SSE probability evaluation
+ calc.rmse = { ## SSE probability evaluation
alpha <- min((p.x/p.old)^(-measurement$N*(1+control$gamma)/2),1)
},
- "4" = { ## Lnp probability evaluation
+ logposterior.density = { ## Lnp probability evaluation
alpha <- min(exp(p.x-p.old),1)
},
- "5" = { ## Similar to 3 but now weighted with Measurement.Sigma
+ calc.weighted.rmse = { ## Similar to 3 but now weighted with Measurement.Sigma
## signs are different because we write -SSR
alpha <- min(exp(-0.5*(-p.x + p.old)/measurement$sigma^2),1);
- })
-
+ },
+ stop("Unrecognised value of control$metrop.opt")
+ )
+
## Generate random numbers
Z <- runif(nr.chains)
## Find which alpha's are greater than Z
idx <- which(Z<alpha)
+
+ ##stopifnot(length(idx)>0) ##Unlikely, but possible
+
## And update these chains
newgen[idx,] <- c(x[idx,],p.x[idx],logp.x[idx])
Modified: pkg/R/offde.R
===================================================================
--- pkg/R/offde.R 2010-02-15 06:48:13 UTC (rev 5)
+++ pkg/R/offde.R 2010-02-17 07:24:18 UTC (rev 6)
@@ -29,7 +29,7 @@
## delta.x. matrix nseq x ndim
## qq. iter. range [1,nseq]
## ii. vector. length nseq-1. range [1,nseq]
- ## rr. vector. length 2. range [1,nseq]
+ ## rr. vector. length [0,2*DEpairs]. range [1,nseq]
## i. vector. length [1,ndim]. range [1,ndim]
## JumpRate. scalar. range (0,~1.683]
## delta. vector. length ndim
@@ -54,7 +54,7 @@
noise.x<-control$eps * (2*rand(nseq,ndim)-1)
## Initialize the delta update
- delta.x<-matrix(NA,nseq,ndim)
+ delta.x<-matrix(0,nseq,ndim)
## Each chain evolves using information from other chains to create offspring
for (qq in 1:nseq){
@@ -63,10 +63,10 @@
ii <- (1:nseq)[-qq]
## randomly select two members of ii
- rr <- ii[tt[1:2*DEversion[qq],qq]]
+ rr <- ii[tt[1:(2*DEversion[qq]),qq]]
## --- WHICH DIMENSIONS TO UPDATE? DO SOMETHING WITH CROSSOVER ----
- i <- which(D[qq,]>(1-CR[qq,1]))
+ i <- which(D[qq,]>(1-CR[qq]))
## Update at least one dimension
if (length(i)==0) i <- sample(ndim,1)
@@ -81,7 +81,7 @@
JumpRate <- Table.JumpRate[NrDim,DEversion[qq]]
## Produce the difference of the pairs used for population evolution
- delta <- colSums(x.old[rr[1:DEversion[qq]],]-x.old[rr[DEversion[qq]+1:2*DEversion[qq]],])
+ delta <- colSums(x.old[rr[1:DEversion[qq]],,drop=FALSE]-x.old[rr[(DEversion[qq]+1):(2*DEversion[qq])],,drop=FALSE])
## Then fill update the dimension
delta.x[qq,i] <- (1+noise.x[qq,i])*JumpRate*delta[i]
@@ -112,11 +112,12 @@
## If delayed rejection step --> generate proposal with DR
## Loop over all chains -- all dimensions are updated
## Generate a new proposal distance using standard procedure
-
## Update x.old with delta_x and eps;
x.new <- x.old+delta.x+eps
x.new <- handleBounds(x.new,lower,upper,control$boundHandling)
+ stopifnot(!any(is.na(x.new)))
+ stopifnot(!any(is.na(CR)))
return(list(x.new=x.new,CR=CR))
}#function offde
More information about the Dream-commits
mailing list