[Dream-commits] r29 - in pkg: . R demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu May 13 05:34:02 CEST 2010
Author: josephguillaume
Date: 2010-05-13 05:34:01 +0200 (Thu, 13 May 2010)
New Revision: 29
Added:
pkg/R/simulate.dream.R
pkg/R/window.dream.R
pkg/man/dreamCalibrate.Rd
pkg/man/simulate.dream.Rd
pkg/man/window.dream.Rd
Removed:
pkg/R/fitted.dream.R
pkg/R/possibility.envelope.R
pkg/man/maxLikPars.Rd
pkg/man/possibility.envelope.Rd
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/coef.dream.R
pkg/R/dream.R
pkg/R/dreamCalibrate.R
pkg/R/plot.dream.R
pkg/R/predict.dream.R
pkg/R/summary.dream.R
pkg/demo/FME.nonlinear.model.R
pkg/demo/FME.nonlinear.model_parallelisation.R
pkg/demo/run_example1.R
pkg/man/coef.dream.Rd
pkg/man/dream.Rd
pkg/man/predict.dream.Rd
Log:
- Working dreamCalibrate and likelihood functions.
- predict now operates on dream_model object
- If thinning is used, Reduced.Seq is returned INSTEAD of Sequences. Made necessary modifications to other functions
- Renamed fitted to window.dream
- coef now uses window.dream
- removed possibility.envelope. More advanced functionality now in simulate and predict
- removed documentation for possibility.envelope, maxLikPars
- added documentation for dreamCalibrate and placeholders for window.dream, simulate.dream
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-05-12 05:58:05 UTC (rev 28)
+++ pkg/DESCRIPTION 2010-05-13 03:34:01 UTC (rev 29)
@@ -1,14 +1,14 @@
Package: dream
-Version: 0.2-1
-Date: 2010-05-11
+Version: 0.3-1
+Date: 2010-05-13
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>
Description: Efficient global MCMC even in high-dimensional spaces.
Based on the original MATLAB code written by Jasper Vrugt.
+ Methods for calibration and prediction using the DREAM algorithm
Depends: coda
-Suggests: multicore, foreach
+Suggests: multicore, foreach, SNOW, doSNOW, doMPI, doMC
Imports: stats, utils
-Enhances: doSNOW, doMPI, doMC
License: file LICENSE
URL: http://dream.r-forge.r-project.org/
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2010-05-12 05:58:05 UTC (rev 28)
+++ pkg/NAMESPACE 2010-05-13 03:34:01 UTC (rev 29)
@@ -4,12 +4,19 @@
export(dream)
export(CovInit)
export(LHSInit)
-export(maxLikPars)
-export(possibility.envelope)
+
S3method(summary,dream)
+S3method(plot,dream)
+S3method(window,dream)
S3method(coef,dream)
-S3method(fitted,dream)
-S3method(plot,dream)
S3method(print,dream)
-S3method(predict,dream)
S3method(simulate,dream)
+
+export(maxLikPars)
+export(maxLikCoda)
+
+export(dreamCalibrate)
+export(calc.rmse)
+export(calc.weighted.rmse)
+export(calc.loglik)
+S3method(predict,dream_model)
Modified: pkg/R/coef.dream.R
===================================================================
--- pkg/R/coef.dream.R 2010-05-12 05:58:05 UTC (rev 28)
+++ pkg/R/coef.dream.R 2010-05-13 03:34:01 UTC (rev 29)
@@ -1,25 +1,13 @@
##' Extract maximum likelihood parameter values
##' @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 uni.mode,mean,median,sample.ml
+##' @param ... arguments to window.dream
##' @return named vector of parameter values
-coef.dream <- function(object,last.prop=.5,use.thinned=FALSE,
- method=c("uni.mode","mean","median","sample.ml"),...)
+coef.dream <- function(object,method=c("uni.mode","mean","median","sample.ml"),...)
{
+
+ sss <- window(object,...)
- stopifnot(last.prop>0)
-
- if (use.thinned & is.null(object$Reduced.Seq)) {
- warning("Attempted to use.thinned when no thinned chains available: setting use.thinned=FALSE")
- use.thinned <- FALSE
- }
-
- if (use.thinned) sss <- object$Reduced.Seq
- else sss <- object$Sequences
-
- stopifnot(!is.null(sss))
-
if (identical(method, "sample.ml")) {
## TODO: make sure ppp corresponds to sss
ppp <- object$hist.logp
@@ -38,9 +26,5 @@
)
}
- if (last.prop==1) return(method(sss))
- else {
- ss <- window(sss, start = end(sss)*(1-last.prop) + 1)
- return(method(ss))
- }
+ return(method(sss))
}##coef.dream
Modified: pkg/R/dream.R
===================================================================
--- pkg/R/dream.R 2010-05-12 05:58:05 UTC (rev 28)
+++ pkg/R/dream.R 2010-05-13 03:34:01 UTC (rev 29)
@@ -144,20 +144,19 @@
## Check INIT and FUN have required extra parameters in INIT.pars & FUN.pars
req.args.init <- names(formals(INIT))
+ if(!all(req.args.init %in% c("pars","nseq",names(INIT.pars))))
+ stop(paste(c("INIT Missing extra arguments:",
+ req.args.init[!req.args.init %in% c("pars","nseq",names(INIT.pars))]),
+ sep=" ",collapse=" "))
+
req.args.FUN <- names(formals(FUN))
-
- if(!all(req.args.init %in% c("pars","nseq",names(INIT.pars))))
- stop(paste(c("INIT Missing extra arguments:",
- req.args.init[!req.args.init %in% c("pars","nseq",names(INIT.pars))]),
- sep=" ",collapse=" "))
-
- if (length(req.args.FUN)<length(FUN.pars)+1) stop("Some FUN.pars are not required by FUN")
- if (length(req.args.FUN)>1){
- req.args.FUN <- req.args.FUN[2:length(req.args.FUN)] ##optional pars only
- if(!all(req.args.FUN %in% names(FUN.pars))) stop(paste(c("FUN Missing extra arguments:",
- req.args.FUN[!req.args.FUN %in% c(names(FUN.pars))]),
- collapse=" "))
- }
+ ## if (length(req.args.FUN)<length(FUN.pars)+1) stop("Some FUN.pars are not required by FUN")
+ ## if (length(req.args.FUN)>1){
+ ## req.args.FUN <- req.args.FUN[2:length(req.args.FUN)] ##optional pars only
+ ## if(!all(req.args.FUN %in% names(FUN.pars))) stop(paste(c("FUN Missing extra arguments:",
+ ## req.args.FUN[!req.args.FUN %in% c(names(FUN.pars))]),
+ ## collapse=" "))
+ ## }
## Update default settings with supplied settings
@@ -339,7 +338,6 @@
##Save history log density of individual chains
hist.logp[1,] <- X[,"logp"]
-
################################
##Start iteration
@@ -481,12 +479,16 @@
counter.report <- counter.report+1
try({
- obj$R.stat[counter.report,] <- c(counter.fun.evals,gelman.diag(
- as.mcmc.list(lapply(1:NSEQ,function(i) as.mcmc(Sequences[1:(counter-1),1:NDIM,i]))),
- autoburnin=TRUE)$psrf[,1])
- if (counter.report == 1)
- message(format(colnames(obj$R.stat), width = 5))
- message(format(obj$R.stat[counter.report,], width = 5, digits = 4))
+ obj$R.stat[counter.report,] <-
+ c(counter.fun.evals,
+ gelman.diag(
+ as.mcmc.list(lapply(1:NSEQ,function(i) as.mcmc(Sequences[1:(counter-1),1:NDIM,i]))),
+ autoburnin=TRUE)$psrf[,1])
+ if (counter.report == 2){
+ message("R.stats:")
+ message(format(colnames(obj$R.stat), width = 10,justify="right"))
+ }
+ message(format(obj$R.stat[counter.report,], width = 10, digits = 4))
})
if (all(!is.na(obj$R.stat[counter.report,])) &&
@@ -496,7 +498,7 @@
## obj$EXITFLAG <- 3
}
- }##counter.report
+ } ##counter.report
## Update the counter.outloop
counter.outloop = counter.outloop + 1
@@ -522,14 +524,15 @@
## Trim outputs to collected data - remove extra rows
## Convert sequences to mcmc objects
- Sequences <- Sequences[1:(counter-1),,]
- obj$Sequences <- as.mcmc.list(lapply(1:NSEQ,function(i) as.mcmc(Sequences[,1:NDIM,i])))
if (!is.na(control$thin.t)){
Reduced.Seq <- Reduced.Seq[1:counter.redseq,,]
- obj$Reduced.Seq <- as.mcmc.list(lapply(1:NSEQ, function(i) {
- mcmc(Reduced.Seq[,1:NDIM,i], start = 1,
- end = counter-1, thin = control$thin.t)
- }))
+ obj$Sequences <- as.mcmc.list(lapply(1:NSEQ, function(i) {
+ mcmc(Reduced.Seq[,1:NDIM,i], start = 1,
+ end = counter-1, thin = control$thin.t)
+ }))
+ } else {
+ Sequences <- Sequences[1:(counter-1),,]
+ obj$Sequences <- as.mcmc.list(lapply(1:NSEQ,function(i) as.mcmc(Sequences[,1:NDIM,i])))
}
## TODO: make these 'ts' objects and sync with Reduced.Seq by thinning
Modified: pkg/R/dreamCalibrate.R
===================================================================
--- pkg/R/dreamCalibrate.R 2010-05-12 05:58:05 UTC (rev 28)
+++ pkg/R/dreamCalibrate.R 2010-05-13 03:34:01 UTC (rev 29)
@@ -1,26 +1,40 @@
## control requires gamma, (N)
-calc.rmse <- function(predicted,observed,control){
- ## TODO: more general case with multidim data?
+calc.rmse <- function(predicted,observed,control=list(gamma=0)){
if (! "N" %in% names(control)) control$N <- length(observed)
+ req.control <- c("N","gamma")
+ if (!all(req.control %in% names(control)))
+ stop(sprintf("Missing arguments to control: %s",
+ paste(req.control[!req.control %in% names(control)],collapse=", ")
+ ))
err <- as.numeric(observed-predicted)
- ## Derive the sum of squared error
SSR <- sum(abs(err)^(2/(1+control$gamma)))
-
- ## TODO: Correctness: is this a p or logp?
- p <- -SSR
- logp <- -0.5*SSR
-
- ## Computational issues
- (-control$N*(1+control$gamma)/2)*log(p)
+ (-control$N*(1+control$gamma)/2)*0.5*SSR
}## calc.rmse
-## TODO: calc.weighted.rmse
+calc.weighted.rmse <- function(predicted,observed,control=list(gamma=0)){
+ ##if (! "sigma" %in% names(control)) control$sigma <- sd(observed)
+ req.control <- c("gamma","sigma")
+ if (!all(req.control %in% names(control)))
+ stop(sprintf("Missing arguments to control: %s",
+ paste(req.control[!req.control %in% names(control)],collapse=", ")
+ ))
+
+ err <- as.numeric(observed-predicted)
+ SSR <- sum(abs(err)^(2/(1+control$gamma)))
+ -0.5*SSR/control$sigma^2
+}
## control: Wb,Cb,gamma,measurement.sigma
-calc.loglik <- function(predicted,observed,control){
+calc.loglik <- function(predicted,observed,control=list(gamma=0)){
+ req.control <- c("gamma")
+ if (!all(req.control %in% names(control)))
+ stop(sprintf("Missing arguments to control: %s",
+ paste(req.control[!req.control %in% names(control)],collapse=", ")
+ ))
+
if (!all(c("Cb","Wb") %in% names(control))) control <- modifyList(control,CalcCbWb(control$gamma))
if (! "sigma" %in% names(control)) control$sigma <- sd(observed)
if (! "N" %in% names(control)) control$N <- length(observed)
@@ -28,29 +42,32 @@
err <- as.numeric(observed-predicted)
## Derive the log likelihood
- logp <- control$N*log(control$Wb/control$sigma)-
- control$Cb*(sum((abs(err/control$sigma))^(2/(1+control$gamma))))
- ## And retain in memory
- p <- logp
+ with(control,N*log(Wb/sigma)-Cb*(sum((abs(err/sigma))^(2/(1+gamma)))))
}## calc.loglik
## Design decisions:
-## lik.fun is log.likelihood=f(predicted,observed,control)
-dreamCalibrate <- function(fun,
+## lik.fun is log.likelihood=f(predicted,observed,control=default.list)
+## must have default, even if it is an empty list
+dreamCalibrate <- function(FUN,
pars,
obs,
- lik.fun=calc.rmse,
- lik.control=list(),
+ lik.fun=calc.loglik,
+ lik.control=NULL,
+ FUN.pars=list(),
... ##Extra arguments to dream
){
-
- FUN <- function(pars,...) lik.fun(fun(pars,...),obs,lik.control)
- dd <- dream(FUN=fun,
+ if (is.null(lik.control)) lik.control <- eval(formals(lik.fun)$control)
+ wrap.lik.fun <- function(pars,...) lik.fun(FUN(pars,...),obs,lik.control)
+ dd <- dream(FUN=wrap.lik.fun,
pars=pars,
func.type="logposterior.density",
- ... ##INIT,control,FUN.pars
+ FUN.pars=FUN.pars,
+ ... ##INIT,control
)
-
- class(dd) <- c("dream-model",class(dd))
+
+ dd$call <- match.call()
+ dd$FUN <- FUN
+ dd$FUN.pars <- FUN.pars
+ class(dd) <- c("dream_model",class(dd))
dd
} ##dreamCalibrate
Deleted: pkg/R/fitted.dream.R
===================================================================
--- pkg/R/fitted.dream.R 2010-05-12 05:58:05 UTC (rev 28)
+++ pkg/R/fitted.dream.R 2010-05-13 03:34:01 UTC (rev 29)
@@ -1,10 +0,0 @@
-
-
-fitted.dream <-
- function(object,
- start = 1+(end(object$Sequences)-1)*(1-fraction),
- fraction = 0.5,
- ...)
-{
- window(object$Sequences, start = start, ...)
-}
Modified: pkg/R/plot.dream.R
===================================================================
--- pkg/R/plot.dream.R 2010-05-12 05:58:05 UTC (rev 28)
+++ pkg/R/plot.dream.R 2010-05-13 03:34:01 UTC (rev 29)
@@ -8,7 +8,7 @@
opar <- devAskNewPage(interactive)
on.exit(devAskNewPage(opar))
- ss <- fitted(x, ...)
+ ss <- window(x, ...)
## Convergence (Gelman plot)
Deleted: pkg/R/possibility.envelope.R
===================================================================
--- pkg/R/possibility.envelope.R 2010-05-12 05:58:05 UTC (rev 28)
+++ pkg/R/possibility.envelope.R 2010-05-13 03:34:01 UTC (rev 29)
@@ -1,35 +0,0 @@
-
-##' Obtain prediction confidence intervals for model, around new input values
-##' @param dd dream object
-##' @param FUN.pars new par values. if missing, use those from dream object
-##' @param ndraw number of new iterations
-##' @param conf % two-sided confidence interval
-## TODO: use caching rather than a second application of FUN to parameter sets
-possibility.envelope <- function(dd,FUN.pars=NULL,ndraw=1000,conf=99){
- stopifnot(is.null(FUN.pars) || is.list(FUN.pars))
-
- ## Generate more results from converged chains
- dd$control$REPORT <- 0
- dd$control$Rthres <- 0
- dd$control$ndraw <- ndraw
- if (is.na(dd$control$thin.t)) dd$control$thin.t <- 10
- dd$call$control <- dd$control
- dd$call$INIT <- function(pars,nseq) dd$X[,1:dd$control$ndim]
-
- print(sprintf("Will require %d function evaluations",ndraw+ndraw/dd$control$thin.t))
-
- ee <- eval(dd$call)
-
- if (is.null(FUN.pars)) FUN.pars <- eval(dd$call$FUN.pars)
- par.name <- names(formals(eval(dd$call$FUN)))[1]
-
- ff <- apply(as.matrix(ee$Reduced.Seq),1,
- function(p) {
- FUN.pars[[par.name]] <- p
- do.call(eval(dd$call$FUN),FUN.pars)
- })
-
- gg <- t(apply(ff,1,quantile,c((100-conf)/200,1-(100-conf)/200)))
-
- return(gg)
-}##possibility.envelope
Modified: pkg/R/predict.dream.R
===================================================================
--- pkg/R/predict.dream.R 2010-05-12 05:58:05 UTC (rev 28)
+++ pkg/R/predict.dream.R 2010-05-13 03:34:01 UTC (rev 29)
@@ -1,70 +1,31 @@
-##' Continue an existing dream MCMC set of chains
-##' @param object. dream object
-##' @param nsim. approximate number of function evaluations. default 1000
-##' @param seed passed to set.seed before continuing
-##' @return a dream object with approximately the requested number of function evaluations
-## TODO. extra parameters to set in control?
-## TODO: does not seem to yield stationary distribution?
-simulate.dream <- function(object,nsim=1000,seed=NULL,...){
-
- ## Generate more results from converged chains
- object$control$REPORT <- 0
- object$control$Rthres <- 0
- object$control$ndraw <- nsim
- object$control$burnin.length <- 0
- if (is.na(object$control$thin.t)) object$control$thin.t <- 10
- object$call$control <- object$control
- object$call$INIT <- function(pars,nseq) object$X[,1:object$control$ndim]
-
- print(sprintf("Will require %d function evaluations",nsim))
-
- if(!is.null(seed)) set.seed(seed)
-
- ee <- eval(object$call)
- return(ee)
-
-}##simulate.dream
-
-##' Predict values using dream object
-##' Predict values using function calibrated by dream, optionally with new data,
+##' Predict values using dream-model object
+##' Predict values using function calibrated using dreamCalibrate, optionally with new data,
##' using various methods of summarising the posterior parameter and output distributions
-##' @param object dream object
+##' @param object dream-model object
##' @param newdata. new FUN.pars list. If NULL, use object's.
-##' @param newFUN. a new function to run with same arguments as the original FUN
##' @param method CI or a \code{\link{method}} of coef
##' @param level. Requested two-sided level of confidence. For CI method.
-##' @param last.prop Proportion of MCMC chains to keep
-##' @param use.thinned Whether to use existing thinned chains
+##' @param ... arguments to window.dream
##' @return whatever FUN returns (either numeric, ts or list). For CI, either a matrix with upper and lower bound or list of matrices.
-predict.dream <- function(object,newdata=NULL,newFUN=NULL,
- method="uni.mode",level=0.99,
- last.prop=0.5,use.thinned=TRUE,...
+predict.dream_model <- function(object,newdata=NULL,
+ method="uni.mode",level=0.99, ...
){
## Check and initialise parameters
stopifnot(is.null(newdata) || is.list(newdata))
stopifnot(!"CI" %in% method || !is.null(level))
- stopifnot(last.prop>0)
-
- if (use.thinned & is.null(object$Reduced.Seq)) {
- warning("Attempted to use.thinned when no thinned chains available: setting use.thinned=FALSE")
- use.thinned <- FALSE
- }
-
###
- ## Fetch function and parameters from dream object
+ ## Fetch function and parameters from dream-model object
- if (!is.null(newFUN) & !identical(formals(eval(object$call$FUN)),formals(newFUN))) stop("FUN used in dream and newFUN must have same parameters")
- if (is.null(newFUN)) newFUN <- eval(object$call$FUN)
+ FUN <- object$FUN
+ if (is.null(newdata)) newdata <- eval(object$FUN.pars)
- if (is.null(newdata)) newdata <- eval(object$call$FUN.pars)
+ par.name <- names(formals(FUN))[1]
- par.name <- names(formals(newFUN))[1]
-
wrap <- function(p) {
newdata[[par.name]] <- p
- do.call(newFUN,newdata)
+ do.call(FUN,newdata)
}
###
@@ -72,9 +33,7 @@
if (method=="CI"){
- if (use.thinned) sss <- object$Reduced.Seq
- else sss <- object$Sequences
- if (last.prop<1) sss <- window(sss, start = end(sss)*(1-last.prop) + 1)
+ sss <- window(object,...)
ff <- apply(as.matrix(sss),1,wrap)
@@ -89,12 +48,10 @@
t(apply(ff.s,1,quantile,c((1-level)/2,1-(1-level)/2)))
})
)
- } else stop("Unexpected output from application of newFUN to matrix of parameters. newFUN should return a numeric or list of numerics")
+ } else stop("Unexpected output from application of FUN to matrix of parameters. FUN should return a numeric or list of numerics")
} else {
- return(wrap(coef(object,method=method,
- use.thinned=use.thinned,last.prop=last.prop)
- ))
+ return(wrap(coef(object,method=method,...)))
}
} ##predict.dream
Added: pkg/R/simulate.dream.R
===================================================================
--- pkg/R/simulate.dream.R (rev 0)
+++ pkg/R/simulate.dream.R 2010-05-13 03:34:01 UTC (rev 29)
@@ -0,0 +1,27 @@
+##' Continue an existing dream MCMC set of chains
+##' @param object. dream object
+##' @param nsim. approximate number of function evaluations. default 1000
+##' @param seed passed to set.seed before continuing
+##' @return a dream object with approximately the requested number of function evaluations
+## TODO. extra parameters to set in control?
+## TODO: does not seem to yield stationary distribution?
+## update method?
+simulate.dream <- function(object,nsim=1000,seed=NULL,...){
+
+ ## Generate more results from converged chains
+ object$control$REPORT <- 0
+ object$control$Rthres <- 0
+ object$control$ndraw <- nsim
+ object$control$burnin.length <- 0
+ if (is.na(object$control$thin.t)) object$control$thin.t <- 10
+ object$call$control <- object$control
+ object$call$INIT <- function(pars,nseq) object$X[,1:object$control$ndim]
+
+ message(sprintf("Will require %d function evaluations",nsim))
+
+ if(!is.null(seed)) set.seed(seed)
+
+ ee <- eval(object$call)
+ return(ee)
+
+}##simulate.dream
Property changes on: pkg/R/simulate.dream.R
___________________________________________________________________
Added: svn:executable
+ *
Modified: pkg/R/summary.dream.R
===================================================================
--- pkg/R/summary.dream.R 2010-05-12 05:58:05 UTC (rev 28)
+++ pkg/R/summary.dream.R 2010-05-13 03:34:01 UTC (rev 29)
@@ -1,7 +1,7 @@
summary.dream <- function(object, fraction = 0.5, ...){
- coda.sum <- summary(fitted(object, fraction = fraction), ...)
+ coda.sum <- summary(window(object, fraction = fraction), ...)
cat(sprintf("
Exit message: %s
Added: pkg/R/window.dream.R
===================================================================
--- pkg/R/window.dream.R (rev 0)
+++ pkg/R/window.dream.R 2010-05-13 03:34:01 UTC (rev 29)
@@ -0,0 +1,10 @@
+
+
+window.dream <-
+ function(object,
+ start = 1+(end(object$Sequences)-1)*(1-fraction),
+ fraction = 0.5,
+ ...)
+{
+ window(object$Sequences, start = start, ...)
+}
Modified: pkg/demo/FME.nonlinear.model.R
===================================================================
--- pkg/demo/FME.nonlinear.model.R 2010-05-12 05:58:05 UTC (rev 28)
+++ pkg/demo/FME.nonlinear.model.R 2010-05-13 03:34:01 UTC (rev 29)
@@ -45,33 +45,31 @@
library(dream)
Model.y <- function(p,x) p[1]*x/(x+p[2])
+pars <- list(p1=c(0,1),p2=c(0,100))
-set.seed(456)
-
control <- list(
nseq=4,
thin.t=10
)
-pars <- list(p1=c(0,1),p2=c(0,100))
+## TODO: check if FUN Missing arguments: is really necessary
+## TODO: header for iteration output
+## TODO: predict.dream-model
+## TODO: dreamCalibrate help
+set.seed(456)
+dd <- dreamCalibrate(FUN=Model.y,
+ pars=pars,
+ obs=obs.all$y,
+ FUN.pars=list(x=obs.all$x),
+ control=control
+ )
-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)
-dd
+print(summary(dd))
-summary(dd)
+print(coef(dd))
-coef(dd)
-
plot(dd)
plotFME()
@@ -98,16 +96,15 @@
}##plotCIs
## Calibrate with Obs
-dd <- dream(
- FUN=Model.y, func.type="calc.rmse",
- pars = pars,
- FUN.pars=list(
- x=Obs$x
- ),
- INIT = LHSInit,
- measurement=list(data=Obs$y),
- control = control
- )
+dd <- dreamCalibrate(
+ FUN=Model.y,
+ pars=pars,
+ obs=Obs$y,
+ FUN.pars=list(
+ x=Obs$x
+ ),
+ control = control
+ )
##Obs1
plotFME()
@@ -125,26 +122,9 @@
### Example with new sample
dd.sim <- simulate(dd)
-predict(dd.sim)
plotFME()
lines(Obs2$x,predict(dd,list(x=Obs2$x)),col="blue")
lines(Obs2$x,predict(dd.sim,list(x=Obs2$x)),col="purple")
-########################
-## Legacy examples
-plotFME()
-lines(Model(p=coef(dd),x=0:375),col="green")
-
-## Naive 95% bounds from residuals
-resid <- Model.y(p=coef(dd),x=Obs$x)-Obs$y
-##densityplot(resid)
-qq <- quantile(resid,c(0.005,.995))
-gg <- t(sapply(Model.y(p=coef(dd),x=Obs$x),function(v) v+qq))
-plotCIs(Obs$x,gg,col="grey")
-
-## Test on Obs2
-cis.2 <- predict(dd,list(x=Obs2$x),out="CI")
-plotCIs(Obs2$x,cis.2,col="red")
-
## TODO: add residual error, using method p6, vrugt. equifinality
Modified: pkg/demo/FME.nonlinear.model_parallelisation.R
===================================================================
--- pkg/demo/FME.nonlinear.model_parallelisation.R 2010-05-12 05:58:05 UTC (rev 28)
+++ pkg/demo/FME.nonlinear.model_parallelisation.R 2010-05-13 03:34:01 UTC (rev 29)
@@ -37,23 +37,24 @@
for (p in c("none","snow","foreach")){
- set.seed(456)
+ print(p)
control$parallel <- p
-
- 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
- )
+
+ set.seed(456)
+ dd <- dreamCalibrate(
+ FUN=Model.y,
+ pars = pars,
+ obs=obs.all$y,
+ FUN.pars=list(
+ x=obs.all$x
+ ),
+ control = control
+ )
- print(dd$control$parallel)
+
+ print("Coefficients:")
print(coef(dd))
- print(dd$time)
+ print(sprintf("Elapsed time %f seconds",dd$time))
}
if (require(doSNOW)) stopCluster(cl)
@@ -83,34 +84,31 @@
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")){
+ print(p)
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)
+ dd <- dreamCalibrate(
+ FUN=Model.y,
+ pars = pars,
+ obs=obs.all$y,
+ FUN.pars=list(
+ x=obs.all$x
+ ),
+ control = control
+ )
+ print(sprintf("Number of function evaluations: %f",dd$fun.evals))
}
if (require(doSNOW)) stopCluster(cl)
Modified: pkg/demo/run_example1.R
===================================================================
--- pkg/demo/run_example1.R 2010-05-12 05:58:05 UTC (rev 28)
+++ pkg/demo/run_example1.R 2010-05-13 03:34:01 UTC (rev 29)
@@ -17,3 +17,5 @@
),
control = control
)
+
+plot(dd)
Modified: pkg/man/coef.dream.Rd
===================================================================
--- pkg/man/coef.dream.Rd 2010-05-12 05:58:05 UTC (rev 28)
+++ pkg/man/coef.dream.Rd 2010-05-13 03:34:01 UTC (rev 29)
@@ -1,18 +1,15 @@
\name{coef.dream}
\alias{coef.dream}
-\title{Extract parameter values from dream object}
+\alias{coef.dream_model}
+\title{Extract parameter values from dream or dream_model object}
\usage{
-\method{coef}{dream}(object, last.prop = 0.5, use.thinned = FALSE,
- method = c("uni.mode", "mean", "median", "sample.ml"), \dots)
+\method{coef}{dream}(object, method = c("uni.mode", "mean", "median", "sample.ml"), \dots)
}
\description{Extract parameter values using a choice of methods (or an
arbitrary function)}
\value{named vector of parameter values}
\arguments{
\item{object}{dream object}
- \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}{method for extracting a parameter set from the MCMC
chains. One of:
\describe{
@@ -36,9 +33,10 @@
}
}
}
- \item{...}{Unused. Matches generic}
+ \item{...}{Passed to \code{\link{window.dream}}}
}
\details{
+
e.g. of using arbitrary function for method: 20\% quantile.
\code{
coef(object,method=function(sss) apply(as.matrix(sss),2,quantile,0.2)
Modified: pkg/man/dream.Rd
===================================================================
--- pkg/man/dream.Rd 2010-05-12 05:58:05 UTC (rev 28)
+++ pkg/man/dream.Rd 2010-05-13 03:34:01 UTC (rev 29)
@@ -67,6 +67,7 @@
}
}
\details{
+
Elements of control are:
\tabular{lll}{
\emph{Element} \tab \emph{Default} \tab \emph{Description} \cr
@@ -138,20 +139,15 @@
There are S3 methods for print, summary,
\code{\link[=plot.dream]{plot}},
\code{\link[=coef.dream]{coef}},
- \code{\link[=predict.dream]{predict}},
- simulate
+ \code{\link[=simulate.dream]{simulate}},
+ \code{\link[=window.dream]{window}}
- 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{
A list with elements:
\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{Sequences}{ \code{\link{mcmc.list}}. \code{nseq} mcmc elements
+ of \code{ndim} variables. The thinned chains 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}
@@ -166,11 +162,10 @@
}
\seealso{
- \code{\link{possibility.envelope}}
-
+ \code{\link{dreamCalibrate}} to calibrate a function using dream.
+
Examples in demo folder:
\itemize{
- \item{FME non linear model: }{Calibrating the non-linear model shown in the FME vignette}
\item{example1: }{Fitting a banana shaped distribution - the first example
in Vrugt matlab code}
}
Added: pkg/man/dreamCalibrate.Rd
===================================================================
--- pkg/man/dreamCalibrate.Rd (rev 0)
+++ pkg/man/dreamCalibrate.Rd 2010-05-13 03:34:01 UTC (rev 29)
@@ -0,0 +1,60 @@
+\name{dreamCalibrate}
+\alias{dreamCalibrate}
+\title{Utility to calibrate a function using dream}
+\usage{
+
+dreamCalibrate(FUN, pars, obs, lik.fun=calc.loglik,lik.control=NULL,FUN.pars = list(), ...)
+
+}
+\description{
+ Calibrate a function using \code{\link{dream}}, a specified
+ likelihood function \code{lik.fun} and observed values \code{obs}
+}
+\arguments{
+ \item{FUN}{
+ model function with first argument a vector of parameter values of
+ length ndim.
+ }
+ \item{pars}{
+ a list of variable ranges. Any names will be propagated to output.
+ }
+ \item{obs}{
+ a numeric vector of observed values, corresponding to the output of \code{FUN}
+ }
+ \item{lik.fun}{
+ A function that returns the log likelihood of model predictions
+ matching observed values. \code{log.lik=f(predicted,observed,control)}
+ }
+ \item{lik.control}{
+ A list of any extra arguments to be passed to \code{lik.fun}
+ }
+ \item{FUN.pars}{
+ A list of any extra arguments to be passed to \code{FUN}.
+ }
+ \item{...}{
+ Extra arguments to be passed to dream, e.g. \code{control}
+ }
+}
+\details{
+
+ There are S3 methods for:
+ \code{\link[=predict.dream_model]{predict}},
+ \code{\link[=coef.dream]{coef}}.
+
+}
+\value{
+ An object inheriting from \code{\link{dream}}, i.e. with the same elements and:
+ \item{FUN}{The function calibrated}
+ \item{FUN.pars}{The extra arguments originally passed to that function}
+}
+
+\seealso{
+ See \code{\link{dream}} for details on the calibration method,
+ visualisation of its results and diagnostics.
+
+ Example in demo folder:
+ \itemize{
+ \item{FME non linear model: }{Calibrating the non-linear model shown
+ in the FME vignette}
+ }
+}
Property changes on: pkg/man/dreamCalibrate.Rd
___________________________________________________________________
Added: svn:executable
+ *
Deleted: pkg/man/maxLikPars.Rd
===================================================================
--- pkg/man/maxLikPars.Rd 2010-05-12 05:58:05 UTC (rev 28)
+++ pkg/man/maxLikPars.Rd 2010-05-13 03:34:01 UTC (rev 29)
@@ -1,10 +0,0 @@
-\name{maxLikPars}
-\alias{maxLikPars}
-\title{Naive maximum density parameter selection...}
-\usage{maxLikPars(ss, ...)}
-\description{Naive maximum density parameter selection}
-\value{named character vector of parameter values
-
-Uses which.max and density function}
-\arguments{\item{ss}{MCMC chains. mcmc or mcmclist object}
-\item{...}{extra parameters to pass to density}}
Deleted: pkg/man/possibility.envelope.Rd
===================================================================
--- pkg/man/possibility.envelope.Rd 2010-05-12 05:58:05 UTC (rev 28)
+++ pkg/man/possibility.envelope.Rd 2010-05-13 03:34:01 UTC (rev 29)
@@ -1,14 +0,0 @@
-\name{possibility.envelope}
-\alias{possibility.envelope}
-\title{Obtain prediction confidence intervals for model, around new input values...}
-\usage{possibility.envelope(dd, FUN.pars, ndraw=1000, conf=99)}
-\description{Obtain prediction confidence intervals for model, around new input values}
-\arguments{\item{dd}{dream object}
-\item{FUN.pars}{new par values. if missing, use those from dream object}
-\item{ndraw}{number of new iterations}
-\item{conf}{Percentage two-sided confidence interval}}
-\details{
-Progresses chains that are assumed to have converged and thins them, to obtain a sample
-of size \code{ndraw} from the posterior parameter distribution. Evaluates the function for
-all parameter sets and reports \code{conf}\% confidence interval.
-}
\ No newline at end of file
Modified: pkg/man/predict.dream.Rd
===================================================================
--- pkg/man/predict.dream.Rd 2010-05-12 05:58:05 UTC (rev 28)
+++ pkg/man/predict.dream.Rd 2010-05-13 03:34:01 UTC (rev 29)
@@ -1,24 +1,23 @@
-\name{predict.dream}
-\alias{predict.dream}
-\title{Predict values using dream object}
+\name{predict.dream_model}
+\alias{predict.dream_model}
+\title{Predict values using dream_model object}
\usage{
-\S3method{predict}{dream}(object, newdata=NULL,newFUN=NULL,method="maxLik", level=0.99, last.prop=0.5, use.thinned=TRUE,\dots)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/dream -r 29
More information about the Dream-commits
mailing list