[Pomp-commits] r871 - in pkg/pomp: . R inst tests

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

Author: kingaa
Date: 2013-12-20 13:14:31 +0100 (Fri, 20 Dec 2013)
New Revision: 871

Modified:
   pkg/pomp/DESCRIPTION
   pkg/pomp/R/pmcmc.R
   pkg/pomp/inst/NEWS
   pkg/pomp/tests/ou2-pmcmc.R
   pkg/pomp/tests/ou2-pmcmc.Rout.save
Log:
- fix bug in pmcmc


Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION	2013-11-11 15:30:43 UTC (rev 870)
+++ pkg/pomp/DESCRIPTION	2013-12-20 12:14:31 UTC (rev 871)
@@ -1,8 +1,8 @@
 Package: pomp
 Type: Package
 Title: Statistical inference for partially observed Markov processes
-Version: 0.45-6
-Date: 2013-11-10
+Version: 0.45-7
+Date: 2013-12-20
 Authors at R: c(
 	   person(given=c("Aaron","A."),family="King",
 		role=c("aut","cre"),email="kingaa at umich.edu"),
@@ -29,7 +29,8 @@
 	 parmat.R logmeanexp.R slice-design.R 
 	 profile-design.R sobol.R bsplines.R sannbox.R
 	 pomp-fun.R pomp-class.R pomp.R pomp-methods.R 
-	 rmeasure-pomp.R rprocess-pomp.R init-state-pomp.R dmeasure-pomp.R dprocess-pomp.R skeleton-pomp.R 
+	 rmeasure-pomp.R rprocess-pomp.R init-state-pomp.R 
+	 dmeasure-pomp.R dprocess-pomp.R skeleton-pomp.R 
 	 simulate-pomp.R trajectory-pomp.R plot-pomp.R 
 	 pfilter.R pfilter-methods.R traj-match.R bsmc.R
 	 mif-class.R particles-mif.R mif.R mif-methods.R compare-mif.R 

Modified: pkg/pomp/R/pmcmc.R
===================================================================
--- pkg/pomp/R/pmcmc.R	2013-11-11 15:30:43 UTC (rev 870)
+++ pkg/pomp/R/pmcmc.R	2013-12-20 12:14:31 UTC (rev 871)
@@ -237,7 +237,7 @@
     gnsi <- FALSE
 
     ## PMCMC update rule (OK because proposal is symmetric)
-    if (runif(1) < exp(pfp.prop at loglik-pfp@loglik+log.prior-log.prior.prop)) {
+    if (runif(1) < exp(pfp.prop at loglik+log.prior.prop-pfp@loglik-log.prior)) {
       pfp <- pfp.prop
       theta <- theta.prop
     }

Modified: pkg/pomp/inst/NEWS
===================================================================
--- pkg/pomp/inst/NEWS	2013-11-11 15:30:43 UTC (rev 870)
+++ pkg/pomp/inst/NEWS	2013-12-20 12:14:31 UTC (rev 871)
@@ -1,7 +1,12 @@
 NEWS
-0.45-5
+0.45-7
+     o	bug fix in 'pmcmc': Metropolis-Hastings ratio was incorrect except for flat priors!
+	
+0.45-6
      o	a fix so that 'pompBuilder' will find 'pomp.h' header file on Windows machines.  Thanks to Dave Hayman for finding the problem.
 
+     o	new 'logmeanexp' function.
+
      o	some reorganization of the source package structure.
 
 0.45-4

Modified: pkg/pomp/tests/ou2-pmcmc.R
===================================================================
--- pkg/pomp/tests/ou2-pmcmc.R	2013-11-11 15:30:43 UTC (rev 870)
+++ pkg/pomp/tests/ou2-pmcmc.R	2013-12-20 12:14:31 UTC (rev 871)
@@ -58,4 +58,23 @@
   acf(conv.rec(f2)[,c("alpha.2","alpha.3")])
 }
 
+dprior.ou2 <- function (params, hyperparams, ..., log) {
+  f <- sum(dnorm(params,mean=hyperparams$mean,sd=hyperparams$sd,log=TRUE))
+  if (log) f else exp(f)
+}
+
+f5 <- pmcmc(
+            ou2,
+            start=coef(ou2),
+            Nmcmc=20,
+            dprior=dprior.ou2,
+            hyperparams=list(mean=coef(ou2),sd=1),
+            rw.sd=c(alpha.2=0.001,alpha.3=0.001),
+            Np=100,
+            max.fail=100, 
+            verbose=FALSE
+            )
+f5 <- continue(f5,Nmcmc=200,max.fail=100)
+plot(f1)
+
 dev.off()

Modified: pkg/pomp/tests/ou2-pmcmc.Rout.save
===================================================================
--- pkg/pomp/tests/ou2-pmcmc.Rout.save	2013-11-11 15:30:43 UTC (rev 870)
+++ pkg/pomp/tests/ou2-pmcmc.Rout.save	2013-12-20 12:14:31 UTC (rev 871)
@@ -1,7 +1,6 @@
 
-R version 2.15.2 (2012-10-26) -- "Trick or Treat"
-Copyright (C) 2012 The R Foundation for Statistical Computing
-ISBN 3-900051-07-0
+R version 3.0.2 (2013-09-25) -- "Frisbee Sailing"
+Copyright (C) 2013 The R Foundation for Statistical Computing
 Platform: x86_64-unknown-linux-gnu (64-bit)
 
 R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -81,10 +80,29 @@
 +   acf(conv.rec(f2)[,c("alpha.2","alpha.3")])
 + }
 > 
+> dprior.ou2 <- function (params, hyperparams, ..., log) {
++   f <- sum(dnorm(params,mean=hyperparams$mean,sd=hyperparams$sd,log=TRUE))
++   if (log) f else exp(f)
++ }
+> 
+> f5 <- pmcmc(
++             ou2,
++             start=coef(ou2),
++             Nmcmc=20,
++             dprior=dprior.ou2,
++             hyperparams=list(mean=coef(ou2),sd=1),
++             rw.sd=c(alpha.2=0.001,alpha.3=0.001),
++             Np=100,
++             max.fail=100, 
++             verbose=FALSE
++             )
+> f5 <- continue(f5,Nmcmc=200,max.fail=100)
+> plot(f1)
+> 
 > dev.off()
 null device 
           1 
 > 
 > proc.time()
    user  system elapsed 
-  6.788   0.052   6.997 
+ 11.480   0.036  11.720