[Dream-commits] r31 - in pkg: R demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon May 17 04:38:57 CEST 2010
Author: josephguillaume
Date: 2010-05-17 04:38:57 +0200 (Mon, 17 May 2010)
New Revision: 31
Modified:
pkg/R/coef.dream.R
pkg/R/dream.R
pkg/R/dreamCalibrate.R
pkg/R/predict.dream.R
pkg/demo/FME.nonlinear.model.R
pkg/man/coef.dream.Rd
pkg/man/dreamCalibrate.Rd
pkg/man/predict.dream.Rd
pkg/man/window.dream.Rd
Log:
- Made sample.ml default coef function - will always give a parameter set with a good likelihood value, regardless of multimodal, density parameters
- Changed thinning to use window.mcmc. Reduced.Seq code is no longer used.
- Synced sample.ml to thinned output of dream by matching mcmc characteristics for window.mcmc
- Added basic text to window.mcmc help
Modified: pkg/R/coef.dream.R
===================================================================
--- pkg/R/coef.dream.R 2010-05-13 07:55:08 UTC (rev 30)
+++ pkg/R/coef.dream.R 2010-05-17 02:38:57 UTC (rev 31)
@@ -3,20 +3,20 @@
##' @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,method=c("uni.mode","mean","median","sample.ml"),...)
+coef.dream <- function(object,method=c("sample.ml","uni.mode","mean","median"),...)
{
-
+
sss <- window(object,...)
+
+ if (!is.function(method) && identical(match.arg(method), "sample.ml")) {
+ ## TODO: make sure ppp corresponds to sss
+ ppp <- window(as.mcmc(object$hist.logp),start=start(sss),thin=thin(sss))
+ maxi <- which.max(ppp)
+ maxchain <- col(ppp)[maxi]
+ maxtime <- row(ppp)[maxi]
+ return(sss[[maxchain]][maxtime,])
+ }
- 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),
@@ -27,4 +27,4 @@
}
return(method(sss))
-}##coef.dream
+} ##coef.dream
Modified: pkg/R/dream.R
===================================================================
--- pkg/R/dream.R 2010-05-13 07:55:08 UTC (rev 30)
+++ pkg/R/dream.R 2010-05-17 02:38:57 UTC (rev 31)
@@ -524,18 +524,16 @@
## Trim outputs to collected data - remove extra rows
## Convert sequences to mcmc objects
- if (!is.na(control$thin.t)){
- Reduced.Seq <- Reduced.Seq[1:counter.redseq,,]
- 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: remove all prior refs to Reduced.Seq - now not needed, given window.mcmc is used
+ 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)){
+ obj$Sequences <- window(obj$Sequences,thin=control$thin.t)
}
## TODO: make these 'ts' objects and sync with Reduced.Seq by thinning
+ ## Would it be better to keep all data, and sync when needed by matching start, end,thin?
+ ## See coef.dream for eg.
obj$X <- X
obj$R.stat <- obj$R.stat[1:counter.report,,drop=FALSE]
obj$hist.logp <- hist.logp[1:(counter-1),,drop=FALSE]
Modified: pkg/R/dreamCalibrate.R
===================================================================
--- pkg/R/dreamCalibrate.R 2010-05-13 07:55:08 UTC (rev 30)
+++ pkg/R/dreamCalibrate.R 2010-05-17 02:38:57 UTC (rev 31)
@@ -56,6 +56,7 @@
... ##Extra arguments to dream
){
if (is.null(lik.control)) lik.control <- eval(formals(lik.fun)$control)
+ ## TODO: remove ... from here and move do.call processing from dream to here.
wrap.lik.fun <- function(pars,...) lik.fun(FUN(pars,...),obs,lik.control)
dd <- dream(FUN=wrap.lik.fun,
pars=pars,
@@ -67,6 +68,7 @@
dd$call <- match.call()
dd$FUN <- FUN
dd$FUN.pars <- FUN.pars
+ dd$lik.fun <- function(p) do.call(wrap.lik.fun,modifyList(FUN.pars,list(pars=p)))
class(dd) <- c("dream_model",class(dd))
dd
} ##dreamCalibrate
Modified: pkg/R/predict.dream.R
===================================================================
--- pkg/R/predict.dream.R 2010-05-13 07:55:08 UTC (rev 30)
+++ pkg/R/predict.dream.R 2010-05-17 02:38:57 UTC (rev 31)
@@ -8,12 +8,12 @@
##' @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_model <- function(object,newdata=NULL,
- method="uni.mode",level=0.99, ...
+ method="sample.ml",level=0.99, ...
){
## Check and initialise parameters
stopifnot(is.null(newdata) || is.list(newdata))
- stopifnot(!"CI" %in% method || !is.null(level))
+ stopifnot(is.function(method) || !("CI" %in% method && is.null(level)))
###
## Fetch function and parameters from dream-model object
@@ -31,7 +31,7 @@
## Predict for desired method(s)
- if (method=="CI"){
+ if (!is.function(method) && method=="CI"){
sss <- window(object,...)
Modified: pkg/demo/FME.nonlinear.model.R
===================================================================
--- pkg/demo/FME.nonlinear.model.R 2010-05-13 07:55:08 UTC (rev 30)
+++ pkg/demo/FME.nonlinear.model.R 2010-05-17 02:38:57 UTC (rev 31)
@@ -48,14 +48,9 @@
pars <- list(p1=c(0,1),p2=c(0,100))
control <- list(
- nseq=4,
- thin.t=10
+ nseq=4
)
-## 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,
@@ -73,9 +68,15 @@
plot(dd)
plotFME()
-lines(Model(p=coef(dd),x=0:375),col="green")
+lines(predict(dd,
+ newdata=list(x=0:375)),
+ col="green")
+## Compare likelihood function for coefficients obtained by dream and FME modfit
+dd$lik.fun(coef(dd))
+dd$lik.fun(P$par)
+
########################
## Calculate bounds around output estimates
Modified: pkg/man/coef.dream.Rd
===================================================================
--- pkg/man/coef.dream.Rd 2010-05-13 07:55:08 UTC (rev 30)
+++ pkg/man/coef.dream.Rd 2010-05-17 02:38:57 UTC (rev 31)
@@ -3,7 +3,7 @@
\alias{coef.dream_model}
\title{Extract parameter values from dream or dream_model object}
\usage{
-\method{coef}{dream}(object, method = c("uni.mode", "mean", "median", "sample.ml"), \dots)
+\method{coef}{dream}(object, method = c("sample.ml","uni.mode", "mean", "median"), \dots)
}
\description{Extract parameter values using a choice of methods (or an
arbitrary function)}
@@ -17,6 +17,8 @@
mode of the univariate density estimate
for each parameter, using settings as in
\code{\link{densityplot.mcmc}}.
+ May not find the optimal parameter combination of the
+ distribution is multi-modal.
}
\item{\code{"mean"}}{
mean of each univariate parameter distribution.
Modified: pkg/man/dreamCalibrate.Rd
===================================================================
--- pkg/man/dreamCalibrate.Rd 2010-05-13 07:55:08 UTC (rev 30)
+++ pkg/man/dreamCalibrate.Rd 2010-05-17 02:38:57 UTC (rev 31)
@@ -50,7 +50,11 @@
\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}
+ \item{FUN.pars}{The extra arguments originally passed to that
+ function}
+ \item{lik.fun}{A convenience function f(pars) that returns the
+ log likelihood for a parameter set, as would be calculated using the
+ given calibration dataset and chosen likelihood function.}
}
\seealso{
Modified: pkg/man/predict.dream.Rd
===================================================================
--- pkg/man/predict.dream.Rd 2010-05-13 07:55:08 UTC (rev 30)
+++ pkg/man/predict.dream.Rd 2010-05-17 02:38:57 UTC (rev 31)
@@ -2,7 +2,7 @@
\alias{predict.dream_model}
\title{Predict values using dream_model object}
\usage{
-\S3method{predict}{dream}(object, newdata=NULL,method="uni.mode", level=0.99, \dots)
+\S3method{predict}{dream}(object, newdata=NULL,method="sample.ml", level=0.99, \dots)
}
\description{
Predict values using function calibrated by \code{\link{dreamCalibrate}}, optionally with new data,
Modified: pkg/man/window.dream.Rd
===================================================================
--- pkg/man/window.dream.Rd 2010-05-13 07:55:08 UTC (rev 30)
+++ pkg/man/window.dream.Rd 2010-05-17 02:38:57 UTC (rev 31)
@@ -1,5 +1,16 @@
\name{window.dream}
\alias{window.dream}
\title{Extract MCMC chains from a DREAM object}
+\usage{
+\method{window}{dream}(object, start = 1 + (end(object$Sequences) - 1) *
+(1 - fraction),fraction = 0.5,\dots)
+}
+\arguments{
+ \item{object}{dream object}
+ \item{start}{the first iteration of interest}
+ \item{fraction}{The fraction of the MCMC chains to keep}
+ \item{\dots}{extra arguments for \code{\link{window.mcmc}}}
+ }
+\value{\code{\link{mcmc.list}} object}
\description{Extract part or all of the MCMC chains from a DREAM object,
-specifying a burn-in period and allowing thinning}
+specifying a burn-in period and allowing thinning}
\ No newline at end of file
More information about the Dream-commits
mailing list