[Dream-commits] r16 - in pkg: . R demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Mar 3 05:55:19 CET 2010
Author: josephguillaume
Date: 2010-03-03 05:55:18 +0100 (Wed, 03 Mar 2010)
New Revision: 16
Added:
pkg/R/predict.dream.R
pkg/demo/
pkg/demo/00Index
pkg/demo/FME.nonlinear.model.R
pkg/demo/example1.R
pkg/man/predict.dream.Rd
Removed:
pkg/demo/FME nonlinear model.R
pkg/demo/example1.R
pkg/demos/
Modified:
pkg/NAMESPACE
pkg/R/CompDensity.R
pkg/R/coef.dream.R
pkg/R/dream.R
pkg/R/possibility.envelope.R
pkg/man/possibility.envelope.Rd
Log:
predict, simulate
outlier removal in burn in period
return to burn in period if outliers are detected later
coef allows choosing mean and median
added help files for coef.dream and predict.dream
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2010-03-01 23:16:27 UTC (rev 15)
+++ pkg/NAMESPACE 2010-03-03 04:55:18 UTC (rev 16)
@@ -9,3 +9,5 @@
S3method(coef,dream)
S3method(plot,dream)
S3method(print,dream)
+S3method(predict,dream)
+S3method(simulate,dream)
Modified: pkg/R/CompDensity.R
===================================================================
--- pkg/R/CompDensity.R 2010-03-01 23:16:27 UTC (rev 15)
+++ pkg/R/CompDensity.R 2010-03-03 04:55:18 UTC (rev 16)
@@ -33,12 +33,11 @@
## modpred. scalar or vector commensurate to measurement$data
## err. vector of same length as modpred
## SSR scalar
+ ## temp. list of length nseq with elements of length 2: p and logp
- p <- rep(NA,control$nseq)
- logp <- rep(NA,control$nseq)
+ ## Ready for multicore
+ temp <- lapply(1:nrow(pars),function (ii){
- ## Sequential evaluation
- for (ii in 1:nrow(pars)){
## Call model to generate simulated data
## TODO: correct use of optional pars?
modpred <- FUN(pars[ii,],...)
@@ -46,18 +45,18 @@
switch(func.type,
## Model directly computes posterior density
posterior.density={
- p[ii] <- modpred
- logp[ii] <- log(modpred)
+ p <- modpred
+ logp <- log(modpred)
},
## Model computes output simulation
calc.loglik={
err <- as.numeric(measurement$data-modpred)
## Derive the log likelihood
- logp[ii] <- measurement$N*log(control$Wb/measurement$sigma)-
+ logp <- measurement$N*log(control$Wb/measurement$sigma)-
control$Cb*(sum((abs(err/measurement$sigma))^(2/(1+control$gamma))))
## And retain in memory
- p[ii] <- logp[ii]
+ p <- logp
},
## Model computes output simulation
## TODO: may need as.numeric
@@ -67,14 +66,14 @@
## Derive the sum of squared error
SSR <- sum(abs(err)^(2/(1+control$gamma)))
## And retain in memory
- p[ii] <- -SSR
- logp[ii] <- -0.5*SSR
+ p <- -SSR
+ logp <- -0.5*SSR
},
## Model directly computes log posterior density
logposterior.density={
- p[ii] <- modpred
- logp[ii] <- modpred
+ p <- modpred
+ logp <- modpred
},
## Similar as 3, but now weights with the Measurement Sigma
## TODO: identical to rmse because difference is in metrop
@@ -84,11 +83,15 @@
## Derive the sum of squared error
SSR <- sum(abs(err)^(2/(1+control$gamma)))
## And retain in memory
- p[ii] <- -SSR
- logp[ii] <- -0.5*SSR
+ p <- -SSR
+ logp <- -0.5*SSR
}) ##switch
- } ## for rows
+ c(p,logp)
+ }) ## lapply
+ p <- sapply(temp,function(x) x[1])
+ logp <- sapply(temp,function(x) x[2])
+
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-01 23:16:27 UTC (rev 15)
+++ pkg/R/coef.dream.R 2010-03-03 04:55:18 UTC (rev 16)
@@ -1,16 +1,27 @@
##' Extract maximum likelihood parameter values
##' @param dream object
-##' @last.prop proportion of total sequence to use (0,1]
+##' @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
##' @return named vector of parameter values
-coef.dream <- function(object,last.prop=.5,use.thinned=FALSE,...){
+coef.dream <- function(object,last.prop=.5,use.thinned=FALSE,
+ method=c("maxLik","mean","median"),...){
stopifnot(last.prop>0)
if (use.thinned) sss <- object$Reduced.Seq
else sss <- object$Sequences
+
+ if (class(method)!="function") {
+ method <- switch(
+ match.arg(method),
+ "maxLik"=maxLikPars,
+ "mean"=function(sss) colMeans(as.matrix(sss)),
+ "median"=function(sss) apply(as.matrix(sss),2,median)
+ )
+ }
- if (last.prop==1) return(maxLikPars(sss))
+ if (last.prop==1) return(method(sss))
else {
ss <- window(object$Sequences, start = end(sss)*(1-last.prop) + 1)
- return(maxLikPars(ss))
+ return(method(ss))
}
}##coef.dream
Modified: pkg/R/dream.R
===================================================================
--- pkg/R/dream.R 2010-03-01 23:16:27 UTC (rev 15)
+++ pkg/R/dream.R 2010-03-03 04:55:18 UTC (rev 16)
@@ -9,6 +9,7 @@
outlierTest = 'IQR_test', ## What kind of test to detect outlier chains?
pCR.Update = TRUE, ## Adaptive tuning of crossover values
boundHandling = 'reflect', ## Boundary handling: "reflect", "bound", "fold", "none"
+ burnin.length=0.1, ## Proportion of iterations considered to be burnin. 0 to turn off.
### Termination criteria. TODO: are the 2nd two valid, given that ndraw is used in adaptive pcr
ndraw = 1e5, ## maximum number of function evaluations
maxtime = Inf, ## maximum duration of optimization in seconds
@@ -99,9 +100,16 @@
## Check INIT and FUN have required extra parameters in INIT.pars & FUN.pars
req.args.init <- names(formals(INIT))
req.args.FUN <- names(formals(FUN))
+ req.args.FUN <- req.args.FUN[2:length(req.args.FUN)] ##optional pars only
+
+ 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(!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=" "))
- if(!all(req.args.FUN[2:length(req.args.FUN)] %in% c(names(FUN.pars)))) stop(paste(c("FUN Missing extra arguments:",req.args.FUN[!req.args.FUN %in% c("x",names(FUN.pars))]),sep=" "))
+ 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
@@ -111,7 +119,7 @@
stop("unrecognised options: ",
toString(names(control)[!isValid]))
- if (!is.null("measurement")){
+ if (!is.null(measurement)){
if (! "sigma" %in% names(measurement)) measurement$sigma <- sd(measurement$data)
if (! "N" %in% names(measurement)) measurement$N <- length(measurement$data)
}
@@ -124,7 +132,8 @@
## Correct to match nseq
control$REPORT <- (control$REPORT%/%control$nseq) * control$nseq
-
+ if (control$burnin.length<1) control$burnin.length <- control$burnin.length*control$ndraw
+
## Check validity of settings
if (control$DEpairs==0) stop("control$DEpairs set to 0. Increase nseq?")
stopifnot(control$DEpairs<=(control$nseq-1)/2) ## Requirement of offde
@@ -149,6 +158,8 @@
## Max number of times through loops
max.counter <- ceiling((control$ndraw+control$steps*NSEQ)/NSEQ)+1
max.counter.outloop <- ceiling((control$ndraw+control$steps*NSEQ)/NSEQ/control$steps)+1
+
+ if (!is.na(control$thin.t) && floor(max.counter/control$thin.t)<2) stop(sprintf("Thin parameter should be much smaller than number of iterations (currently %d,%d)",control$thin.t,max.counter))
## Calculate the parameters in the exponential power density function of Box and Tiao (1973)
cbwb <- CalcCbWb(control$gamma)
@@ -176,7 +187,9 @@
lCR <- NSEQ*control$steps
} ##pCR.Update
+ end.burnin <- control$burnin.length
+
############################
## Initialise output object
@@ -185,8 +198,10 @@
obj$call <- match.call()
obj$control <- control
- EXITFLAG <- NA
- EXITMSG <- NULL
+ obj$in.burnin <- TRUE
+
+ obj$EXITFLAG <- NA
+ obj$EXITMSG <- NULL
## counter.fun.evals + AR at each step
obj$AR<-matrix(NA,max.counter,2)
@@ -214,7 +229,7 @@
## Check whether will save a reduced sample
if (!is.na(control$thin.t)){
iloc.2 <- 0
- Reduced.Seq <- array(NA,c(floor(max.counter/control$thin.t),NDIM+2,NSEQ))
+ Reduced.Seq <- array(NA,c(ceiling(max.counter/control$thin.t),NDIM+2,NSEQ))
} else Reduced.Seq <- NULL
############################
@@ -353,32 +368,42 @@
## Do this to get rounded iteration numbers
if (counter.outloop == 2) control$steps <- control$steps + 1
+ outliers <- RemOutlierChains(X,hist.logp[1:(counter-1),],control)
+
## Check whether to update individual pCR values
- if (counter.fun.evals <= 0.1 * control$ndraw) {
+ if (counter.fun.evals <= end.burnin) {
if (control$pCR.Update) {
## Update pCR values
tmp <- AdaptpCR(CR, delta.tot, lCR, control)
pCR <- tmp$pCR
lCR <- tmp$lCR
}
- } else {
- ## See whether there are any outlier chains, and remove them to current best value of X
- tmp <- RemOutlierChains(X,hist.logp[1:(counter-1),],control)
+
+ ## Change any outlier chains to current best value of X
+ ## TODO: matlab code doesn't match paper. Outlier removal should be within burnin period
+
## Loop over each outlier chain (if length>0)
- for (out.id in tmp$chain.id){
+ for (out.id in outliers$chain.id){
## Draw random other chain -- cannot be the same as current chain
- r.idx <- which.max(tmp$mean.hist.logp)
+ r.idx <- which.max(outliers$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[1:(counter-1),r.idx]
- ## Jump outlier chain to r_idx -- Sequences
+ ## Jump outlier chain to r_idx -- Sequences, X
Sequences[iloc,1:(NDIM+2),out.id] <- X[r.idx,]
- ## Jump outlier chain to r_idx -- X
X[out.id,1:(NDIM+2)] <- X[r.idx,]
- ## Add to chainoutlier
+ ## Add to outlier tracker
obj$outlier <- rbind(obj$outlier,c(counter.fun.evals,out.id))
} ##for remove outliers
- } ##else
+ } else if (control$burnin.length!=0){
+ obj$in.burnin <- FALSE
+ if (length(outliers)>0){
+ warning("Outliers detected outside burn-in period, returning to burn in")
+ obj$in.burnin <- TRUE
+ end.burnin <- counter.fun.evals+control$burnin.length
+ }
+ } ##in burn in.
+
if (control$pCR.Update) {
## Generate CR values based on current pCR values
Modified: pkg/R/possibility.envelope.R
===================================================================
--- pkg/R/possibility.envelope.R 2010-03-01 23:16:27 UTC (rev 15)
+++ pkg/R/possibility.envelope.R 2010-03-03 04:55:18 UTC (rev 16)
@@ -12,8 +12,12 @@
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)
Added: pkg/R/predict.dream.R
===================================================================
--- pkg/R/predict.dream.R (rev 0)
+++ pkg/R/predict.dream.R 2010-03-03 04:55:18 UTC (rev 16)
@@ -0,0 +1,78 @@
+##' 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,
+##' using various methods of summarising the posterior parameter and output distributions
+##' @param object dream object
+##' @param newdata. new FUN.pars list. If NULL, use object's.
+##' @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
+##' @return data frame with each column corresponding to a returned vector
+predict.dream <- function(object,newdata=NULL,
+ method="maxLik",level=0.99,
+ last.prop=0.5,use.thinned=TRUE
+ ){
+
+ ## Check and initialise parameters
+ stopifnot(is.null(newdata) || is.list(newdata))
+ stopifnot(!"CI" %in% method || !is.null(level))
+ stopifnot(last.prop>0)
+
+###
+ ## Fetch function and parameters from dream object
+
+ if (is.null(newdata)) newdata <- eval(object$call$FUN.pars)
+ par.name <- names(formals(eval(object$call$FUN)))[1]
+
+ wrap <- function(p) {
+ newdata[[par.name]] <- p
+ as.numeric(do.call(eval(object$call$FUN),newdata))
+ }
+###
+
+ ## Predict for desired method(s)
+
+ 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)
+
+ ff <- apply(as.matrix(sss),1,wrap)
+
+ return(t(apply(ff,1,quantile,c((1-level)/2,1-(1-level)/2))))
+
+ } else {
+ return(wrap(coef(object,method=method,
+ use.thinned=use.thinned,last.prop=last.prop)
+ ))
+ }
+
+} ##predict.dream
Property changes on: pkg/R/predict.dream.R
___________________________________________________________________
Name: svn:executable
+ *
Copied: pkg/demo (from rev 12, pkg/demos)
Added: pkg/demo/00Index
===================================================================
--- pkg/demo/00Index (rev 0)
+++ pkg/demo/00Index 2010-03-03 04:55:18 UTC (rev 16)
@@ -0,0 +1,2 @@
+FME.nonlinear.model Nonlinear model calibration as with FME vignette.
+run_example1 n-dimensional banana shaped gaussian distribution. From Vrugt's matlab code.
Property changes on: pkg/demo/00Index
___________________________________________________________________
Name: svn:executable
+ *
Deleted: pkg/demo/FME nonlinear model.R
===================================================================
--- pkg/demos/FME nonlinear model.R 2010-02-26 04:31:57 UTC (rev 12)
+++ pkg/demo/FME nonlinear model.R 2010-03-03 04:55:18 UTC (rev 16)
@@ -1,140 +0,0 @@
-
-## 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
- 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)
-
-
-###########################
-### FME results
-
-Model <- function(p,x) return(data.frame(x=x,y=p[1]*x/(x+p[2])))
-##Model(c(0.1,1),obs.all$x)
-
-
-library(FME)
-Residuals <- function(p) {
- cost<-modCost(model=Model(p,Obs$x),obs=Obs,x="x")
- modCost(model=Model(p,Obs2$x),obs=Obs2,cost=cost,x="x")
-}
-
-P <- modFit(f=Residuals,p=c(0.1,1))
-print(P$par)
-
-plotFME <- function(){
-plot(Obs,xlab="mg COD/l",ylab="1/hour", pch=16, cex=1.5,
- xlim=c(25,400),ylim=c(0,0.15))
-points(Obs2,pch=18,cex=1.5, col="red")
-lines(Model(p=P$par,x=0:375))
-}
-
-
-##########################
-## DREAM results
-
-
-unloadNamespace("dream")
-library(dream)
-
-
-Model.y <- function(p,x) p[1]*x/(x+p[2])
-
-control <- list(
- nseq=5,
- gamma=0,
- nCR=3,
- ndraw=1e5,
- steps=10,
- eps=5e-2,
- outlierTest='IQR_test',
- pCR.Update=TRUE,
- thin=TRUE,
- thin.t=10,
- boundHandling="fold",
- Rthres=1+1e-3
- )
-
-pars <- list(p1=c(0,1),p2=c(0,100))
-
-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
- )
-
-dd
-
-summary(dd)
-
-coef(dd)
-
-plot(dd)
-
-plotFME()
-lines(Model(p=coef(dd),x=0:375),col="green")
-
-
-########################
-## Calculate bounds around output estimates
-
-plotCIs <- function(x,cis,...){
- em <- strwidth("M")/2
- segments(x0=x,y0=cis[,1],
- x1=x,y1=cis[,2],
- ...
- )
- segments(x0=x-em,y0=cis[,1],
- x1=x+em,y1=cis[,1],
- ...
- )
- segments(x0=x-em,y0=cis[,2],
- x1=x+em,y1=cis[,2],
- ...
- )
-}##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.all$y),
- control = control
- )
-
-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")
-
-## Bounds on Obs
-cis.1 <- possibility.envelope(dd)
-plotCIs(Obs$x,cis.1,col="black")
-
-## Test on Obs2
-cis.2 <- possibility.envelope(dd,list(x=Obs2$x))
-plotCIs(Obs2$x,cis.2,col="red")
-
-
-## TODO: add residual error, using method p6, vrugt. equifinality
Copied: pkg/demo/FME.nonlinear.model.R (from rev 15, pkg/demos/FME.nonlinear.model.R)
===================================================================
--- pkg/demo/FME.nonlinear.model.R (rev 0)
+++ pkg/demo/FME.nonlinear.model.R 2010-03-03 04:55:18 UTC (rev 16)
@@ -0,0 +1,154 @@
+
+## 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
+ 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)
+
+
+###########################
+### FME results
+
+Model <- function(p,x) return(data.frame(x=x,y=p[1]*x/(x+p[2])))
+##Model(c(0.1,1),obs.all$x)
+
+
+library(FME)
+Residuals <- function(p) {
+ cost<-modCost(model=Model(p,Obs$x),obs=Obs,x="x")
+ modCost(model=Model(p,Obs2$x),obs=Obs2,cost=cost,x="x")
+}
+
+P <- modFit(f=Residuals,p=c(0.1,1))
+print(P$par)
+
+plotFME <- function(){
+plot(Obs,xlab="mg COD/l",ylab="1/hour", pch=16, cex=1.5,
+ xlim=c(25,400),ylim=c(0,0.15))
+points(Obs2,pch=18,cex=1.5, col="red")
+lines(Model(p=P$par,x=0:375))
+}
+
+
+##########################
+## DREAM results
+
+
+unloadNamespace("dream")
+library(dream)
+
+
+Model.y <- function(p,x) as.ts(p[1]*x/(x+p[2]))
+
+set.seed(456)
+
+control <- list(
+ nseq=4
+## REPORT=0
+## ndraw=1000
+ ## Rthres=1+1e-3
+ )
+
+pars <- list(p1=c(0,1),p2=c(0,100))
+
+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
+ )
+
+dd
+
+summary(dd)
+
+coef(dd)
+
+plot(dd)
+
+plotFME()
+lines(Model(p=coef(dd),x=0:375),col="green")
+
+
+########################
+## Calculate bounds around output estimates
+
+plotCIs <- function(x,cis,...){
+ em <- strwidth("M")/2
+ segments(x0=x,y0=cis[,1],
+ x1=x,y1=cis[,2],
+ ...
+ )
+ segments(x0=x-em,y0=cis[,1],
+ x1=x+em,y1=cis[,1],
+ ...
+ )
+ segments(x0=x-em,y0=cis[,2],
+ x1=x+em,y1=cis[,2],
+ ...
+ )
+}##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
+ )
+
+##Obs1
+plotFME()
+lines(Obs$x,predict(dd),col="blue")
+lines(Obs$x,predict(dd,method="mean"),col="red")
+lines(Obs$x,predict(dd,method="median"),col="orange")
+plotCIs(Obs$x,predict(dd,method="CI"),col="black")
+
+
+##Obs2
+plotFME()
+lines(Obs2$x,predict(dd,list(x=Obs2$x)),col="blue")
+lines(Obs2$x,predict(dd,list(x=Obs2$x),method="mean"),col="red")
+lines(Obs2$x,predict(dd,list(x=Obs2$x),method="median"),col="orange")
+plotCIs(Obs2$x,predict(dd,list(x=Obs2$x),method="CI"),col="red")
+
+### 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
Deleted: pkg/demo/example1.R
===================================================================
--- pkg/demos/example1.R 2010-02-26 04:31:57 UTC (rev 12)
+++ pkg/demo/example1.R 2010-03-03 04:55:18 UTC (rev 16)
@@ -1,47 +0,0 @@
-unloadNamespace("dream")
-library(dream)
-
-## n-dimensional banana shaped gaussian distribution
-## Hyperbolic shaped posterior probability distribution
-## @param x vector of length ndim
-## @param bpar banana-ness. scalar.
-## @param imat matrix ndim x ndim
-## @return SSR. scalar
-## Cursorily verified to match matlab version
-Banshp <- function(x,bpar,imat){
- x[2] <- x[2]+ bpar * x[1]^2 - 100*bpar
- S <- -0.5 * (x %*% imat) %*% as.matrix(x)
- return(S)
-}
-
-control <- list(
- ndim=10,
- DEpairs=3,
- gamma=0,
- nCR=3,
- ndraw=1e5,
- steps=10,
- eps=5e-2,
- outlierTest='IQR_test',
- pCR.Update=TRUE,
- thin=TRUE,
- thin.t=10,
- boundHandling='none'
- )
-
-
-pars <- lapply(1:control$ndim,function(x) c(-Inf,Inf))
-names(pars) <- paste("b",1:control$ndim,sep="")
-cmat <- diag(control$ndim)
-cmat[1,1] <- 100
-muX=rep(0,control$ndim)
-qcov=diag(control$ndim)*5
-
-FUN=Banshp
-func.type="logposterior.density"
-FUN.pars=list(
- imat=solve(cmat),
- bpar=0.1
- )
-
-
Copied: pkg/demo/example1.R (from rev 14, pkg/demos/example1.R)
===================================================================
--- pkg/demo/example1.R (rev 0)
+++ pkg/demo/example1.R 2010-03-03 04:55:18 UTC (rev 16)
@@ -0,0 +1,46 @@
+unloadNamespace("dream")
+library(dream)
+
+## n-dimensional banana shaped gaussian distribution
+## Hyperbolic shaped posterior probability distribution
+## @param x vector of length ndim
+## @param bpar banana-ness. scalar.
+## @param imat matrix ndim x ndim
+## @return SSR. scalar
+## Cursorily verified to match matlab version
+Banshp <- function(x,bpar,imat){
+ x[2] <- x[2]+ bpar * x[1]^2 - 100*bpar
+ S <- -0.5 * (x %*% imat) %*% as.matrix(x)
+ return(S)
+}
+
+control <- list(
+ ndim=10,
+ DEpairs=3,
+ gamma=0,
+ nCR=3,
+ ndraw=1e5,
+ steps=10,
+ eps=5e-2,
+ outlierTest='IQR_test',
+ pCR.Update=TRUE,
+ thin.t=10,
+ boundHandling='none'
+ )
+
+
+pars <- lapply(1:control$ndim,function(x) c(-Inf,Inf))
+names(pars) <- paste("b",1:control$ndim,sep="")
+cmat <- diag(control$ndim)
+cmat[1,1] <- 100
+muX=rep(0,control$ndim)
+qcov=diag(control$ndim)*5
+
+FUN=Banshp
+func.type="logposterior.density"
+FUN.pars=list(
+ imat=solve(cmat),
+ bpar=0.1
+ )
+
+
Modified: pkg/man/possibility.envelope.Rd
===================================================================
--- pkg/man/possibility.envelope.Rd 2010-03-01 23:16:27 UTC (rev 15)
+++ pkg/man/possibility.envelope.Rd 2010-03-03 04:55:18 UTC (rev 16)
@@ -7,3 +7,8 @@
\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
Added: pkg/man/predict.dream.Rd
===================================================================
--- pkg/man/predict.dream.Rd (rev 0)
+++ pkg/man/predict.dream.Rd 2010-03-03 04:55:18 UTC (rev 16)
@@ -0,0 +1,14 @@
+\name{predict.dream}
+\alias{predict.dream}
+\title{Predict values using dream object}
+\usage{predict.dream(object, newdata, method="maxLik", level=0.99, last.prop=0.5, use.thinned=TRUE)}
+\description{
+Predict values using function calibrated by dream, optionally with new data,
+using various methods of summarising the posterior parameter and output distributions}
+\value{data frame with each column corresponding to a returned vector}
+\arguments{\item{object}{dream object}
+\item{newdata.}{new FUN.pars list. If NULL, use object's.}
+\item{method}{CI or a \code{method} of \code{\link\{coef.dream}}
+\item{level.}{Requested two-sided level of confidence. For CI method.}
+\item{last.prop}{Proportion of MCMC chains to keep}
+\item{use.thinned}{Whether to use existing thinned chains}}
Property changes on: pkg/man/predict.dream.Rd
___________________________________________________________________
Name: svn:executable
+ *
More information about the Dream-commits
mailing list