[Dream-commits] r32 - in pkg: . R

[noreply at r-forge.r-project.org]

Author: josephguillaume
Date: 2010-06-11 03:42:19 +0200 (Fri, 11 Jun 2010)
New Revision: 32

Modified:
   pkg/DESCRIPTION
   pkg/R/CompDensity.R
   pkg/R/dream.R
Log:
- Changed returned hist.logp so that it now matches Sequences for each iteration - hist.logp did not match because of the way it was used for outlier detection
- Enforcing values for logposterior.density must be <=0

Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION	2010-05-17 02:38:57 UTC (rev 31)
+++ pkg/DESCRIPTION	2010-06-11 01:42:19 UTC (rev 32)
@@ -1,6 +1,6 @@
 Package: dream
-Version: 0.3-1
-Date: 2010-05-13
+Version: 0.3-2
+Date: 2010-06-11
 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/R/CompDensity.R
===================================================================
--- pkg/R/CompDensity.R	2010-05-17 02:38:57 UTC (rev 31)
+++ pkg/R/CompDensity.R	2010-06-11 01:42:19 UTC (rev 32)
@@ -45,7 +45,7 @@
            ## Model directly computes posterior density
            posterior.density={
              p <- modpred
-		 if (any(modpred<=0)) stop("Posterior density returned by FUN should be strictly positive. Otherwise use logposterior.density?")
+             if (any(modpred<0)) stop("Posterior density returned by FUN should be positive. Otherwise use logposterior.density?")
              logp <- log(modpred)
            },
            ## Model computes output simulation           
@@ -73,6 +73,7 @@
            logposterior.density={
              p <- modpred
              logp <- modpred
+             stopifnot(all(logp<=0))
            },
            ## Similar as 3, but now weights with the Measurement Sigma
            ## TODO: identical to rmse because difference is in metrop

Modified: pkg/R/dream.R
===================================================================
--- pkg/R/dream.R	2010-05-17 02:38:57 UTC (rev 31)
+++ pkg/R/dream.R	2010-06-11 01:42:19 UTC (rev 32)
@@ -233,6 +233,7 @@
   
   ## Initialize the array that contains the history of the log_density of each chain
   hist.logp <- matrix(NA_real_,max.counter,NSEQ)
+  real.hist.logp <- matrix(NA_real_,max.counter,NSEQ)
   
   if (control$pCR.Update){
     ## Calculate multinomial probabilities of each of the nCR CR values
@@ -338,6 +339,7 @@
   
   ##Save history log density of individual chains
   hist.logp[1,] <- X[,"logp"]
+  real.hist.logp[1,] <- X[,"logp"]
 
 ################################
   ##Start iteration
@@ -410,6 +412,7 @@
 
       ## Update hist.logp
       hist.logp[counter,] <- X[,NDIM+2]
+      real.hist.logp[counter,] <- X[,NDIM+2]
       
       ## Save Acceptance Rate
       obj$AR[counter,] <- c(counter.fun.evals,100 * sum(accept) / NSEQ)
@@ -447,6 +450,7 @@
         r.idx <- which.max(outliers$mean.hist.logp)
         ## Added -- update hist_logp -- chain will not be considered as an outlier chain then
         hist.logp[1:(counter-1),out.id] <- hist.logp[1:(counter-1),r.idx]
+        real.hist.logp[(counter-1),out.id] <- real.hist.logp[(counter-1),r.idx]
         ## Jump outlier chain to r_idx -- Sequences, X
         Sequences[(counter-1),1:(NDIM+2),out.id] <- X[r.idx,]
         X[out.id,1:(NDIM+2)] <- X[r.idx,]
@@ -536,7 +540,7 @@
   ##   See coef.dream for eg.
   obj$X <- X
   obj$R.stat <- obj$R.stat[1:counter.report,,drop=FALSE]
-  obj$hist.logp <- hist.logp[1:(counter-1),,drop=FALSE]
+  obj$hist.logp <- real.hist.logp[1:(counter-1),,drop=FALSE]
   obj$AR <- obj$AR[1:(counter-1),,drop=FALSE]
   obj$CR <- obj$CR[1:(counter.outloop-1),,drop=FALSE]