[Dream-commits] r5 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Feb 15 07:48:14 CET 2010
Author: josephguillaume
Date: 2010-02-15 07:48:13 +0100 (Mon, 15 Feb 2010)
New Revision: 5
Added:
pkg/R/rand_utils.R
Modified:
pkg/R/CompDensity.R
pkg/R/CovInit.R
pkg/R/RemOutlierChains.R
pkg/R/dream.R
pkg/R/metrop.R
Log:
Fixed R CMD check fails. Added dream parameter func.type to specify type of function output
Modified: pkg/R/CompDensity.R
===================================================================
--- pkg/R/CompDensity.R 2010-02-14 23:42:16 UTC (rev 4)
+++ pkg/R/CompDensity.R 2010-02-15 06:48:13 UTC (rev 5)
@@ -1,32 +1,38 @@
##' Computes the density of each x value
##'
##' @param x matrix nseq x ndim
-##' @param control. list containing gamma,Wb,Cb
+##' @param control. list containing gamma,Wb,Cb,nseq
##' @param FUN - the model to run
##' R function with first argument a vector of length ndim.
-##' returning a single scalar corresponding to one of the options:
-##' @param option unambiguous abbrev of:
-##' posterior.density, calc.loglik, calc.rmse, logposterior.density,calc.weighted.rmse
+##' returns a scalar or vector corresponding to one of the options below.
+##' @param option 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 measurement list containing TODO: not sure
##' data: vector of observations corresponding to model output
##' sigma: scalar
##' @param ... additional arguments to FUN
-##' @return ... list with components
+##' @return list with elements
##' p vector of length nseq
##' logp vector of length nseq
##
## TODO: p may be erroneously equal to logp?
-CompDensity <- function(x,control,FUN,option,
- measurement,...){
+## TODO: more appropriate naming of options?
+## TODO: allow shortenings of option?
+CompDensity <- function(x,control,FUN,func.type,
+ measurement=NULL,...){
+ stopifnot(!is.null(measurement) || func.type%in% c("posterior.density","logposterior.density"))
+
## dimensions:
## i. iter 1:nseq
## modpred. scalar or vector commensurate to measurement$data
## err. vector of same length as modpred
## SSR scalar
- p <- rep(NA,nseq)
- logp <- rep(NA,nseq)
+ p <- rep(NA,control$nseq)
+ logp <- rep(NA,control$nseq)
## Sequential evaluation
for (ii in 1:nrow(x)){
@@ -34,7 +40,7 @@
## TODO: correct use of optional pars?
modpred <- FUN(x[ii,],...)
- switch(option,
+ switch(func.type,
## Model directly computes posterior density
posterior.density={
p[ii] <- modpred
Modified: pkg/R/CovInit.R
===================================================================
--- pkg/R/CovInit.R 2010-02-14 23:42:16 UTC (rev 4)
+++ pkg/R/CovInit.R 2010-02-15 06:48:13 UTC (rev 5)
@@ -17,7 +17,7 @@
{
##[x] = repmat(Extra.muX,MCMCPar.seq,1) + randn(MCMCPar.seq,MCMCPar.n) * chol(Extra.qcov);
- x <- t(matrix(rep(a,nseq),length(a)))+randn(nseq,length(pars)) %*% chol(qcov)
+ x <- t(matrix(rep(muX,nseq),length(muX)))+randn(nseq,length(pars)) %*% chol(qcov)
lower <- sapply(pars, function(x) min(x[[1]]))
upper <- sapply(pars, function(x) max(x[[1]]))
Modified: pkg/R/RemOutlierChains.R
===================================================================
--- pkg/R/RemOutlierChains.R 2010-02-14 23:42:16 UTC (rev 4)
+++ pkg/R/RemOutlierChains.R 2010-02-15 06:48:13 UTC (rev 5)
@@ -40,7 +40,7 @@
## 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]
+ iqr <- q13[1]-q13[2]
## Compute the upper range -- to detect outliers
## TODO: shouldn't upper.range be Q3+3*IQR?
upper.range <- q13[2]-2*iqr
Modified: pkg/R/dream.R
===================================================================
--- pkg/R/dream.R 2010-02-14 23:42:16 UTC (rev 4)
+++ pkg/R/dream.R 2010-02-15 06:48:13 UTC (rev 5)
@@ -25,11 +25,13 @@
## 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)),
- measurement,
+dream <- function(FUN, func.type,
+ pars = list(x = range(0, 1e6)),
...,
INIT = LHSInit,
- control = list()) #, do.SSE = FALSE)
+ control = list(),
+ measurement=NULL
+ )
{
@@ -44,12 +46,14 @@
## delta.tot vector of length nCR
- If (is.character(FUN))
- FUN <- get(FUN, mode = "function")
+ if (is.character(FUN))
+ FUN <- get(FUN, mode = "function")
stopifnot(is.function(FUN))
stopifnot(is.list(par))
stopifnot(length(par) > 0)
+ stopifnot(!is.null(measurement) || func.type %in% c("posterior.density","logposterior.density"))
+
pars <- lapply(pars, function(x) if (is.list(x)) x else list(x))
############################
@@ -144,8 +148,8 @@
## Sequences[1,] <- sapply(pars, mean) ## TODO: include?
## Check whether will save a reduced sample
- if (thin){
- Reduced.Seq <- array(NA,c(floor(n.elem/thin.t),NDIM+2,NSEQ))
+ if (control$thin){
+ Reduced.Seq <- array(NA,c(floor(n.elem/control$thin.t),NDIM+2,NSEQ))
} else Reduced.Seq <- NULL
############################
@@ -160,14 +164,14 @@
################################
## Step 1: Sample s points in the parameter space
- x <- INIT(pars, nseq,...)
+ x <- INIT(pars, 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, ...)
+ tmp<-CompDensity(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")
@@ -187,7 +191,7 @@
##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])))))
+ obj$R.stat[1,] <- c(Iter,gelman.diag(as.mcmc.list(lapply(1:NSEQ,function(i) as.mcmc(Sequences[1:iloc,1:NDIM,i])))))
################################
##Start iteration
@@ -234,7 +238,7 @@
Sequences[iloc,,] <- t(newgen)
## Check reduced sample collection
- if (thin && new_teller == thin.t){
+ if (control$thin && new_teller == control$thin.t){
## Update iloc_2 and new_teller
iloc.2 <- iloc.2+1
new_teller <- 0
@@ -259,7 +263,7 @@
hist.logp[counter,] <- X[,NDIM+2]
## Save some important output -- Acceptance Rate
- obj$AR[counter] <- 100 * sum(accept) / nseq
+ obj$AR[counter] <- 100 * sum(accept) / NSEQ
## Update Iteration and counter
Iter <- Iter + NSEQ
@@ -309,7 +313,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[,1:NDIM,i]))),autoburnin=TRUE)
## break if maximum time exceeded
toc <- as.numeric(Sys.time()) - tic
@@ -331,7 +335,7 @@
obj$Sequences <- Sequences[1:i,,]
}
## Remove extra rows from Reduced.Seq
- if (thin){
+ if (control$thin){
i <- which(rowSums(Reduced.Seq)==0)
if (length(i)>0) {
i <- i[1]-1
@@ -358,7 +362,8 @@
}
## store number of function evaluations
- obj$counts <- funevals
+ ## TODO: funevals is not calculated atm
+ ## obj$counts <- funevals
## store number of iterations
obj$iterations <- i
## store the amount of time taken
Modified: pkg/R/metrop.R
===================================================================
--- pkg/R/metrop.R 2010-02-14 23:42:16 UTC (rev 4)
+++ pkg/R/metrop.R 2010-02-15 06:48:13 UTC (rev 5)
@@ -7,13 +7,14 @@
##' @param p.old vector of length nseq
##' @param logp.old vector of length nseq
##' @param measurement. needs N and sigma
-##' @param control. needs gamma,metrop.opt range [1,5]
-
+##' @param control list with elements:
+##' gamma
+##' metrop.opt range [1,5]
##' @return ... list with elements
##' 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,
measurement,control
@@ -34,7 +35,7 @@
## And initialize accept with false
accept <- rep(FALSE,nr.chains)
- switch(option,
+ switch(control$metrop.opt,
"1" = {
alpha <- min(p.x/p.old,1)
},
Added: pkg/R/rand_utils.R
===================================================================
--- pkg/R/rand_utils.R (rev 0)
+++ pkg/R/rand_utils.R 2010-02-15 06:48:13 UTC (rev 5)
@@ -0,0 +1,2 @@
+rand<-function(m,n) matrix(runif(m*n),m,n)
+randn<-function(m,n) matrix(rnorm(m*n),m,n)
\ No newline at end of file
Property changes on: pkg/R/rand_utils.R
___________________________________________________________________
Name: svn:executable
+ *
More information about the Dream-commits
mailing list