[Dream-commits] r23 - in pkg: R demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Apr 28 07:09:21 CEST 2010
Author: josephguillaume
Date: 2010-04-28 07:09:21 +0200 (Wed, 28 Apr 2010)
New Revision: 23
Added:
pkg/R/maxLikCoda.R
pkg/man/plot.dream.Rd
Modified:
pkg/R/CompDensity.R
pkg/R/coef.dream.R
pkg/R/predict.dream.R
pkg/demo/FME.nonlinear.model.R
pkg/demo/FME.nonlinear.model_parallelisation.R
pkg/man/dream.Rd
Log:
Changed multicore error message
Added maxLikCoda, made default for coef - matches par prob plots and is therefore the intuitive result.
Allowed output from newFUN to be a single value for "CI"
Added help file plot.dream
Modified: pkg/R/CompDensity.R
===================================================================
--- pkg/R/CompDensity.R 2010-03-19 02:30:25 UTC (rev 22)
+++ pkg/R/CompDensity.R 2010-04-28 05:09:21 UTC (rev 23)
@@ -112,7 +112,10 @@
p <- sapply(temp,function(x) x[1])
logp <- sapply(temp,function(x) x[2])
- if(class(p)!="numeric") stop(sprintf("Expected class numeric, got class %s. Error with multicore? Set control$use.multicore=FALSE to turn it off",class(p)))
+ if(class(p)!="numeric") {
+ print(p)
+ stop(sprintf("Expected class numeric, got class %s. Error with multicore? Set control$parallel='none' to not use parallelisation",class(p)))
+ }
stopifnot(!any(is.na(p)))
##stopifnot(!any(is.na(logp))) ##Not used anyway
return(list(p=p,logp=logp))
Modified: pkg/R/coef.dream.R
===================================================================
--- pkg/R/coef.dream.R 2010-03-19 02:30:25 UTC (rev 22)
+++ pkg/R/coef.dream.R 2010-04-28 05:09:21 UTC (rev 23)
@@ -22,7 +22,7 @@
if (class(method)!="function") {
method <- switch(
match.arg(method),
- "maxLik"=maxLikPars,
+ "maxLik"=maxLikCoda,
"mean"=function(sss) colMeans(as.matrix(sss)),
"median"=function(sss) apply(as.matrix(sss),2,median)
)
Added: pkg/R/maxLikCoda.R
===================================================================
--- pkg/R/maxLikCoda.R (rev 0)
+++ pkg/R/maxLikCoda.R 2010-04-28 05:09:21 UTC (rev 23)
@@ -0,0 +1,40 @@
+maxLikCoda <- function (x)
+{
+ xx <- as.matrix(x)
+ pars.maxp <- rep(NA, ncol(xx))
+ for (i in 1:nvar(x)) {
+ y <- xx[, i, drop = TRUE]
+ bwf <- function(x) {
+ x <- x[!is.na(as.vector(x))]
+ return(1.06 * min(sd(x), IQR(x)/1.34) * length(x)^-0.2)
+ }
+ bw <- bwf(y)
+ width <- 4 * bw
+ scale <- "open"
+ if (max(y) <= 1 && 1 - max(y) < 2 * bw) {
+ if (min(y) >= 0 && min(y) < 2 * bw) {
+ scale <- "proportion"
+ y <- c(y, -y, 2 - y)
+ }
+ }
+ else if (min(y) >= 0 && min(y) < 2 * bw) {
+ scale <- "positive"
+ y <- c(y, -y)
+ }
+ else scale <- "open"
+ dens <- density(y, width = width)
+ if (scale == "proportion") {
+ dens$y <- 3 * dens$y[dens$x >= 0 & dens$x <=
+ 1]
+ dens$x <- dens$x[dens$x >= 0 & dens$x <= 1]
+ }
+ else if (scale == "positive") {
+ dens$y <- 2 * dens$y[dens$x >= 0]
+ dens$x <- dens$x[dens$x >= 0]
+ }
+ ii <- which.max(dens$y)
+ pars.maxp[i] <- dens$x[ii]
+
+ }##for
+ return(pars.maxp)
+}##maxLikCoda
Property changes on: pkg/R/maxLikCoda.R
___________________________________________________________________
Name: svn:executable
+ *
Modified: pkg/R/predict.dream.R
===================================================================
--- pkg/R/predict.dream.R 2010-03-19 02:30:25 UTC (rev 22)
+++ pkg/R/predict.dream.R 2010-04-28 05:09:21 UTC (rev 23)
@@ -77,8 +77,9 @@
if (last.prop<1) sss <- window(sss, start = end(sss)*(1-last.prop) + 1)
ff <- apply(as.matrix(sss),1,wrap)
-
+
if (inherits(ff,"matrix")) return(t(apply(ff,1,quantile,c((1-level)/2,1-(1-level)/2))))
+ else if (inherits(ff,"numeric")) return(quantile(ff,c((1-level)/2,1-(1-level)/2)))
else if (inherits(ff,"list")) {
## Calculate CI for each series separately
## list is of format list[[run.number]][[series.number]]=numeric
Modified: pkg/demo/FME.nonlinear.model.R
===================================================================
--- pkg/demo/FME.nonlinear.model.R 2010-03-19 02:30:25 UTC (rev 22)
+++ pkg/demo/FME.nonlinear.model.R 2010-04-28 05:09:21 UTC (rev 23)
@@ -50,6 +50,7 @@
control <- list(
nseq=4,
+ thin.t=10,
parallel="none"
## REPORT=0
## ndraw=1000
Modified: pkg/demo/FME.nonlinear.model_parallelisation.R
===================================================================
--- pkg/demo/FME.nonlinear.model_parallelisation.R 2010-03-19 02:30:25 UTC (rev 22)
+++ pkg/demo/FME.nonlinear.model_parallelisation.R 2010-04-28 05:09:21 UTC (rev 23)
@@ -1,38 +1,35 @@
+## Example of parallelisation
+## Uses FME data (see other example)
-## Example from FME vignette, p21, section 7. Soetaert & Laine
-## Nonlinear model parameter estimation
-## TODO: document more clearly, better outputs
-##
-
-## http://r-forge.r-project.org/plugins/scmsvn/viewcvs.php/pkg/FME/inst/doc/FMEmcmc.Rnw?rev=96&root=fme&view=markup
-
-Obs <- data.frame(x=c( 28, 55, 83, 110, 138, 225, 375), # mg COD/l
+obs.all <- data.frame(x=c( 28, 55, 83, 110, 138, 225, 375), # mg COD/l
y=c(0.053,0.06,0.112,0.105,0.099,0.122,0.125)) # 1/hour
-Obs2<- data.frame(x=c( 20, 55, 83, 110, 138, 240, 325), # mg COD/l
- y=c(0.05,0.07,0.09,0.10,0.11,0.122,0.125)) # 1/hour
-obs.all <- rbind(Obs,Obs2)
-##########################
-## DREAM results
unloadNamespace("dream")
library(dream)
-Model.y <- function(p,x) p[1]*x/(x+p[2])
-
control <- list(
nseq=4
)
-
pars <- list(p1=c(0,1),p2=c(0,100))
+
+##########################
+## Demonstration of parallelisation packages
+## Comparison of run times for known fixed termination problem
+
+Model.y <- function(p,x) {
+ p[1]*x/(x+p[2])
+
+}
+
+system.time(sapply(1:5e4,function(x) Model.y(c(0,0),obs.all$x)))[["elapsed"]]/5e4
+
## Optional parallelisation of function evaluations
## Uses foreach, which allows Rmpi, SNOW or multicore to be registered,
## or just normal sequential evaluation
## Example here is for SNOW with a socket cluster, which doesn't require any other software or setup
## Parallelisation incurs an overhead and is not worthwhile for fast functions
-
-
if (require(doSNOW)){
cl <- makeCluster(2, type = "SOCK")
registerDoSNOW(cl)
@@ -64,3 +61,65 @@
## none: 6.234, 6.219
## snow: 12.234, 12.125
## foreach: 107.999, 107.893
+
+## For function with 1e-3 sleep
+## [1] "none"
+## p1 p2
+## 0.1511961 50.7765344
+## [1] 109.483
+## [1] "snow"
+## p1 p2
+## 0.1511961 50.7765344
+## [1] 54.858
+## [1] "foreach"
+## p1 p2
+## 0.1511961 50.7765344
+## [1] 109.593
+
+
+##################
+## Test of number of function evaluations in fixed time for a time-expensive function
+
+Model.y <- function(p,x) {
+ Sys.sleep(1e-3)
+ p[1]*x/(x+p[2])
+
+}
+
+system.time(sapply(1:5e1,function(x) Model.y(c(0,0),obs.all$x)))[["elapsed"]]/5e1
+
+
+if (require(doSNOW)){
+ cl <- makeCluster(2, type = "SOCK")
+ registerDoSNOW(cl)
+}
+
+for (p in c("none","snow","foreach")){
+ set.seed(456)
+ control$parallel <- p
+ control$maxtime <- 20
+
+ dd <- dream(
+ FUN=Model.y, func.type="calc.rmse",
+ pars = pars,
+ FUN.pars=list(
+ x=obs.all$x
+ ),
+ INIT = LHSInit,
+ measurement=list(data=obs.all$y),
+ control = control
+ )
+ print(dd$control$parallel)
+ print(dd$fun.evals)
+}
+
+if (require(doSNOW)) stopCluster(cl)
+
+## Number of function evaluations
+## TODO: appears to be an error in foreach evaluation
+## [1] "none"
+## [1] 1280
+## [1] "snow"
+## [1] 2560
+## [1] "foreach"
+## [1] 1280
Modified: pkg/man/dream.Rd
===================================================================
--- pkg/man/dream.Rd 2010-03-19 02:30:25 UTC (rev 22)
+++ pkg/man/dream.Rd 2010-04-28 05:09:21 UTC (rev 23)
@@ -90,12 +90,15 @@
If dream frequently warns that it is reentering the burn-period,
\code{control$burnin.length} may need to be increased.
- There are S3 methods for print, summary, plot,
+ There are S3 methods for print, summary,
+ \code{\link[=plot.dream]{plot}},
\code{\link[=coef.dream]{coef}},
\code{\link[=predict.dream]{predict}},
simulate
- coef uses \code{\link{maxLikPars}}, a naive approach to extracting the most likely parameter value.
+ coef uses \code{\link{maxLikPars}}, a naive approach to extracting the
+ most likely parameter value, using the bandwith parameters used for
+ CODA density plots.
}
\value{
Added: pkg/man/plot.dream.Rd
===================================================================
--- pkg/man/plot.dream.Rd (rev 0)
+++ pkg/man/plot.dream.Rd 2010-04-28 05:09:21 UTC (rev 23)
@@ -0,0 +1,58 @@
+\name{plot.dream}
+\alias{plot.dream}
+\title{Visualising DREAM results}
+\usage{
+\S3method{plot}{dream}(x,interactive=T,...)
+}
+\arguments{
+ \item{x}{A \code{\link{dream}} object}
+ \item{interactive}{whether to use devAskNewPage (the default) or print
+ all at once}
+ \item{...}{Placeholder for S3method}
+}
+\description{
+ Plots various visualisation of MCMC chain results
+}
+\details{
+ Extracts the second half of the MCMC chains as a coda object
+
+ Calls:
+ \itemize{
+ \item coda's \code{\link[=plot.mcmc]{plot}}
+ \item coda's \code{\link[=xyplot.mcmc]{xyplot}}
+ \item coda's \code{\link[=densityplot.mcmc]{densityplot}}
+ \item Plots the distribution of acceptance rates
+ \item coda's \code{\link{gelman.plot}} (may fail depending on characteristics of
+ MCMC chains)
+ \item Plots the evolution of the Gelman R statistic, using previously
+ internally compueted results
+ \item Plots the multi-variate density of the first chain
+ }
+}
+\examples{
+\dontrun{
+ ##To run the equivalent coda visualisations:
+
+ ss <- window(dream.obj$Sequences, start = end(dream.obj$Sequences)/2 + 1)
+ plot(ss)
+ print(xyplot(ss))
+
+ ## Lattice graphics parameter density plots
+ ## Each chain is superposed
+ print(densityplot(ss))
+
+ ## All chains combined (treated as one)
+ print(densityplot(as.mcmc(as.matrix(ss))))
+
+ ## Alternatively, specify burn-in within the call
+ print(densityplot(dream.obj$Sequences, start = end(dream.obj$Sequences)/2 + 1))
+
+ ## Base graphics parameter density plots
+ plot(ss,trace=F)
+
+ ## Uses densplot behind the scenes, e.g.
+ densplot(ss[,1])
+
+ }
+}
+
Property changes on: pkg/man/plot.dream.Rd
___________________________________________________________________
Name: svn:executable
+ *
More information about the Dream-commits
mailing list