[Pomp-commits] r1174 - pkg

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

Author: kingaa
Date: 2015-06-04 14:22:34 +0200 (Thu, 04 Jun 2015)
New Revision: 1174

Modified:
   pkg/mif2.R
Log:
- mif2 -> pomp: stage 4 (ok?)

Modified: pkg/mif2.R
===================================================================
--- pkg/mif2.R	2015-06-04 12:22:32 UTC (rev 1173)
+++ pkg/mif2.R	2015-06-04 12:22:34 UTC (rev 1174)
@@ -525,7 +525,7 @@
           }
 )
 
-## mifList class
+## mif2List class
 setClass(
   'mif2List',
   contains='list',
@@ -606,6 +606,132 @@
   }
 )
 
+mif2.diagnostics <- function (z) {
+  ## assumes that z is a list of mif2d.pomps with identical structure
+  mar.multi <- c(0,5.1,0,2.1)
+  oma.multi <- c(6,0,5,0)
+  xx <- z[[1]]
+  ivpnames <- xx at ivps
+  estnames <- c(xx at pars,ivpnames)
+  parnames <- names(coef(xx,transform=xx at transform))
+  unestnames <- parnames[-match(estnames,parnames)]
+
+  ## plot filter means
+  filt.diag <- rbind("eff. sample size"=xx at eff.sample.size,filter.mean(xx))
+  filtnames <- rownames(filt.diag)
+  plotnames <- if(length(unestnames)>0) filtnames[-match(unestnames,filtnames)] else filtnames
+  lognames <- filtnames[1] # eff. sample size
+  nplots <- length(plotnames)
+  n.per.page <- min(nplots,10)
+  if(n.per.page<=4) nc <- 1 else nc <- 2
+  nr <- ceiling(n.per.page/nc)
+  oldpar <- par(mar=mar.multi,oma=oma.multi,mfcol=c(nr,nc),ask=dev.interactive(orNone=TRUE))
+  on.exit(par(oldpar))
+  low <- 1
+  hi <- 0
+  time <- time(xx)
+  while (hi<nplots) {
+    hi <- min(low+n.per.page-1,nplots)
+    for (i in seq(from=low,to=hi,by=1)) {
+      n <- i-low+1
+      logplot <- if (plotnames[i]%in%lognames) "y" else ""
+      dat <- sapply(
+                    z,
+                    function(po, label) {
+                      if (label=="eff. sample size")
+                        po at eff.sample.size
+                      else
+                        filter.mean(po,label)
+                    },
+                    label=plotnames[i]
+                    )
+      matplot(
+              y=dat,
+              x=time,
+              axes = FALSE,
+              xlab = "",
+              log=logplot,
+              ylab = "",
+              type = "l"
+              )
+      box()
+      y.side <- 2
+      axis(y.side, xpd = NA)
+      mtext(plotnames[i], y.side, line = 3)
+      do.xax <- (n%%nr==0||n==n.per.page)
+      if (do.xax) axis(1,xpd=NA)
+      if (do.xax) mtext("time",side=1,line=3)
+    }
+    low <- hi+1
+    mtext("Filter diagnostics (last iteration)",3,line=2,outer=TRUE)
+  }
+
+  ## plot mif convergence diagnostics
+  other.diagnostics <- c("loglik", "nfail")
+  plotnames <- c(other.diagnostics,estnames)
+  nplots <- length(plotnames)
+  n.per.page <- min(nplots,10)
+  nc <- if (n.per.page<=4) 1 else 2
+  nr <- ceiling(n.per.page/nc)
+  par(mar=mar.multi,oma=oma.multi,mfcol=c(nr,nc))
+  ## on.exit(par(oldpar))
+  low <- 1
+  hi <- 0
+  iteration <- seq(0,xx at Nmif)
+  while (hi<nplots) {
+    hi <- min(low+n.per.page-1,nplots)
+    for (i in seq(from=low,to=hi,by=1)) {
+      n <- i-low+1
+      dat <- sapply(z,function(po,label) conv.rec(po,label),label=plotnames[i])
+      matplot(
+              y=dat,
+              x=iteration,
+              axes = FALSE,
+              xlab = "",
+              ylab = "",
+              type = "l"
+              )
+      box()
+      y.side <- 2
+      axis(y.side,xpd=NA)
+      mtext(plotnames[i],y.side,line=3)
+      do.xax <- (n%%nr==0||n==n.per.page)
+      if (do.xax) axis(1,xpd=NA)
+      if (do.xax) mtext("MIF iteration",side=1,line=3)
+    }
+    low <- hi+1
+    mtext("MIF convergence diagnostics",3,line=2,outer=TRUE)
+  }
+  invisible(NULL)
+}
+
+
+setMethod(
+          "plot",
+          "mif2d.pomp",
+          function (x, y, ...) {
+            if (!missing(y)) {
+              y <- substitute(y)
+              warning(sQuote(y)," is ignored")
+            }
+            pomp:::mif.diagnostics(list(x))
+          }
+          )
+
+setMethod(
+          "plot",
+          signature=signature(x='mif2List'),
+          definition=function (x, y, ...) {
+            if (!missing(y)) {
+              y <- substitute(y)
+              warning(sQuote(y)," is ignored")
+            }
+            pomp:::mif.diagnostics(x)
+          }
+          )
+
+
+
 require(ggplot2)
 require(plyr)
 require(reshape2)
@@ -710,3 +836,6 @@
 
 ## ----first-mif-results-table,echo=FALSE,cache=FALSE----------------------
 print(results.table)
+
+# plot(do.call(c,lapply(mf,getElement,"mif")))
+