[Dream-commits] r11 - in pkg: . R demos
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Feb 26 00:44:36 CET 2010
Author: josephguillaume
Date: 2010-02-26 00:44:36 +0100 (Fri, 26 Feb 2010)
New Revision: 11
Added:
pkg/R/coef.dream.R
pkg/R/plot.dream.R
pkg/R/print.dream.R
pkg/R/summary.dream.R
Modified:
pkg/NAMESPACE
pkg/R/CompDensity.R
pkg/R/dream.R
pkg/R/maxLikPars.R
pkg/R/metrop.R
pkg/demos/FME nonlinear model.R
Log:
Added print,summary,coef,plot functions. Fixed errors so other func.types work.
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2010-02-25 06:57:41 UTC (rev 10)
+++ pkg/NAMESPACE 2010-02-25 23:44:36 UTC (rev 11)
@@ -4,3 +4,7 @@
export(CovInit)
export(LHSInit)
export(maxLikPars)
+S3method(summary,dream)
+S3method(coef,dream)
+S3method(plot,dream)
+S3method(print,dream)
\ No newline at end of file
Modified: pkg/R/CompDensity.R
===================================================================
--- pkg/R/CompDensity.R 2010-02-25 06:57:41 UTC (rev 10)
+++ pkg/R/CompDensity.R 2010-02-25 23:44:36 UTC (rev 11)
@@ -51,11 +51,8 @@
calc.loglik={
err <- measurement$data-modpred
- ## Compute the number of measurement data
- N <- length(measurement$data)
-
## Derive the log likelihood
- logp[ii] <- N*log(control$Wb/measurement$sigma)-
+ logp[ii] <- 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]
Added: pkg/R/coef.dream.R
===================================================================
--- pkg/R/coef.dream.R (rev 0)
+++ pkg/R/coef.dream.R 2010-02-25 23:44:36 UTC (rev 11)
@@ -0,0 +1,14 @@
+##' Extract maximum likelihood parameter values
+##' @param dream object
+##' @last.prop proportion of total sequence to use (0,1]
+##' if 1, use whole sequence
+##' @return named vector of parameter values
+## TODO: potentially use reduced.seq instead
+coef.dream <- function(dream.obj,last.prop=.5){
+ stopifnot(last.prop>0)
+ if (last.prop==1) return(maxLikPars(dream.obj$Sequences))
+ else {
+ ss <- window(dream.obj$Sequences, start = end(dd$Sequences)*(1-last.prop) + 1)
+ return(maxLikPars(ss))
+ }
+}##coef.dream
Property changes on: pkg/R/coef.dream.R
___________________________________________________________________
Name: svn:executable
+ *
Modified: pkg/R/dream.R
===================================================================
--- pkg/R/dream.R 2010-02-25 06:57:41 UTC (rev 10)
+++ pkg/R/dream.R 2010-02-25 23:44:36 UTC (rev 11)
@@ -19,10 +19,10 @@
## Parameters with auto-set values
ndim=NA, ## number of parameters (automatically set from length of pars)
DEpairs = NA, ## Number of DEpairs. defaults to max val floor((nseq-1)/2)
- nseq = NA, ## Number of Markov Chains / sequences (defaults to N)
+ nseq = NA ## Number of Markov Chains / sequences (defaults to N)
## Currently unused parameters
- trace = 0, ## level of user feedback
- REPORT = 10 ## number of iterations between reports when trace >= 1
+## trace = 0, ## level of user feedback
+## REPORT = 10 ## number of iterations between reports when trace >= 1
)
library(coda)
@@ -30,7 +30,8 @@
##' @param FUN model function with first argument a vector of parameter values of length ndim
-##' @param func.type. one of posterior.density, logposterior.density
+##' @param func.type type of value FUN returns.
+##' one of: posterior.density, logposterior.density,calc.loglik,calc.rmse,calc.weighted.rmse
##' @param pars a list of variable ranges
##' @param INIT f(pars,nseq,...) returns nseq x ndim matrix of initial parameter values
##' @param control
@@ -109,6 +110,11 @@
stop("unrecognised options: ",
toString(names(control)[!isValid]))
+ if (!is.null("measurement")){
+ if (! "sigma" %in% names(measurement)) measurement$sigma <- sd(measurement$data)
+ if (! "N" %in% names(measurement)) measurement$N <- length(measurement$data)
+ }
+
## determine number of variables to be optimized
control$ndim<-length(pars)
if (is.na(control$nseq)) control$nseq <- control$ndim
@@ -135,8 +141,8 @@
## Calculate the parameters in the exponential power density function of Box and Tiao (1973)
cbwb <- CalcCbWb(control$gamma)
- control$Cb <- cbwb$cb
- control$Wb <- cbwb$wb
+ control$Cb <- cbwb$Cb
+ control$Wb <- cbwb$Wb
## Generate the Table with JumpRates (dependent on number of dimensions and number of pairs)
Table.JumpRate<-matrix(NA,NDIM,control$DEpairs)
@@ -178,12 +184,16 @@
## Iter + AR at each step
obj$AR<-matrix(NA,n.elem,2)
obj$AR[1,2]<-NSEQ-1 ##Number if only one rejected
+ colnames(obj$AR) <- c("fun.evals","AR")
##Iter + R statistic for each variable at each step
obj$R.stat<-matrix(NA,ceiling(n.elem/control$steps),NDIM+1)
+ if (!is.null(names(pars))) colnames(obj$R.stat) <- c("fun.evals",names(pars))
+ else colnames(obj$R.stat) <- c("fun.evals",paste("p",1:length(pars),sep=""))
##Iter + pCR for each CR
obj$CR <- matrix(NA,ceiling(n.elem/control$steps),length(pCR)+1)
+ colnames(obj$CR) <- c("fun.evals",paste("CR",1:length(pCR),sep=""))
obj$outlier<-NULL
Modified: pkg/R/maxLikPars.R
===================================================================
--- pkg/R/maxLikPars.R 2010-02-25 06:57:41 UTC (rev 10)
+++ pkg/R/maxLikPars.R 2010-02-25 23:44:36 UTC (rev 11)
@@ -1,19 +1,20 @@
##' Naive maximum density parameter selection
##' @param ss MCMC chains. mcmc or mcmclist object
+##' @param ... extra parameters to pass to density
##' @return named character vector of parameter values
##'
-##' Uses which.max and density function with default parameters
-maxLikPars <- function(ss){
+##' Uses which.max and density function
+maxLikPars <- function(ss,...){
xx <- as.matrix(ss)
pars.maxp <- rep(NA,ncol(xx))
names(pars.maxp) <- colnames(xx)
maxp.res <- rep(NA,ncol(xx))
names(maxp.res) <- colnames(xx)
for (n in colnames(xx)){
- den <- density(xx[,n])
+ den <- density(xx[,n],...)
ii <- which.max(den$y)
pars.maxp[n] <- den$x[ii]
maxp.res[n] <- mean(diff(den$x))
}
return(pars.maxp)
-} ##maxLikPars
\ No newline at end of file
+} ##maxLikPars
Modified: pkg/R/metrop.R
===================================================================
--- pkg/R/metrop.R 2010-02-25 06:57:41 UTC (rev 10)
+++ pkg/R/metrop.R 2010-02-25 23:44:36 UTC (rev 11)
@@ -50,7 +50,7 @@
alpha <- pmin(exp(p.x-p.old),1)
},
calc.rmse = { ## SSE probability evaluation
- alpha <- pmin((p.x/p.old)^(-length(measurement$data)*(1+control$gamma)/2),1)
+ alpha <- pmin((p.x/p.old)^(-measurement$N*(1+control$gamma)/2),1)
},
logposterior.density = { ## Lnp probability evaluation
alpha <- pmin(exp(p.x-p.old),1)
@@ -59,7 +59,7 @@
## signs are different because we write -SSR
alpha <- pmin(exp(-0.5*(-p.x + p.old)/measurement$sigma^2),1);
},
- stop("Unrecognised value of control$metrop.opt")
+ stop("Unrecognised value of func.type")
)
## Generate random numbers
Added: pkg/R/plot.dream.R
===================================================================
--- pkg/R/plot.dream.R (rev 0)
+++ pkg/R/plot.dream.R 2010-02-25 23:44:36 UTC (rev 11)
@@ -0,0 +1,42 @@
+##' Plot characteristics of dream object
+##' @param dream object
+##' @param interactive - stop for each plot
+
+##' Uses second half of sequences
+
+plot.dream <- function(dream.obj,interactive=TRUE){
+ devAskNewPage(interactive)
+
+ ss <- window(dream.obj$Sequences, start = end(dd$Sequences)/2 + 1)
+
+ ## Trace and parameter density
+
+ plot(ss)
+
+ xyplot(ss)
+ densityplot(ss)
+
+ ## Acceptance rate
+ plot(table(dd$AR[,2]),main="Distribution of % acceptance rate")
+
+ ##Convergence
+
+ try(gelman.plot(ss))
+
+ plot(dd$R.stat[,1],dd$R.stat[,2],type="l",ylim=c(0,2))
+ for (i in 2:dd$control$ndim) lines(dd$R.stat[,1],dd$R.stat[,i+1],ylim=c(0,2))
+ title(main="Evolution of R.stat",sub="Equivalent to gelman.plot")
+
+ ## Multi-variate density for first chain
+
+ splom(as.data.frame(dd$Sequences[[1]]),
+ upper.panel = panel.smoothScatter, nrpoints = 0,
+ lower.panel = function(x, y, ...) {
+ panel.grid(-1, -1)
+ panel.loess(x, y, span = 1/3, lwd = 1)
+ panel.loess(x, y, span = 2/3, lwd = 2)
+ grid::grid.text(paste("cor =", round(cor(x, y),2)),
+ y = 0.1)
+ })
+
+}##plot.dream
Property changes on: pkg/R/plot.dream.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/print.dream.R
===================================================================
--- pkg/R/print.dream.R (rev 0)
+++ pkg/R/print.dream.R 2010-02-25 23:44:36 UTC (rev 11)
@@ -0,0 +1,18 @@
+print.dream <- function(dream.obj){
+ cat("
+Call:
+")
+
+ print(dream.obj$call)
+
+ cat("
+Control:
+")
+ for (i in names(dream.obj$control)){
+ v <- dream.obj$control[[i]]
+ cat(sprintf("%15s: ",i))
+ cat(v,"\n")
+ }
+
+ cat("\nExit condition:",dream.obj$EXITMSG,"\n")
+}
Property changes on: pkg/R/print.dream.R
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/R/summary.dream.R
===================================================================
--- pkg/R/summary.dream.R (rev 0)
+++ pkg/R/summary.dream.R 2010-02-25 23:44:36 UTC (rev 11)
@@ -0,0 +1,33 @@
+
+summary.dream <- function(dream.obj){
+
+ cat(sprintf("
+Exit message: %s
+Num fun evals: %d
+Time (secs): %.1f
+Final R.stats:
+",
+ dream.obj$EXITMSG,
+ dream.obj$fun.evals,
+ dream.obj$time
+ ))
+
+ R.stat.last <- tail(dream.obj$R.stat,1)
+ for (i in 2:ncol(dream.obj$R.stat)){
+ cat(sprintf("\t%s:\t%f\n",
+ colnames(dream.obj$R.stat)[i],
+ R.stat.last[i]
+ ))
+ } ##for
+
+ cat("
+CODA summary for last 50% of MCMC chains:
+")
+ print(summary(window(dream.obj$Sequences, start = end(dd$Sequences)/2 + 1)))
+
+ cat("
+Acceptance Rate
+")
+ summary(dream.obj$AR[,2])
+
+} ##summary.dream
Property changes on: pkg/R/summary.dream.R
___________________________________________________________________
Name: svn:executable
+ *
Modified: pkg/demos/FME nonlinear model.R
===================================================================
--- pkg/demos/FME nonlinear model.R 2010-02-25 06:57:41 UTC (rev 10)
+++ pkg/demos/FME nonlinear model.R 2010-02-25 23:44:36 UTC (rev 11)
@@ -16,12 +16,8 @@
y=c(0.05,0.07,0.09,0.10,0.11,0.122,0.125)) # 1/hour
obs.all <- rbind(Obs,Obs2)
-Model <- function(p,x) return(data.frame(x=x,y=p[1]*x/(x+p[2])))
-##Model(c(0.1,1),obs.all$x)
-
Model.y <- function(p,x) p[1]*x/(x+p[2])
-
control <- list(
nseq=5,
gamma=0,
@@ -33,7 +29,8 @@
pCR.Update=TRUE,
thin=TRUE,
thin.t=10,
- boundHandling="fold"
+ boundHandling="fold",
+ Rthres=1+1e-3
)
pars <- list(p1=c(0,1),p2=c(0,100))
@@ -49,28 +46,22 @@
control = control
)
-cat(sprintf("
-Exit message: %s
-Num fun evals: %d
-Time (secs): %.1f
-",
- dd$EXITMSG,
- dd$fun.evals,
- dd$time
- ))
-tail(dd$R.stat,1)
-maxLikPars(window(dd$Sequences, start = 0.75*end(dd$Sequences) + 1))
+dd
-ss <- window(dd$Sequences, start = end(dd$Sequences)/2 + 1)
-pars.maxp <- maxLikPars(ss)
-print(pars.maxp)
+summary(dd)
-plot(ss)
-gelman.plot(ss)
+coef(dd)
+plot(dd)
-### FME
+
+### Comparison to 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")
@@ -89,4 +80,4 @@
plot(Obs,xlab="mg COD/l",ylab="1/hour", pch=16, cex=1.5)
points(Obs2,pch=18,cex=1.5, col="red")
lines(Model(p=P$par,x=0:375))
-lines(Model(p=pars.maxp,x=0:375),col="green")
+lines(Model(p=coef(dd),x=0:375),col="green")
More information about the Dream-commits
mailing list