[Dream-commits] r3 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Feb 10 00:57:50 CET 2010
Author: josephguillaume
Date: 2010-02-10 00:57:49 +0100 (Wed, 10 Feb 2010)
New Revision: 3
Added:
pkg/R/ACR.R
pkg/R/AdaptpCR.R
pkg/R/CalcCbWb.R
pkg/R/CalcDelta.R
pkg/R/CompDensity.R
pkg/R/CovInit.R
pkg/R/DEStrategy.R
pkg/R/GenCR.R
pkg/R/LHSInit.R
pkg/R/RemOutlierChains.R
pkg/R/handleBounds.R
pkg/R/metrop.R
pkg/R/offde.R
Modified:
pkg/R/dream.R
Log:
Relatively complete first commit. May be broken due to dimension and data type errors.
Added: pkg/R/ACR.R
===================================================================
--- pkg/R/ACR.R (rev 0)
+++ pkg/R/ACR.R 2010-02-09 23:57:49 UTC (rev 3)
@@ -0,0 +1,75 @@
+## ACR Upper percentiles critical value for test of single multivariate normal outlier.
+## From the method given by Wilks (1963) and approaching to a F distribution
+## function by the Yang and Lee (1987) formulation, we provide an m-file to
+## get the critical value of the maximun squared Mahalanobis distance to detect
+## outliers from a normal multivariate sample.
+##
+## Syntax: function x = ACR(p,n,alpha)
+## $$ The function's name is giving as a gratefull to Dr. Alvin C. Rencher for his
+## unvaluable contribution to multivariate statistics with his text 'Methods of
+## Multivariate Analysis'.$$
+##
+## Inputs:
+## p - number of independent variables.
+## n - sample size.
+## alpha - significance level (default = 0.05).
+##
+## Output:
+## x - critical value of the maximun squared Mahalanobis distance.
+##
+## We can generate all the critical values of the maximun squared Mahalanobis
+## distance presented on the Table XXXII of by Barnett and Lewis (1978) and
+## Table A.6 of Rencher (2002). Also with any given significance level (alpha).
+##
+## Example: For p = 3; n = 25; alpha=0.01;
+##
+## Calling on Matlab the function:
+## ACR(p,n,alpha)
+##
+## Answer is:
+##
+## 13.1753
+##
+## Created by A. Trujillo-Ortiz, R. Hernandez-Walls, A. Castro-Perez and K. Barba-Rojo
+## Facultad de Ciencias Marinas
+## Universidad Autonoma de Baja California
+## Apdo. Postal 453
+## Ensenada, Baja California
+## Mexico.
+## atrujo at uabc.mx
+##
+## Copyright. August 20, 2006.
+##
+## To cite this file, this would be an appropriate format:
+## Trujillo-Ortiz, A., R. Hernandez-Walls, A. Castro-Perez and K. Barba-Rojo. (2006).
+## ACR:Upper percentiles critical value for test of single multivariate normal outlier.
+## A MATLAB file. [WWW document]. URL http://www.mathworks.com/matlabcentral/
+## fileexchange/loadFile.do?objectId=12161
+##
+## References:
+## Barnett, V. and Lewis, T. (1978), Outliers on Statistical Data.
+## New-York:John Wiley & Sons.
+## Rencher, A. C. (2002), Methods of Multivariate Analysis. 2nd. ed.
+## New-Jersey:John Wiley & Sons. Chapter 13 (pp. 408-450).
+## Wilks, S. S. (1963), Multivariate Statistical Outliers. Sankhya,
+## Series A, 25: 407-426.
+## Yang, S. S. and Lee, Y. (1987), Identification of a Multivariate
+## Outlier. Presented at the Annual Meeting of the American
+## Statistical Association, San Francisco, August 1987.
+##
+ACR <- function(p,n,alpha){
+
+if (missing(alpha)) alpha <- 0.05
+
+if (alpha<=0 | alpha>=1) stop("Significance level must be between 0 and 1")
+
+if (missing(p)|missing(n)) stop("Requires args p and n")
+
+a <- alpha
+## F distribution critical value with p and n-p-1 degrees of freedom using the Bonferroni correction
+Fc <- qf(1-a/n,p,n-p-1)
+ACR <- (p*(n-1)^2*Fc)/(n*(n-p-1)+(n*p*Fc))
+
+return(ACR)
+
+}##ACR
Property changes on: pkg/R/ACR.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/AdaptpCR.R
===================================================================
--- pkg/R/AdaptpCR.R (rev 0)
+++ pkg/R/AdaptpCR.R 2010-02-09 23:57:49 UTC (rev 3)
@@ -0,0 +1,27 @@
+## Updates the probabilities of the various crossover values
+## in:
+## lCR.old vector of length nCR
+AdaptpCR <- function(CR,delta.tot,lCR.old,control){
+ ## Make CR to be a single vector
+ CR <- c(CR)
+
+ ## Determine lCR
+ lCR <- rep(NA,control$nCR)
+ for (zz in 1:control$nCR){
+
+ ## Determine how many times a particular CR value is used
+ ## TODO: shouldn't this be which(CR==CR[z]])?
+ idx <- which(CR==zz/control$nCR)
+
+ ## This is used to weight delta.tot
+ lCR[zz] <- lCR.old[zz]+length(idx)
+ } #for CRs
+
+ ## Adapt pCR using information from averaged normalized jumping distance
+ pCR <- control$nseq*delta.tot/lCR / sum(delta.tot)
+
+ ## Normalize pCR
+ pCR <- pCR/sum(pCR)
+
+ return(list(pCR=pCR,lCR=lCR))
+} ##AdaptpCR
Property changes on: pkg/R/AdaptpCR.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/CalcCbWb.R
===================================================================
--- pkg/R/CalcCbWb.R (rev 0)
+++ pkg/R/CalcCbWb.R 2010-02-09 23:57:49 UTC (rev 3)
@@ -0,0 +1,13 @@
+CalcCbWb <- function(beta){
+ ## This function calculates the parameters for the exponential power density
+ ## Equation [20] paper by Thiemann et al. WRR 2001, Vol 37, No 10, 2521-2535
+
+ ## First calculate some dummy variables
+ A1 <- gamma(3*(1+beta)/2)
+ A2 <- gamma((1+beta)/2)
+ ## 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))
+}
+
Property changes on: pkg/R/CalcCbWb.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/CalcDelta.R
===================================================================
--- pkg/R/CalcDelta.R (rev 0)
+++ pkg/R/CalcDelta.R 2010-02-09 23:57:49 UTC (rev 3)
@@ -0,0 +1,18 @@
+## Calculate total normalized Euclidean distance for each crossover value
+## in:
+## delta.tot vector of length nCR ? TODO
+## delta.normX vector of length ndim ? TODO
+## CR vector of length nseq ? TODO
+CalcDelta <- function(control,delta.tot,delta.normX,CR){
+
+ ## Derive sum_p2 for each different CR value
+ for (zz in 1:control$nCR){
+ ## Find which chains are updated with zz/MCMCPar.nCR
+ ## TODO: possible that floating point error prevents exact comparison?
+ idx <- which(CR==zz/control$NCR)
+
+ ## Add the normalized squared distance tot the current delta_tot;
+ delta.tot[zz] <- delta.tot[zz]+sum(delta.normX[idx])
+ } ## for CRs
+ return(delta.tot)
+} ##CalcDelta
Property changes on: pkg/R/CalcDelta.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/CompDensity.R
===================================================================
--- pkg/R/CompDensity.R (rev 0)
+++ pkg/R/CompDensity.R 2010-02-09 23:57:49 UTC (rev 3)
@@ -0,0 +1,61 @@
+CompDensity <- function(x,control,FUN,option,
+ measurement,...){
+ ## This function computes the density of each x value
+ ## option is unambiguous abbrev of:
+ ## posterior.density, calc.loglik, calc.rmse, logposterior.density,calc.weighted.rmse
+
+ ## Sequential evaluation
+ for (ii in 1:nrow(x)){
+ ## Call model to generate simulated data
+ ## TODO: correct use of optional pars?
+ modpred <- FUN(...)
+
+ switch(option,
+ ## Model directly computes posterior density
+ posterior.density={
+ p[ii,1:2] <- c(ModPred,ii)
+ logp <- log(p[ii,1])
+ },
+ ## Model computes output simulation
+ calc.loglik={
+ err <- measurement$data-modpred
+
+ ## Compute the number of measurement data
+ N <- nrow(measurement)
+
+ ## Derive the log likelihood
+ logp[ii,1] <- N*log(control$Wb/measurement$sigma)-
+ control$CB*(sum((abs(Err/measurement$sigma))^(2/(1+control$gamma))))
+ ## And retain in memory
+ p[ii,1:2] <- c(logp[ii,1],ii)
+ },
+ ## Model computes output simulation
+ calc.rmse={
+ err <- measurement$data-modpred
+ ## Derive the sum of squared error
+ SSR <- sum(abs(err)^(2/(1+control$gamma)))
+ ## And retain in memory
+ p[ii,1:2] <- c(-SSR,ii)
+ logp[ii,1] <- -0.5*SSR
+
+ },
+ ## Model directly computes log posterior density
+ logposterior.density={
+ p[ii,1:2] <- [modpred ii]
+ logp[ii,1] <- p[ii,1]
+ },
+ ## Similar as 3, but now weights with the Measurement Sigma
+ ## TODO: appears to be no difference to calc.rmse
+ calc.weighted.rmse={
+ ## Define the error
+ err <- measurement$data-modpred
+ ## Derive the sum of squared error
+ SSR <- sum(abs(err)^(2/(1+control$gamma)))
+ ## And retain in memory
+ p[ii,2] <- c(-SSR,ii)
+ logp[ii,1] <- -0.5*SSR
+ }
+ )
+ }## for rows
+return(list(p=p,logp=logp))
+}##CompDensity
Property changes on: pkg/R/CompDensity.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/CovInit.R
===================================================================
--- pkg/R/CovInit.R (rev 0)
+++ pkg/R/CovInit.R 2010-02-09 23:57:49 UTC (rev 3)
@@ -0,0 +1,16 @@
+## in:
+## muX. vector of length ndim
+## out:
+## x. dim nseq x ndim
+CovInit <- function(pars, nseq,muX,qcov,bound.handling)
+{
+
+ ##[x] = repmat(Extra.muX,MCMCPar.seq,1) + randn(MCMCPar.seq,MCMCPar.n) * chol(Extra.qcov);
+ x <- t(matrix(rep(a,nseq),length(a)))+randn(control$nseq,length(pars)) %*% chol(qcov)
+
+ lower <- sapply(pars, function(x) min(x[[1]]))
+ upper <- sapply(pars, function(x) max(x[[1]]))
+
+ x <- handleBounds(x,lower,upper,bound.handling)
+ return(x)
+}
\ No newline at end of file
Property changes on: pkg/R/CovInit.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/DEStrategy.R
===================================================================
--- pkg/R/DEStrategy.R (rev 0)
+++ pkg/R/DEStrategy.R 2010-02-09 23:57:49 UTC (rev 3)
@@ -0,0 +1,20 @@
+DEStrategy<-function(control){
+ ## Determine which sequences to evolve with what DE strategy
+ ## in: control$ DEpairs, nseq
+ ## out: DEversion vector of length control$nseq
+
+ ## Determine probability of selecting a given number of pairs
+ p.pair <- (1/control$DEpairs)*rep(1,control$DEpairs)
+ p.pair <- c(0,cumsum(p.pair))
+
+ ## Generate a random number between 0 and 1
+ Z <- runif(control$nseq)
+
+ ## Select number of pairs
+ DEversion<-rep(0,control$nseq)
+ for (qq in 1:control$nseq){
+ z <- which(Z[qq]>p.pair)
+ DEversion[qq]<-tail(z,1)
+ }
+ return(DEversion)
+}#DEStrategy
Property changes on: pkg/R/DEStrategy.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/GenCR.R
===================================================================
--- pkg/R/GenCR.R (rev 0)
+++ pkg/R/GenCR.R 2010-02-09 23:57:49 UTC (rev 3)
@@ -0,0 +1,28 @@
+GenCR <- function(pCR,control){
+ ## Generates CR values based on current probabilities
+
+ ## How many candidate points for each crossover value?
+ # TODO: rmultinom may not be equivalent to multrnd
+ L <- rmultinom(1,size=control$nseq*control$steps,p=pCR)
+ L2 <- c(0,cumsum(L))
+
+ ## Then select which candidate points are selected with what CR
+ r <- sample(control$nseq*control$steps)
+
+ ## Then generate CR values for each chain
+ for (zz in 1:control$nCR){
+ ## Define start and end
+ i.start <- L2[zz]+1
+ i.end <- L2[zz+1]
+
+ ## Take the appropriate elements of r
+ idx <- r[i.start:i.end]
+
+ ## Assign these indices MCMCPar.CR(zz)
+ CR[idx] <- zz/control$nCR
+
+ ## Now reshape CR
+ CR <- array(CR,control$nseq,control$steps)\
+ } ## for nCR
+return(CR)
+} ## GenCR
Property changes on: pkg/R/GenCR.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/LHSInit.R
===================================================================
--- pkg/R/LHSInit.R (rev 0)
+++ pkg/R/LHSInit.R 2010-02-09 23:57:49 UTC (rev 3)
@@ -0,0 +1,28 @@
+## Latin Hypercube sampling
+LHSInit <- function(pars, nseq){
+ ##TODO: R method?
+ ## lapply(pars, function(r)
+ ## sample(seq(min(r), max(r), length = nseq))
+ ## )
+
+ xmin <- sapply(pars, function(x) min(x[[1]]))
+ xmax <- sapply(pars, function(x) max(x[[1]]))
+
+ ## Define the size of xmin
+ nvar <- length(xmin)
+ ## Initialize array ran with random numbers
+ ran <- rand(nseq,nvar)
+
+ ## Initialize array s with zeros
+ s <- matrix(0,nseq,nvar)
+
+ ## Now fill s
+ for (j in 1:nvar){
+ ## Random permutation
+ idx <- sample(nseq)
+ P <- idx-ran[,j]/nseq
+ s[,j] <- xmin[j]+P*(xmax[j]-xmin[j])
+ } ##for pars
+
+ return(x)
+} ## LHSInit
\ No newline at end of file
Property changes on: pkg/R/LHSInit.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/RemOutlierChains.R
===================================================================
--- pkg/R/RemOutlierChains.R (rev 0)
+++ pkg/R/RemOutlierChains.R 2010-02-09 23:57:49 UTC (rev 3)
@@ -0,0 +1,86 @@
+## Finds outlier chains and removes them when needed
+## in:
+## hist.logp matrix dim ndraw/nseq x nseq = length ndraw
+RemOutlierChains <- function(X,Sequences,hist.logp,Iter,outlier,control){
+ ## Determine the number of elements of L_density
+ idx.end <- nrow(hist.logp)
+ idx.start <- floor(0.5*idx.end)
+
+ ## Then determine the mean log density of the active chains
+ ## vector length nseq
+ mean.hist.logp <- colMeans(hist.logp[idx.start:idx.end,])
+
+ ## Initialize chain_id and Nid
+ chain.id <- NULL
+ Nid <- 0
+
+ ## Check whether any of these active chains are outlier chains
+ switch(control$outlierTest,
+ 'IQR_test'={
+ ## TODO: shouldn't upper.range be Q3+3*IQR?
+
+ ## Derive the upper and lower quantile of the data
+ q13<-quantile(mean.hist.logp,c(0.75,0.25))
+ ## Derive the Inter quartile range
+ iqr <- q13[1]-q12[2]
+ ## Compute the upper range -- to detect outliers
+ upper.range <- q13[2]-2*iqr
+ ## See whether there are any outlier chains
+ chain.id <- which(mean.hist.logp<upper.range)
+ Nid <- length(chain.id)
+ },
+ 'Grubbs_test'={
+ ## Test whether minimum log_density is outlier
+ G <- (mean(mean.hist.logp)-min(mean.hist.logp))/sd(mean.hist.logp)
+
+ ## Determine t-value of one-sided interval
+ t2 = tinv(1 - 0.01/control$nseq,control$nseq-2)^2; ## 95% interval
+
+ ## Determine the critical value
+ Gcrit <- ((control$nseq-1)/sqrt(control$nseq))*sqrt(t2/(control$nseq-2+t2))
+
+ ## Then check this
+ if (G > Gcrit) { ## Reject null-hypothesis
+ chain.id <- which.min(mean.hist.logp)
+ Nid <- 1
+ }
+ },
+ 'Mahal_test'={
+ ## Use the Mahalanobis distance to find outlier chains
+ alpha <- 0.01
+ upper.range <- ACR(control$ndim,control$nseq-1,alpha)
+ ## Find which chain has minimum log_density
+ idx <- which.min(mean.hist.logp)
+ ## Then check the Mahalanobis distance
+ d1 <- mahalanobis(X[idx,control$ndim],X[-idx,control$ndim])
+ ## Then see whether idx is an outlier in X
+ if (d1>upper.range) {
+ chain.id <- idx
+ Nid <- 1
+ }
+ }
+ )
+
+ if (Nid>0){
+ ## Loop over each outlier chain
+ for (qq in 1:Nid){
+ ## Draw random other chain -- cannot be the same as current chain
+ r.idx <- which.max(mean.hist.logp)
+ ## Added -- update hist_logp -- chain will not be considered as an outlier chain then
+ hist.logp[,chain.id[qq]] <- hist.logp[,r.idx]
+ ## Jump outlier chain to r_idx -- Sequences
+ Sequences[1,1:(control$nseq+2),chain.id[qq]] <- X[r.idx,1:(control$nseq+2)]
+ ## Jump outlier chain to r_idx -- X
+ X[chain.id[qq],1:(control$nseq+2)]] <- X[r.idx,1:(control$nseq+2)]
+ ## Add to chainoutlier
+ outlier <- rbind(outlier,c(Iter,chain.id[qq]))
+ }
+
+return(list(
+X=X,
+Sequences=Sequences,
+hist.logp=hist.logp,
+outlier=outlier
+))
+
+} ##RemOutlierChains
Property changes on: pkg/R/RemOutlierChains.R
___________________________________________________________________
Name: svn:executable
+ *
Modified: pkg/R/dream.R
===================================================================
--- pkg/R/dream.R 2010-02-08 23:41:34 UTC (rev 2)
+++ pkg/R/dream.R 2010-02-09 23:57:49 UTC (rev 3)
@@ -10,215 +10,320 @@
outlierTest = 'IQR_test', ## What kind of test to detect outlier chains?
pCR.Update = TRUE, ## Adaptive tuning of crossover values
boundHandling = 'reflect', ## Boundary handling: "reflect", "bound", "fold", "none"
- maxit = Inf, ## maximum number of iterations
+ ndraw = Inf, ## maximum number of iterations
maxeval = Inf, ## maximum number of function evaluations
maxtime = 60, ## maximum duration of optimization in seconds
trace = 0, ## level of user feedback
- REPORT = 10) ## number of iterations between reports when trace >= 1
+ REPORT = 10, ## number of iterations between reports when trace >= 1
+ thin=FALSE, ## do reduced sample collection
+ thin.t=NA ## parameter for reduced sample collection
+ )
-
-LHSInit <- function(pars, nseq)
-{
- lapply(pars, function(r)
- sample(seq(min(r), max(r), length = nseq))
- )
-}
-CovInit <- function(pars, nseq)
-{
-
-}
-
## MATLAB function:
# function [Sequences,Reduced_Seq,X,output,hist_logp] = dream(MCMCPar,ParRange,Measurement,ModelName,Extra,option)
dream <- function(FUN, pars = list(x = range(0, 1e6)),
...,
INIT = LHSInit,
- thin = 0, control = list()) #, do.SSE = FALSE)
+ control = list()) #, do.SSE = FALSE)
{
- if (is.character(FUN))
- FUN <- get(FUN, mode = "function")
- stopifnot(is.function(FUN))
- stopifnot(is.list(par))
- stopifnot(length(par) > 0)
- stopifnot(is.numeric(thin))
+ If (is.character(FUN))
+ FUN <- get(FUN, mode = "function")
+ stopifnot(is.function(FUN))
+ stopifnot(is.list(par))
+ stopifnot(length(par) > 0)
- ## determine number of variables to be optimized
- NDIM <- length(par)
+ ## update default options with supplied options
+ stopifnot(is.list(control))
+ control <- modifyList(dreamDefaults(), control)
+ isValid <- names(control) %in% names(dreamDefaults())
+ if (any(!isValid))
+ stop("unrecognised options: ",
+ toString(names(control)[!isValid]))
- ## update default options with supplied options
- stopifnot(is.list(control))
- control <- modifyList(dreamDefaults(), control)
- isValid <- names(control) %in% names(dreamDefaults())
- if (any(!isValid))
- stop("unrecognised options: ",
- toString(names(control)[!isValid]))
+ ## Initialize variables
+
+ ## determine number of variables to be optimized
+ control$ndim<-length(par)
- nseq <- control$nseq
- nCR <- control$nCR
+ ## for each iteration...
+ Iter <- nseq #? 1
+ counter <- 2
+ iloc <- 1
+ teller <- 2
+ new_teller <- 1
+
+ ## Calculate the parameters in the exponential power density function of Box and Tiao (1973)
+ cbwb <- CalcCbWb(control$gamma)
+ control$Cb<-cbwb$cb
+ control$Wb <- cbwb$wb
+
+ ## Initialize the array that contains the history of the log_density of each chain
+ hist.logp<-matrix(NA,floor(control$ndraw/control$nseq),control$nseq)
+
+ ## Derive the number of elements in the output file
+ n.elem<-floor(control$ndraw/control$seq)+1
- ## Sample s points in the parameter space
- x <- INIT(pars, nseq)
+## the output object
+ obj <- list()
+ class(obj) <- c("dream", class(obj))
+ obj$call <- match.call()
+ obj$control <- control
- ## make each element of pars a list and extract lower / upper
- pars <- lapply(pars, function(x) if (is.list(x)) x else list(x))
- lower <- sapply(pars, function(x) min(x[[1]]))
- upper <- sapply(pars, function(x) max(x[[1]]))
+ EXITFLAG <- NA
+ EXITMSG <- NULL
+
+ ## Initialize output information -- AR
+ obj$AR<-matrix(0,n.elem,2)
+ obj$AR[1,1:2]<-control$nseq-1
+
+ ## Initialize output information -- Outlier chains
+ obj$outlier<-c()
+
+ ## Initialize output information -- R statistic
+ obj$R.stat<-matrix(floor(n.elem/control$steps),control$ndim+1)
- x <- handleBounds(x, lower, upper, control$bound.handling)
+ if (control$pCR.Update){
+ ## Calculate multinomial probabilities of each of the nCR CR values
+ pCR <- (1/control$nCR)*rep(1,control$nCR)
+
+ ## Calculate the actual CR values based on p
+ CR <- GenCR(control,pCR)
+ lCR <- rep(0,control$nCR)
+ } else {
+ pCR <- 1/control$NCR
+ ## Define
+ CR <- pCR*matrix(1,control$nseq,control$steps)
+ lCR <- control$nseq*control$steps
+ } ##pCR.Update
+
+ ## Initialize output information -- N_CR
+ obj$CR <- matrix(0,floor(n.elem/control$steps),length(pCR)+1)
-
- ## initialize variables
- if (isTRUE(control$saveInMemory)) {
- ## support this case?
- } else {
- Sequences <- matrix(as.numeric(NA), nrow = nseq, ncol = NDIM+2)
- }
- if (!is.null(names(pars)))
- colnames(Sequences) <- names(pars)
+ if (isTRUE(control$saveInMemory)) {
+ ## TODO: support this case?
+ Sequences <- array(NA, c(floor(1.25*n.elem),control$ndim+1,control$nseq))
+ ## TODO: embellishments. accurate?
+ if (!is.null(names(pars))) colnames(Sequences) <- names(pars)
Sequences[1,] <- sapply(pars, mean)
+ } else {
+ Sequences<-NULL
+ }
- ## the output object
- obj <- list()
- class(obj) <- c("dream", class(obj))
- obj$call <- match.call()
- obj$control <- control
+ ## Generate the Table with JumpRates (dependent on number of dimensions and number of pairs)
+ Table.JumpRate<-matrix(NA,control$ndim,control$DEpairs)
+ for (zz in 1:control$DEpairs) Table.JumpRate[,zz] <- 2.38/sqrt(2*zz*1:control$ndim)
- EXITFLAG <- NA
- EXITMSG <- NULL
+ ## Check whether will save a reduced sample
+ # TODO
+
+ ## Change MCMCPar.steps to make sure to get nice iteration numbers in first loop
+ control$steps<-control$steps-1
+
+ ## initialize timer
+ tic <- as.numeric(Sys.time())
+ toc <- 0
- ## initialize timer
- tic <- as.numeric(Sys.time())
- toc <- 0
+
+ ## Step 1: Sample s points in the parameter space
+ x <- INIT(pars, nseq,...)
- ## for each iteration...
- Iter <- nseq #? 1
- counter <- 2
- iloc <- 1
- teller <- 2
- new_teller <- 1
+ ## make each element of pars a list and extract lower / upper
+ pars <- lapply(pars, function(x) if (is.list(x)) x else list(x))
+ 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
+ X<-CompDensity(x, control = control, FUN = FUN, ...)
+ #Save the initial population, density and log density in one list X
+ X<-cbind(x=x,p=X$p[,1],logp=X$logp)
- while (Iter < MAXIT) {
+ #Initialise the sequences
+ for (qq in 1:control$nseq){
+ Sequences[1,1:(control$ndim+2),qq] <- X[qq,]
+ }
- for (gen_number in 1:control$steps) {
- new_teller <- new_teller + 1 ## counter for thinning
+ #Save N_CR in memory and initialize delta_tot
+ obj$CR[1,1:(length(pCR)+1)] <- c(Iter,pCR)
+ delta.tot <- rep(NA,control$nCR)
+
+ #Save history log density of individual chains
+ hist.logp[1,1:control$nseq] <- c(Iter, X$logp)
+
+ #Compute R-statistic. Using coda package
+ # TODO: gelman.diag takes mcmc.list as input
+ obj$R.stat[1,1:(control$ndim+1)] <- c(Iter,gelman.diag(Sequences(1:iloc,1:control$ndim,1:control$nseq)))
+
+ ##Start iteration
+ while (Iter < control$ndraw) {
- ## Define the current locations and associated posterior densities
- xold <- GetLocation(X, control)
+ for (gen.number in 1:control$steps) {
- ## Now generate candidate in each sequence using current point and members of X
- xnew <- offde(xold, X, control = control, lower = lower, upper = upper)
+ ## Initialize DR properties
+ new_teller <- new_teller + 1 ## counter for thinning
- ## Now compute the likelihood of the new points
- xnew <- CompDensity(xnew, control = control, FUN = FUN, ...)
+ ## Define the current locations and associated posterior densities
+ x.old <- X[,1:control$ndim]
+ p.old <- X[,1:(control$ndim+1)]
+ logp.old <- X[,1:(control$ndim+2)]
- ## Now apply the acceptance/rejectance rule for the chain itself
- tmp <- metrop(xnew, xold, ETC)
+ ## Now generate candidate in each sequence using current point and members of X
+ tmp <- offde(x.old, control = control,
+ CR=CR,
+ lower = lower, upper = upper,bound.handling=bound.handling,
+ Table.JumpRate=Table.JumpRate)
+ x.new <- tmp$x.new
+ CR[,gen.number] <- tmp$CR
- ## (NOTE: original MATLAB code had option for DR Delayed Rejection here)
+ ## Now compute the likelihood of the new points
+ tmp <- CompDensity(xnew, control = control, FUN = FUN, ...)
+ p.new <- tmp$p
+ logp.new <- tmp$logp
- ## Now update the locations of the Sequences with the current locations
- if (isTRUE(control$saveInMemory))
- iloc <- iloc + 1
- #Sequences[iloc, 1:(NDIM+2), 1:nseq] <- matrix(newgen, NDIM+2, nseq)
- Sequences[1:(NDIM+2), 1:nseq] <- matrix(newgen, NDIM+2, nseq)
+ ## Now apply the acceptance/rejectance rule for the chain itself
+ tmp <- metrop(x.new,p.new,logp.new,
+ x.old,p.old,logp.old,
+ N,sigma,control,option
+ )
+ newgen <- tmp$newgen
+ alpha12 <- tmp$alpha
+ accept <- tmp$accept
+
- ## Check reduced sample collection
- if (thin > 0) {
- ## TODO
- }
+ ## NOTE: original MATLAB code had option for DR Delayed Rejection here)
+ ## accept2,ItExtra not required
- ## And update X using current members of Sequences
- X <- newgen; rm(newgen)
+ ## Update location in sequence and update the locations of the Sequences with the current locations
+ if (isTRUE(control$saveInMemory)) iloc <- iloc + 1
+ ##Sequences[iloc, 1:(NDIM+2), 1:nseq] <- matrix(newgen, NDIM+2, nseq)
+ Sequences[1:(NDIM+2), 1:nseq] <- matrix(newgen, NDIM+2, nseq)
- if (control$pCR.Update) {
- ## Calculate the standard deviation of each dimension of X
+ ## Check reduced sample collection
+ if (thin && new_teller == thin.t){
+ ## Update iloc_2 and new_teller
+ iloc.2 <- iloc.2+1
+ new_teller <- 0
+ ## Reduced sample collection
+ Reduced.Seq[iloc.2,1:(control$ndim+2),1:control$nseq] <-
+ array(newgen,c(1,control$ndim+2,control$nseq))
+ }
- ## Compute the Euclidean distance between new X and old X
+ ## And update X using current members of Sequences
+ X <- newgen; rm(newgen)
- ## Use this information to update sum_p2 to update N_CR
+ 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:control$ndim]))
+ ## Compute the Euclidean distance between new X and old X
+ delta.normX <- colSums(((x.old-X[,1:control$ndim])/r)^2,2)
+ ## Use this information to update sum_p2 to update N_CR
+ delta.tot <- CalcDelta(control,delta.tot,delta.normX,CR[,gen.number])
+ }
- }
+ ## Update hist.logp
+ hist.logp[counter,1:(control$nseq+1)] <- c(Iter,X[,control$ndim+2])
+
+ ## Save some important output -- Acceptance Rate
+ obj$AR[counter] <- 100 * sum(accept) / nseq
- ## Update hist_logp
+ ## Update Iteration and counter
+ Iter <- Iter + control$nseq
+ counter <- counter + 1
+ } ##for gen.number steps
- ## Save some important output -- Acceptance Rate
- obj$AR[counter] <- 100 * sum(accept) / nseq
+ ## ---------------------------------------------------------------------
- ## Update Iteration and counter
- Iter <- Iter + nseq
- counter <- counter + 1
- }
+ ## Store Important Diagnostic information -- Probability of individual crossover values
+ obj$CR[teller, ] <- pCR
- ## ---------------------------------------------------------------------
+ ## Do this to get rounded iteration numbers
+ if (teller == 2) steps <- steps + 1
- ## Store Important Diagnostic information -- Probability of individual crossover values
- obj$CR[teller, ] <- pCR
+ ## Check whether to update individual pCR values
+ if (Iter <= 0.1 * control$ndraw) {
+ if (control$pCR.Update) {
+ ## Update pCR values
+ tmp <- AdaptpCR(CR, delta_tot, lCR, control)
+ pCR <- tmp$pCR
+ lCR <- tmp$lCR
+ }
+ } else {
+ ## See whether there are any outlier chains, and remove them to current best value of X
+ tmp <- RemOutlierChains(X,
+ Sequences[iloc,,drop=FALSE],
+ hist.logp[1:(counter-1),2:(control$nseq+1)],Iter,outlier,control
+ )
+ Sequences[iloc,,] <- tmp$Sequences
+ X <- tmp$X
+ hist.logp[1:(counter-1),2:(control$nseq+1)] <- tmp$hist.logp
+ obj$outlier <- tmp$outlier
+ }
- ## Do this to get rounded iteration numbers
- if (teller == 2)
- steps <- steps + 1
+ if (control$pCR.Update) {
+ ## Generate CR values based on current pCR values
+ CR <- GenCR(pCR, control = control)
+ } else {
+ CR <- matrix(pCR, nrow = nseq, ncol = steps)
+ }
- ## Check whether to update individual pCR values
- if (Iter <= 0.1 * MAXIT) {
- if (control$pCR.Update) {
- ## Update pCR values
- tmp <- AdaptpCR(CR, delta_tot, lCR, control)
- pCR <- tmp$pCR
- lCR <- tmp$lCR
- }
- } else {
- ## See whether there are any outlier chains, and remove them to current best value of X
- outliers <- OutlierChains(X, Sequences[iloc,,drop=FALSE],
- ..., ETC)
- Sequences[iloc,,] <- Sequences[iloc,,] outliers
- ETC
- }
+ ## Calculate Gelman and Rubin convergence diagnostic
+ ## Compute the R-statistic using 50% burn-in from Sequences
+ ## TODO: NQR
+ obj$R.stat <- gelman.diag(Sequences[seq(iloc %/% 2, iloc), control$ndim, control$nseq))
- if (control$pCR.Update) {
- ## Generate CR values based on current pCR values
- CR <- GenCR(pCR, control = control)
- } else {
- CR <- matrix(pCR, nrow = nseq, ncol = steps)
- }
+ ## break if maximum time exceeded
+ toc <- as.numeric(Sys.time()) - tic
+ if (toc > control$maxtime) {
+ EXITMSG <- 'Exceeded maximum time.'
+ EXITFLAG <- 2
+ break
+ }
- ## Calculate Gelman and Rubin convergence diagnostic
- ## Compute the R-statistic using 50% burn-in from Sequences
- obj$Rstat <- Gelman(Sequences[seq(iloc %/% 2, iloc), NDIM, nseq, control))
+ ## Update the teller
+ teller = teller + 1
+ } ##while
- ## break if maximum time exceeded
- toc <- as.numeric(Sys.time()) - tic
- if (toc > control$maxtime) {
- EXITMSG <- 'Exceeded maximum time.'
- EXITFLAG <- 2
- break
- }
+ ## Postprocess output from DREAM before returning arguments
+ ## Remove extra rows from Sequences
+ i <- which(rowSums(Sequences[,,1])==0)
+ if (length(i)>0) {
+ i <- i[1]-1
+ obj$Sequences <- Sequences[1:i,,]
+ }
+ ## Remove extra rows from Reduced.Seq
+ i <- which(rowSums(Reduced.Seq)==0)
+ if (length(i)>0) {
+ i <- i[1]-1
+ obj$Reduced.Seq <- Reduced.Seq[1:i,,]
+ }
+ ##Remove extra rows R.stat
+ i <- which(rowSums(obj$R.stat[,,1])==0)
+ if (length(i)>0) {
+ i <- i[1]-1
+ obj$R.stat <- obj$R.stat[1:i,,]
+ }
+ ##Remove extra rows AR
+ i <- which(rowSums(obj$AR)==0)
+ if (length(i)>0) {
+ i <- i[1]-1
+ obj$AR <- obj$AR[1:i,]
+ }
+ ## Remove extra rows from CR
+ i <- which(rowSums(obj$CR)==0)
+ if (length(i)>0) {
+ i <- i[1]-1
+ obj$CR <- obj$CR[1:i,]
+ }
+
+ ## store number of function evaluations
+ obj$counts <- funevals
+ ## store number of iterations
+ obj$iterations <- i
+ ## store the amount of time taken
+ obj$time <- toc
- ## Update the teller
- teller = teller + 1
- }
+ obj
+}##dream
- ## store number of function evaluations
- obj$counts <- funevals
- ## store number of iterations
- obj$iterations <- i
- ## store the amount of time taken
- obj$time <- toc
- obj
-}
-
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/dream -r 3
More information about the Dream-commits
mailing list