[Dream-commits] r7 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Feb 19 06:39:19 CET 2010
Author: josephguillaume
Date: 2010-02-19 06:39:19 +0100 (Fri, 19 Feb 2010)
New Revision: 7
Modified:
pkg/R/AdaptpCR.R
pkg/R/CalcDelta.R
pkg/R/CovInit.R
pkg/R/dream.R
pkg/R/metrop.R
pkg/R/offde.R
Log:
Fixed numerous errors in several files. Now terminates successfully on vrugt's example 1 :-D
Modified: pkg/R/AdaptpCR.R
===================================================================
--- pkg/R/AdaptpCR.R 2010-02-17 07:24:18 UTC (rev 6)
+++ pkg/R/AdaptpCR.R 2010-02-19 05:39:19 UTC (rev 7)
@@ -10,7 +10,7 @@
##' lCR vector of length nCR
AdaptpCR <- function(CR,delta.tot,lCR.old,control){
- stopifnot(sum(delta.tot)>0)
+ ##stopifnot(sum(delta.tot)>0)
## dimensions:
## CR vector of length nseq*steps
Modified: pkg/R/CalcDelta.R
===================================================================
--- pkg/R/CalcDelta.R 2010-02-17 07:24:18 UTC (rev 6)
+++ pkg/R/CalcDelta.R 2010-02-19 05:39:19 UTC (rev 7)
@@ -8,7 +8,7 @@
##' @return delta.tot vector of length nCR
CalcDelta <- function(nCR,delta.tot,delta.normX,CR){
- stopifnot(sum(delta.tot)>0 || sum(delta.normX)>0)
+ ##stopifnot(sum(delta.tot)>0 || sum(delta.normX)>0)
## Dimensions:
## zz. iter 1:nCR
@@ -25,6 +25,6 @@
} ## for CRs
- stopifnot(!any(is.na(delta.tot)) && sum(delta.tot)>0)
+ ##stopifnot(!any(is.na(delta.tot)) && sum(delta.tot)>0)
return(delta.tot)
} ##CalcDelta
Modified: pkg/R/CovInit.R
===================================================================
--- pkg/R/CovInit.R 2010-02-17 07:24:18 UTC (rev 6)
+++ pkg/R/CovInit.R 2010-02-19 05:39:19 UTC (rev 7)
@@ -17,6 +17,11 @@
{
##[x] = repmat(Extra.muX,MCMCPar.seq,1) + randn(MCMCPar.seq,MCMCPar.n) * chol(Extra.qcov);
+
+ ## Components verified to match matlab
+ ## print(t(matrix(rep(muX,nseq),length(muX))))
+ ## print(chol(qcov))
+
x <- t(matrix(rep(muX,nseq),length(muX)))+randn(nseq,length(pars)) %*% chol(qcov)
lower <- sapply(pars, function(x) min(x[[1]]))
Modified: pkg/R/dream.R
===================================================================
--- pkg/R/dream.R 2010-02-17 07:24:18 UTC (rev 6)
+++ pkg/R/dream.R 2010-02-19 05:39:19 UTC (rev 7)
@@ -255,7 +255,7 @@
newgen <- tmp$newgen
alpha12 <- tmp$alpha
accept <- tmp$accept
- ## stopifnot(any(accept))
+ ## stopifnot(any(accept)) #Unlikely, but possible
## NOTE: original MATLAB code had option for DR Delayed Rejection here)
## accept2,ItExtra not required
@@ -288,7 +288,7 @@
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))
+ stopifnot(any(delta.tot!=0) | any(delta.normX!=0))
}
@@ -327,13 +327,13 @@
## Draw random other chain -- cannot be the same as current chain
r.idx <- which.max(tmp$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[,r.idx]
+ hist.logp[1:(counter-1),out.id] <- hist.logp[1:(counter-1),r.idx]
## Jump outlier chain to r_idx -- Sequences
Sequences[iloc,1:(NSEQ+2),out.id] <- X[r.idx,]
## Jump outlier chain to r_idx -- X
X[out.id,1:(NSEQ+2)] <- X[r.idx,]
## Add to chainoutlier
- outlier <- rbind(outlier,c(Iter,out.id))
+ obj$outlier <- rbind(obj$outlier,c(Iter,out.id))
} ##for remove outliers
} ##else
@@ -346,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[i,1:NDIM,]))),autoburnin=TRUE)
+ obj$R.stat[teller,] <- c(Iter,gelman.diag(as.mcmc.list(lapply(1:NSEQ,function(i) as.mcmc(Sequences[i,1:NDIM,]))),autoburnin=TRUE)$psrf[,1])
## break if maximum time exceeded
toc <- as.numeric(Sys.time()) - tic
@@ -360,26 +360,30 @@
teller = teller + 1
} ##while
+
+ obj$Sequences <- obj$Sequences
+ obj$Reduced.Seq <- obj$Reduced.Seq
+
## 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,,]
+ obj$Sequences <- obj$Sequences[1:i,,]
}
## Remove extra rows from Reduced.Seq
if (control$thin){
- i <- which(rowSums(Reduced.Seq)==0)
+ i <- which(rowSums(Reduced.Seq[,,1])==0)
if (length(i)>0) {
i <- i[1]-1
- obj$Reduced.Seq <- Reduced.Seq[1:i,,]
+ obj$Reduced.Seq <- obj$Reduced.Seq[1:i,,]
}
}
##Remove extra rows R.stat
- i <- which(rowSums(obj$R.stat[,,1])==0)
+ i <- which(is.na(rowSums(obj$R.stat)))
if (length(i)>0) {
i <- i[1]-1
- obj$R.stat <- obj$R.stat[1:i,,]
+ obj$R.stat <- obj$R.stat[1:i,]
}
##Remove extra rows AR
i <- which(rowSums(obj$AR)==0)
Modified: pkg/R/metrop.R
===================================================================
--- pkg/R/metrop.R 2010-02-17 07:24:18 UTC (rev 6)
+++ pkg/R/metrop.R 2010-02-19 05:39:19 UTC (rev 7)
@@ -18,7 +18,6 @@
##' newgen matrix nseq x ndim+2 (same as X in dream.R)
##' alpha scalar probability of acceptance. range [0,1]
##' accept vector indicating whether each sequences was accepted. length nseq
-## TODO: names for control$metrop.opt.
metrop<-function(x,p.x,logp.x,
x.old,p.old,logp.old,
func.type,control,
@@ -45,20 +44,20 @@
switch(func.type,
posterior.density = {
- alpha <- min(p.x/p.old,1)
+ alpha <- pmin(p.x/p.old,1)
},
calc.loglik = { ## Lnp probability evaluation
- alpha <- min(exp(p.x-p.old),1)
+ alpha <- pmin(exp(p.x-p.old),1)
},
calc.rmse = { ## SSE probability evaluation
- alpha <- min((p.x/p.old)^(-measurement$N*(1+control$gamma)/2),1)
+ alpha <- pmin((p.x/p.old)^(-measurement$N*(1+control$gamma)/2),1)
},
logposterior.density = { ## Lnp probability evaluation
- alpha <- min(exp(p.x-p.old),1)
+ alpha <- pmin(exp(p.x-p.old),1)
},
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);
+ alpha <- pmin(exp(-0.5*(-p.x + p.old)/measurement$sigma^2),1);
},
stop("Unrecognised value of control$metrop.opt")
)
@@ -66,7 +65,7 @@
## Generate random numbers
Z <- runif(nr.chains)
## Find which alpha's are greater than Z
- idx <- which(Z<alpha)
+ idx <- which(alpha>Z)
##stopifnot(length(idx)>0) ##Unlikely, but possible
Modified: pkg/R/offde.R
===================================================================
--- pkg/R/offde.R 2010-02-17 07:24:18 UTC (rev 6)
+++ pkg/R/offde.R 2010-02-19 05:39:19 UTC (rev 7)
@@ -6,7 +6,7 @@
##' @param lower see handleBounds
##' @param upper see handleBounds
##' @param Table.JumpRate. defined in dream. matrix ndim x DEpairs. range (0,~1.683]
-##' @return ...
+##' @return list with elments
## x.new, CR. same dim and range as input
##
## depends:
@@ -85,7 +85,7 @@
## Then fill update the dimension
delta.x[qq,i] <- (1+noise.x[qq,i])*JumpRate*delta[i]
-
+
} else {
## Set the JumpRate to 1 and overwrite CR and DEversion
JumpRate <- 1
@@ -96,16 +96,17 @@
## Now jumprate to facilitate jumping from one mode to the other in all dimensions
delta.x[qq,] <- JumpRate*delta
-
- ## Avoid that jump = 0 and xnew is similar to xold
- if (sum(delta.x[qq,]^2)==0){
- ## Compute the Cholesky Decomposition of x.old
- R <- (2.38/sqrt(nseq))*chol(cov(x.old)+1e-5*diag(ndim))
+ }##runif
+
+ ## Avoid that jump = 0 and xnew is similar to xold
+ if (sum(delta.x[qq,]^2)==0){
+ ## Compute the Cholesky Decomposition of x.old
+ R <- (2.38/sqrt(nseq))*chol(cov(x.old)+1e-5*diag(ndim))
- ## Generate jump using multinormal distribution
- delta.x[qq,] <- c( rnorm(ndim) %*% R )
- }
- }#runif
+ ## Generate jump using multinormal distribution
+ delta.x[qq,] <- c( rnorm(ndim) %*% R )
+ }
+
}#for qq
## TODO?:
More information about the Dream-commits
mailing list