[Dream-commits] r28 - in pkg: R demo

[noreply at r-forge.r-project.org]

Author: josephguillaume
Date: 2010-05-12 07:58:05 +0200 (Wed, 12 May 2010)
New Revision: 28

Added:
   pkg/R/dreamCalibrate.R
Modified:
   pkg/demo/FME.nonlinear.model.R
Log:
Added preliminary dreamCalibrate function and split out each likelihood fn

Added: pkg/R/dreamCalibrate.R
===================================================================
--- pkg/R/dreamCalibrate.R	                        (rev 0)
+++ pkg/R/dreamCalibrate.R	2010-05-12 05:58:05 UTC (rev 28)
@@ -0,0 +1,56 @@
+## control requires gamma, (N)
+calc.rmse <- function(predicted,observed,control){
+  ## TODO: more general case with multidim data?
+  if (! "N" %in% names(control)) control$N <- length(observed)
+  
+  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)
+}## calc.rmse
+
+## TODO: calc.weighted.rmse
+
+## control: Wb,Cb,gamma,measurement.sigma
+calc.loglik <- function(predicted,observed,control){
+
+  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)
+
+  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
+}## calc.loglik
+
+## Design decisions:
+##  lik.fun is log.likelihood=f(predicted,observed,control)
+dreamCalibrate <- function(fun,
+                           pars,
+                           obs,
+                           lik.fun=calc.rmse,
+                           lik.control=list(),
+                           ... ##Extra arguments to dream
+                           ){
+  
+  FUN <- function(pars,...) lik.fun(fun(pars,...),obs,lik.control)
+  dd <- dream(FUN=fun,
+              pars=pars,
+              func.type="logposterior.density",
+              ... ##INIT,control,FUN.pars
+              )
+  
+  class(dd) <- c("dream-model",class(dd))
+  dd
+} ##dreamCalibrate


Property changes on: pkg/R/dreamCalibrate.R
___________________________________________________________________
Added: svn:executable
   + *

Modified: pkg/demo/FME.nonlinear.model.R
===================================================================
--- pkg/demo/FME.nonlinear.model.R	2010-05-12 04:50:28 UTC (rev 27)
+++ pkg/demo/FME.nonlinear.model.R	2010-05-12 05:58:05 UTC (rev 28)
@@ -50,11 +50,7 @@
 
 control <- list(
                 nseq=4,
-                thin.t=10,
-                parallel="none"
-                ##                REPORT=0
-                ##                ndraw=1000
-                ##                Rthres=1+1e-3
+                thin.t=10
                 )
 
 pars <- list(p1=c(0,1),p2=c(0,100))