[Pomp-commits] r979 - in pkg/pomp: R tests

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

Author: kingaa
Date: 2014-06-23 16:12:41 +0200 (Mon, 23 Jun 2014)
New Revision: 979

Modified:
   pkg/pomp/R/abc.R
   pkg/pomp/tests/abc.Rout.save
Log:
- fix bug with inadmissible proposals
- more efficient algorithm

Modified: pkg/pomp/R/abc.R
===================================================================
--- pkg/pomp/R/abc.R	2014-06-23 14:12:35 UTC (rev 978)
+++ pkg/pomp/R/abc.R	2014-06-23 14:12:41 UTC (rev 979)
@@ -95,9 +95,12 @@
   }
 
   theta <- start
-  log.prior <- dprior(object,params=theta,log=TRUE,
-                      .getnativesymbolinfo=gnsi)
-  ## we suppose that theta is a "match", which does the right thing for continue() and
+  log.prior <- dprior(object,params=theta,log=TRUE,.getnativesymbolinfo=gnsi)
+  if (!is.finite(log.prior))
+    stop("inadmissible value of ",sQuote("dprior"),
+         " at parameters ",sQuote("start"))
+  ## we suppose that theta is a "match",
+  ## which does the right thing for continue() and
   ## should have negligible effect unless doing many short calls to continue()
 
   conv.rec <- matrix(
@@ -132,38 +135,39 @@
   for (n in seq_len(Nabc)) { # main loop
 
     theta.prop <- theta
-    theta.prop[pars] <- rnorm(n=length(pars),mean=theta.prop[pars],sd=rw.sd)
+    theta.prop[pars] <- rnorm(n=length(pars),mean=theta[pars],sd=rw.sd)
+    log.prior.prop <- dprior(object,params=theta.prop,log=TRUE)
 
-    gnsi <- FALSE
+    if (is.finite(log.prior.prop) &&
+        runif(1) < exp(log.prior.prop-log.prior)) {
 
-    ## compute the probes for the proposed new parameter values
+      ## compute the probes for the proposed new parameter values
 
-    simval <- try(
-                  .Call(
-                        apply_probe_sim,
-                        object=object,
-                        nsim=1,
-                        params=theta.prop,
-                        seed=NULL,
-                        probes=probes,
-                        datval=datval
-                        ),
-                  silent=FALSE
-                  )
+      simval <- try(
+                    .Call(
+                          apply_probe_sim,
+                          object=object,
+                          nsim=1,
+                          params=theta.prop,
+                          seed=NULL,
+                          probes=probes,
+                          datval=datval
+                          ),
+                    silent=FALSE
+                    )
 
-    if (inherits(simval,'try-error'))
-      stop("abc error: error in ",sQuote("apply_probe_sim"),call.=FALSE)
+      if (inherits(simval,'try-error'))
+        stop("abc error: error in ",sQuote("apply_probe_sim"),call.=FALSE)
 
-    ## ABC update rule
-    distance <- sum(((datval-simval)/scale)^2)
-    if( (is.finite(distance)) && (distance<epssq) ){ 
-      log.prior.prop <- dprior(object,params=theta.prop,log=TRUE)
-      if (runif(1) < exp(log.prior.prop-log.prior)) {
+      ## ABC update rule
+      distance <- sum(((datval-simval)/scale)^2)
+      if( (is.finite(distance)) && (distance<epssq) ){ 
         theta <- theta.prop
         log.prior <- log.prior.prop
       }
+
     }
-
+    
     ## store a record of this iteration
     conv.rec[n+1,names(theta)] <- theta
     if (verbose && (n%%5==0))

Modified: pkg/pomp/tests/abc.Rout.save
===================================================================
--- pkg/pomp/tests/abc.Rout.save	2014-06-23 14:12:35 UTC (rev 978)
+++ pkg/pomp/tests/abc.Rout.save	2014-06-23 14:12:41 UTC (rev 979)
@@ -64,7 +64,7 @@
 > runs <- rle(as.vector(conv.rec(abc1)[, "alpha.1"]))
 > hist(runs$lengths)
 > mean(runs$length)
-[1] 7.31602
+[1] 6.949965
 > 
 > abc2 <- abc(po,
 +             Nabc=2000,
@@ -134,4 +134,4 @@
 > 
 > proc.time()
    user  system elapsed 
-  8.616   0.052   8.991 
+ 11.122   0.067  11.192