[Dream-commits] r25 - in pkg: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue May 11 07:15:28 CEST 2010
Author: felix
Date: 2010-05-11 07:15:27 +0200 (Tue, 11 May 2010)
New Revision: 25
Modified:
pkg/DESCRIPTION
pkg/R/coef.dream.R
pkg/R/dream.R
pkg/man/coef.dream.Rd
pkg/man/dream.Rd
Log:
store hist.logp; coef method "sample.ml"
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-05-07 08:13:54 UTC (rev 24)
+++ pkg/DESCRIPTION 2010-05-11 05:15:27 UTC (rev 25)
@@ -1,6 +1,6 @@
Package: dream
-Version: 0.2-0
-Date: 2010-02-19
+Version: 0.2-1
+Date: 2010-05-11
Title: DiffeRential Evolution Adaptive Metropolis
Author: Jasper Vrugt, CJF ter Braak, et al. (R port by Joseph Guillaume)
Maintainer: Joseph Guillaume <joseph.guillaume at anu.edu.au>
Modified: pkg/R/coef.dream.R
===================================================================
--- pkg/R/coef.dream.R 2010-05-07 08:13:54 UTC (rev 24)
+++ pkg/R/coef.dream.R 2010-05-11 05:15:27 UTC (rev 25)
@@ -2,10 +2,10 @@
##' @param dream object
##' @param last.prop proportion of total sequence to use (0,1]
##' if 1, use whole sequence
-##' @param method. either a function or one of maxLik,mean,median
+##' @param method. either a function or one of uni.mode,mean,median,sample.ml
##' @return named vector of parameter values
coef.dream <- function(object,last.prop=.5,use.thinned=FALSE,
- method=c("maxLik","mean","median"),...){
+ method=c("uni.mode","mean","median","sample.ml"),...){
stopifnot(last.prop>0)
@@ -19,10 +19,19 @@
stopifnot(!is.null(sss))
- if (class(method)!="function") {
+ if (identical(method, "sample.ml")) {
+ ## TODO: make sure ppp corresponds to sss
+ ppp <- object$hist.logp
+ maxi <- which.max(ppp)
+ maxchain <- col(ppp)[maxi]
+ maxtime <- row(ppp)[maxi]
+ return(sss[[maxchain]][maxtime,])
+ }
+
+ if (!is.function(method)) {
method <- switch(
match.arg(method),
- "maxLik"=maxLikCoda,
+ "uni.mode"=maxLikCoda,
"mean"=function(sss) colMeans(as.matrix(sss)),
"median"=function(sss) apply(as.matrix(sss),2,median)
)
Modified: pkg/R/dream.R
===================================================================
--- pkg/R/dream.R 2010-05-07 08:13:54 UTC (rev 24)
+++ pkg/R/dream.R 2010-05-11 05:15:27 UTC (rev 25)
@@ -92,6 +92,8 @@
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 = "")
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"))
@@ -188,7 +190,7 @@
for (zz in 1:control$DEpairs) Table.JumpRate[,zz] <- 2.38/sqrt(2*zz*1:NDIM)
## Initialize the array that contains the history of the log_density of each chain
- hist.logp<-matrix(NA,max.counter,NSEQ)
+ hist.logp <- matrix(NA_real_,max.counter,NSEQ)
if (control$pCR.Update){
## Calculate multinomial probabilities of each of the nCR CR values
@@ -223,15 +225,14 @@
## counter.fun.evals + AR at each step
obj$AR<-matrix(NA,max.counter,2)
obj$AR[1,2]<-NSEQ-1 ##Number if only one rejected
- colnames(obj$AR) <- c("fun.evals","AR")
+ colnames(obj$AR) <- c("fun.evals", "AR")
##counter.fun.evals + R statistic for each variable at each step
## TODO: now using counter.report
obj$R.stat<-matrix(NA,max.counter.outloop,NDIM+1)
## n<10 matlab: -2 * ones(1,MCMCPar.n);
obj$R.stat[1,] <- c(counter.fun.evals,rep(-2,NDIM))
- if (!is.null(names(pars))) colnames(obj$R.stat) <- c("fun.evals",names(pars))
- else colnames(obj$R.stat) <- c("fun.evals",paste("p",1:length(pars),sep=""))
+ colnames(obj$R.stat) <- c("fun.evals", names(pars))
##counter.fun.evals + pCR for each CR
obj$CR <- matrix(NA,max.counter.outloop,length(pCR)+1)
@@ -239,14 +240,14 @@
obj$outlier<-NULL
- Sequences <- array(NA, c(max.counter,NDIM+2,NSEQ))
- if (!is.null(names(pars))) colnames(Sequences) <- c(names(pars),"p","logp")
+ Sequences <- array(NA_real_, c(max.counter,NDIM+2,NSEQ))
+ colnames(Sequences) <- c(names(pars), "p", "logp")
## Sequences[1,] <- sapply(pars, mean) ## TODO: include?
## Check whether will save a reduced sample
if (!is.na(control$thin.t)){
counter.redseq <- 0
- Reduced.Seq <- array(NA,c(ceiling(max.counter/control$thin.t),NDIM+2,NSEQ))
+ Reduced.Seq <- array(NA_real_,c(ceiling(max.counter/control$thin.t),NDIM+2,NSEQ))
} else Reduced.Seq <- NULL
############################
@@ -280,8 +281,8 @@
##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")
-
+ colnames(X) <- c(names(pars), "p", "logp")
+
##Initialise the sequences
for (qq in 1:NSEQ){
Sequences[1,,qq] <- X[qq,]
@@ -486,11 +487,13 @@
thin=control$thin.t)
))
}
-
+
+ ## TODO: make these 'ts' objects and sync with Reduced.Seq by thinning
obj$X <- X
obj$R.stat <- obj$R.stat[1:counter.report,,drop=FALSE]
- obj$AR <- obj$AR[1:(counter-1),]
- obj$CR <- obj$CR[1:(counter.outloop-1),]
+ obj$hist.logp <- 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]
## store number of iterations
obj$iterations <- counter.outloop
Modified: pkg/man/coef.dream.Rd
===================================================================
--- pkg/man/coef.dream.Rd 2010-05-07 08:13:54 UTC (rev 24)
+++ pkg/man/coef.dream.Rd 2010-05-11 05:15:27 UTC (rev 25)
@@ -2,7 +2,8 @@
\alias{coef.dream}
\title{Extract parameter values from dream object}
\usage{
-\S3method{coef}{dream}(object, last.prop=0.5, use.thinned=FALSE, method=c("maxLik", "mean", "median"),\dots)
+\method{coef}{dream}(object, last.prop = 0.5, use.thinned = FALSE,
+ method = c("uni.mode", "mean", "median", "sample.ml"), \dots)
}
\description{Extract parameter values using a choice of methods (or an
arbitrary function)}
@@ -12,8 +13,29 @@
\item{last.prop}{proportion of total sequence to use (0,1]
if 1, use whole sequence}
\item{use.thinned}{Use thinned MCMC chains}
- \item{method}{either one of maxLik,mean,median or a function parameter
- vector=f(dream object)}
+ \item{method}{method for extracting a parameter set from the MCMC
+ chains. One of:
+ \describe{
+ \item{\code{"uni.mode"}}{
+ mode of the univariate density estimate
+ for each parameter, using settings as in
+ \code{\link{densityplot.mcmc}}.
+ }
+ \item{\code{"mean"}}{
+ mean of each univariate parameter distribution.
+ }
+ \item{\code{"median"}}{
+ median of each univariate parameter distribution.
+ }
+ \item{\code{"sample.ml"}}{
+ parameter set with maximum likelihood (according to the chosen
+ likelihood function) from the generated MCMC chains.
+ }
+ \item{\code{function(object)}}{
+ a function of the dream object which returns a parameter vector.
+ }
+ }
+ }
\item{...}{Unused. Matches generic}
}
\details{
Modified: pkg/man/dream.Rd
===================================================================
--- pkg/man/dream.Rd 2010-05-07 08:13:54 UTC (rev 24)
+++ pkg/man/dream.Rd 2010-05-11 05:15:27 UTC (rev 25)
@@ -148,17 +148,18 @@
}
\value{
A list with elements:
- \item{X}{converged nseq points in parameter space. matrix nseq x ndim}
- \item{Sequences}{ \code{\link{mcmc.list}}. nseq mcmc elements of ndim variables}
- \item{Reduced.Seq}{ \code{\link{mcmc.list}}. nseq mcmc elements of ndim
- variables, if \code{control$thin.t!=NA}}
- \item{AR}{Acceptance rate for each draw. matrix n.elem x 2}
+ \item{X}{converged nseq points in parameter space. matrix \code{nseq x ndim}.}
+ \item{Sequences}{ \code{\link{mcmc.list}}. \code{nseq} mcmc elements of \code{ndim} variables.}
+ \item{Reduced.Seq}{ \code{\link{mcmc.list}}. \code{nseq} mcmc elements of \code{ndim}
+ variables, if \code{control$thin.t!=NA}.}
+ \item{hist.logp}{History of log(\var{p}) values. matrix \code{nseq x n.elem}.}
+ \item{AR}{Acceptance rate for each draw. matrix \code{n.elem x 2}.}
\item{outlier}{Iterations at which chains were replaced. vector of variable length}
- \item{R.stat}{gelman.diag statistic for each variable at each step. matrix n.elem/steps x 1+ndim}
- \item{CR}{Probability of crossover at each step. matrix n.elem/steps x
- 1+length(pCR)}
+ \item{R.stat}{gelman.diag statistic for each variable at each step. matrix \code{n.elem/steps x 1+ndim}.}
+ \item{CR}{Probability of crossover at each step. matrix \code{n.elem/steps x
+ 1+length(pCR)}.}
\item{in.burnin}{Boolean showing whether dream was in the burn-in
- period at termination (see description above)}
+ period at termination (see description above).}
\item{fun.evals}{Total number of function evaluations.}
\item{time}{running time (wall time) in seconds.}
\item{EXITMSG}{a message about the stopping condition.}
More information about the Dream-commits
mailing list