[Pomp-commits] r826 - in pkg/pomp: . R man src tests

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Tue Feb 26 14:20:27 CET 2013


Author: kingaa
Date: 2013-02-26 14:20:27 +0100 (Tue, 26 Feb 2013)
New Revision: 826

Modified:
   pkg/pomp/DESCRIPTION
   pkg/pomp/NAMESPACE
   pkg/pomp/R/bsmc.R
   pkg/pomp/R/dmeasure-pomp.R
   pkg/pomp/R/dprocess-pomp.R
   pkg/pomp/R/init-state-pomp.R
   pkg/pomp/R/mif-methods.R
   pkg/pomp/R/mif.R
   pkg/pomp/R/particles-mif.R
   pkg/pomp/R/pfilter.R
   pkg/pomp/R/plot-pomp.R
   pkg/pomp/R/plugins.R
   pkg/pomp/R/pmcmc.R
   pkg/pomp/R/pomp-fun.R
   pkg/pomp/R/pomp-methods.R
   pkg/pomp/R/probe.R
   pkg/pomp/R/rmeasure-pomp.R
   pkg/pomp/R/rprocess-pomp.R
   pkg/pomp/R/simulate-pomp.R
   pkg/pomp/R/skeleton-pomp.R
   pkg/pomp/R/spect-match.R
   pkg/pomp/R/trajectory-pomp.R
   pkg/pomp/man/pmcmc.Rd
   pkg/pomp/man/probe.Rd
   pkg/pomp/src/SSA_wrapper.c
   pkg/pomp/src/dmeasure.c
   pkg/pomp/src/dprocess.c
   pkg/pomp/src/euler.c
   pkg/pomp/src/partrans.c
   pkg/pomp/src/pomp_fun.c
   pkg/pomp/src/pomp_internal.h
   pkg/pomp/src/rmeasure.c
   pkg/pomp/src/rprocess.c
   pkg/pomp/src/simulate.c
   pkg/pomp/src/skeleton.c
   pkg/pomp/src/trajectory.c
   pkg/pomp/tests/dimchecks.Rout.save
   pkg/pomp/tests/ou2-mif.Rout.save
   pkg/pomp/tests/ou2-mif2.Rout.save
   pkg/pomp/tests/ou2-pmcmc.R
   pkg/pomp/tests/ou2-pmcmc.Rout.save
   pkg/pomp/tests/ricker-probe.Rout.save
   pkg/pomp/tests/sir.Rout.save
Log:
- new mechanism to prevent unnecessary symbol-table lookups when native routines are used
- named functions are now used for most methods, to make error messages more sensible
- some changes to the implementation of pmcmc
- fix bug in pmcmc method on pmcmc objects


Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION	2013-02-05 22:49:05 UTC (rev 825)
+++ pkg/pomp/DESCRIPTION	2013-02-26 13:20:27 UTC (rev 826)
@@ -2,7 +2,7 @@
 Type: Package
 Title: Statistical inference for partially observed Markov processes
 Version: 0.44-1
-Date: 2013-02-04
+Date: 2013-02-14
 Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall, Helen Wearing, Matthew J. Ferrari, Michael Lavine, Daniel C. Reuman
 Maintainer: Aaron A. King <kingaa at umich.edu>
 URL: http://pomp.r-forge.r-project.org

Modified: pkg/pomp/NAMESPACE
===================================================================
--- pkg/pomp/NAMESPACE	2013-02-05 22:49:05 UTC (rev 825)
+++ pkg/pomp/NAMESPACE	2013-02-26 13:20:27 UTC (rev 826)
@@ -1,6 +1,5 @@
 useDynLib(			
           pomp,
-          get_pomp_fun,
           bspline_basis,
           periodic_bspline_basis,
           bspline_basis_function,

Modified: pkg/pomp/R/bsmc.R
===================================================================
--- pkg/pomp/R/bsmc.R	2013-02-05 22:49:05 UTC (rev 825)
+++ pkg/pomp/R/bsmc.R	2013-02-26 13:20:27 UTC (rev 826)
@@ -31,327 +31,368 @@
 
 setGeneric("bsmc",function(object,...)standardGeneric("bsmc"))
 
-setMethod(
-          "bsmc",
-          "pomp",
-          function (object, params, Np, est,
-                    smooth = 0.1,
-                    ntries = 1,
-                    tol = 1e-17,
-                    lower = -Inf, upper = Inf,
-                    seed = NULL,
-                    verbose = getOption("verbose"),
-                    max.fail = 0,
-                    transform = FALSE,
-                    ...) {
+bsmc.internal <- function (object, params, Np, est,
+                           smooth = 0.1,
+                           ntries = 1,
+                           tol = 1e-17,
+                           lower = -Inf, upper = Inf,
+                           seed = NULL,
+                           verbose = getOption("verbose"),
+                           max.fail = 0,
+                           transform = FALSE,
+                           .getnativesymbolinfo = TRUE,
+                           ...) {
 
-            transform <- as.logical(transform)
+  gnsi.rproc <- gnsi.dmeas <- as.logical(.getnativesymbolinfo)
+  ptsi.inv <- ptsi.for <- TRUE
+  transform <- as.logical(transform)
 
-            if (missing(seed)) seed <- NULL
-            if (!is.null(seed)) {
-              if (!exists(".Random.seed",where=.GlobalEnv))
-                runif(n=1L) ## need to initialize the RNG
-              save.seed <- get(".Random.seed",pos=.GlobalEnv)
-              set.seed(seed)
-            }
+  if (missing(seed)) seed <- NULL
+  if (!is.null(seed)) {
+    if (!exists(".Random.seed",where=.GlobalEnv))
+      runif(n=1L) ## need to initialize the RNG
+    save.seed <- get(".Random.seed",pos=.GlobalEnv)
+    set.seed(seed)
+  }
 
-            error.prefix <- paste(sQuote("bsmc"),"error: ")
+  error.prefix <- paste(sQuote("bsmc"),"error: ")
 
-            if (missing(params)) {
-              if (length(coef(object))>0) {
-                params <- coef(object)
-              } else {
-                stop(error.prefix,sQuote("params")," must be supplied",call.=FALSE)
-              }
-            }
+  if (missing(params)) {
+    if (length(coef(object))>0) {
+      params <- coef(object)
+    } else {
+      stop(error.prefix,sQuote("params")," must be supplied",call.=FALSE)
+    }
+  }
 
-            if (missing(Np)) Np <- NCOL(params)
-            else if (is.matrix(params)&&(Np!=ncol(params)))
-              warning(sQuote("Np")," is ignored when ",sQuote("params")," is a matrix")
-            
-            if (transform)
-              params <- partrans(object,params,dir="inverse")
+  if (missing(Np)) Np <- NCOL(params)
+  else if (is.matrix(params)&&(Np!=ncol(params)))
+    warning(sQuote("Np")," is ignored when ",sQuote("params")," is a matrix")
+  
+  if (transform)
+    params <- partrans(object,params,dir="inverse",
+                       .getnativesymbolinfo=ptsi.inv)
+  ptsi.inv <- FALSE
+  
+  ntimes <- length(time(object))
+  if (is.null(dim(params))) {
+    params <- matrix(
+                     params,
+                     nrow=length(params),
+                     ncol=Np,
+                     dimnames=list(
+                       names(params),
+                       NULL
+                       )
+                     )
+  }
 
-            ntimes <- length(time(object))
-            if (is.null(dim(params))) {
-              params <- matrix(
-                               params,
-                               nrow=length(params),
-                               ncol=Np,
-                               dimnames=list(
-                                 names(params),
-                                 NULL
-                                 )
-                               )
-            }
+  npars <- nrow(params)
+  paramnames <- rownames(params)
+  prior <- params
 
-            npars <- nrow(params)
-            paramnames <- rownames(params)
-            prior <- params
+  if (missing(est))
+    est <- paramnames[apply(params,1,function(x)diff(range(x))>0)]
+  estind <- match(est,paramnames)
+  npars.est <- length(estind)
+  
+  if (npars.est<1)
+    stop(error.prefix,"no parameters to estimate",call.=FALSE)
 
-            if (missing(est))
-              est <- paramnames[apply(params,1,function(x)diff(range(x))>0)]
-            estind <- match(est,paramnames)
-            npars.est <- length(estind)
-            
-            if (npars.est<1)
-              stop(error.prefix,"no parameters to estimate",call.=FALSE)
+  if (is.null(paramnames))
+    stop(error.prefix,sQuote("params")," must have rownames",call.=FALSE)
 
-            if (is.null(paramnames))
-              stop(error.prefix,sQuote("params")," must have rownames",call.=FALSE)
+  if ((length(smooth)!=1)||(smooth>1)||(smooth<=0))
+    stop(error.prefix,sQuote("smooth")," must be a scalar in [0,1)",call.=FALSE)
 
-            if ((length(smooth)!=1)||(smooth>1)||(smooth<=0))
-              stop(error.prefix,sQuote("smooth")," must be a scalar in [0,1)",call.=FALSE)
+  hsq <- smooth^2             #  see Liu & West eq(3.6) p10
+  shrink <- sqrt(1-hsq)
 
-            hsq <- smooth^2             #  see Liu & West eq(3.6) p10
-            shrink <- sqrt(1-hsq)
+  if (
+      ((length(lower)>1)&&(length(lower)!=npars.est))||
+      ((length(upper)>1)&&(length(upper)!=npars.est))
+      ) {
+    stop(
+         error.prefix,
+         sQuote("lower")," and ",sQuote("upper"),
+         " must each have length 1 or length equal to that of ",sQuote("est"),
+         call.=FALSE                   
+         )
+  }
 
-            if (
-                ((length(lower)>1)&&(length(lower)!=npars.est))||
-                ((length(upper)>1)&&(length(upper)!=npars.est))
-                ) {
-              stop(
-                   error.prefix,
-                   sQuote("lower")," and ",sQuote("upper"),
-                   " must each have length 1 or length equal to that of ",sQuote("est"),
-                   call.=FALSE                   
-                   )
-            }
+  for (j in seq_len(Np)) {
+    if (any((params[estind,j]<lower)|(params[estind,j]>upper))) {
+      ind <- which((params[estind,j]<lower)|(params[estind,j]>upper))
+      stop(
+           error.prefix,
+           "parameter(s) ",paste(paramnames[estind[ind]],collapse=","),
+           " in column ",j," in ",sQuote("params"),
+           " is/are outside the box defined by ",
+           sQuote("lower")," and ",sQuote("upper"),
+           call.=FALSE
+           )
+    }
+  }
 
-            for (j in seq_len(Np)) {
-              if (any((params[estind,j]<lower)|(params[estind,j]>upper))) {
-                ind <- which((params[estind,j]<lower)|(params[estind,j]>upper))
-                stop(
-                     error.prefix,
-                     "parameter(s) ",paste(paramnames[estind[ind]],collapse=","),
-                     " in column ",j," in ",sQuote("params"),
-                     " is/are outside the box defined by ",
-                     sQuote("lower")," and ",sQuote("upper"),
-                     call.=FALSE
-                     )
-              }
-            }
+  xstart <- init.state(
+                       object,
+                       params=if (transform) {
+                         partrans(object,params,dir="forward",
+                                  .getnativesymbolinfo=ptsi.for)
+                       } else {
+                         params
+                       }
+                       )
+  statenames <- rownames(xstart)
+  nvars <- nrow(xstart)
+  ptsi.for <- FALSE
+  
+  times <- time(object,t0=TRUE)
+  x <- xstart
 
-            xstart <- init.state(
-                                 object,
-                                 params=if (transform) {
-                                   partrans(object,params,dir="forward")
-                                 } else {
-                                   params
-                                 }
-                                 )
-            statenames <- rownames(xstart)
-            nvars <- nrow(xstart)
-            
-            times <- time(object,t0=TRUE)
-            x <- xstart
+  evidence <- rep(NA,ntimes)
+  eff.sample.size <- rep(NA,ntimes)
+  nfail <- 0
+  
+  mu <- array(data=NA,dim=c(nvars,Np,1)) 
+  rownames(mu) <- rownames(xstart)
+  m  <- array(data=NA,dim=c(npars,Np))
+  rownames(m) <- rownames(params)
+  
+  for (nt in seq_len(ntimes)) {
+    
+    ## calculate particle means ; as per L&W AGM (1)
+    params.mean <- apply(params,1,mean) 
+    ## calculate particle covariances : as per L&W AGM (1)
+    params.var  <- cov(t(params[estind,,drop=FALSE]))
 
-            evidence <- rep(NA,ntimes)
-            eff.sample.size <- rep(NA,ntimes)
-            nfail <- 0
-            
-            mu <- array(data=NA,dim=c(nvars,Np,1)) 
-            rownames(mu) <- rownames(xstart)
-            m  <- array(data=NA,dim=c(npars,Np))
-            rownames(m) <- rownames(params)
-            
-            for (nt in seq_len(ntimes)) {
-              
-              ## calculate particle means ; as per L&W AGM (1)
-              params.mean <- apply(params,1,mean) 
-              ## calculate particle covariances : as per L&W AGM (1)
-              params.var  <- cov(t(params[estind,,drop=FALSE]))
+    if (verbose) {
+      cat("at step",nt,"(time =",times[nt+1],")\n")
+      print(
+            rbind(
+                  prior.mean=params.mean[estind],
+                  prior.sd=sqrt(diag(params.var))
+                  )
+            )
+    }
 
-              if (verbose) {
-                cat("at step",nt,"(time =",times[nt+1],")\n")
-                print(
-                      rbind(
-                            prior.mean=params.mean[estind],
-                            prior.sd=sqrt(diag(params.var))
-                            )
+    ## update mean of states at time nt as per L&W AGM (1) 
+    tries <- rprocess(
+                      object,
+                      xstart=parmat(x,nrep=ntries),
+                      times=times[c(nt,nt+1)],
+                      params=if (transform) {
+                        partrans(object,params,dir="forward",
+                                 .getnativesymbolinfo=ptsi.for)
+                      } else {
+                        params
+                      },
+                      offset=1,
+                      .getnativesymbolinfo=gnsi.rproc
                       )
-              }
+    dim(tries) <- c(nvars,Np,ntries,1)
+    mu <- apply(tries,c(1,2,4),mean)
+    rownames(mu) <- statenames
+    ## shrink parameters towards mean as per Liu & West eq (3.3) and L&W AGM (1)
+    m <- shrink*params+(1-shrink)*params.mean
+    gnsi.rproc <- FALSE
+    
+    ## evaluate probability of obervation given mean value of parameters and states (used in L&W AGM (5) below)
+    g <- dmeasure( 
+                  object,
+                  y=object at data[,nt,drop=FALSE],
+                  x=mu,
+                  times=times[nt+1],
+                  params=if (transform) {
+                    partrans(object,m,dir="forward",
+                             .getnativesymbolinfo=ptsi.for)
+                  } else {
+                    m
+                  },
+                  .getnativesymbolinfo=gnsi.dmeas
+                  )	
+    gnsi.dmeas <- FALSE
+    storeForEvidence1 <- log(sum(g))
+    ## sample indices -- From L&W AGM (2)
+    ##              k <- .Call(systematic_resampling,g)
+    k <- sample.int(n=Np,size=Np,replace=TRUE,prob=g)
+    params <- params[,k]
+    m <- m[,k]
+    g <- g[k]
 
-              ## update mean of states at time nt as per L&W AGM (1) 
-              tries <- rprocess(
-                                object,
-                                xstart=parmat(x,nrep=ntries),
-                                times=times[c(nt,nt+1)],
-                                params=if (transform) {
-                                  partrans(object,params,dir="forward")
-                                } else {
-                                  params
-                                },
-                                offset=1
-                                )
-              dim(tries) <- c(nvars,Np,ntries,1)
-              mu <- apply(tries,c(1,2,4),mean)
-              rownames(mu) <- statenames
-              ## shrink parameters towards mean as per Liu & West eq (3.3) and L&W AGM (1)
-              m <- shrink*params+(1-shrink)*params.mean
-              
-              ## evaluate probability of obervation given mean value of parameters and states (used in L&W AGM (5) below)
-              g <- dmeasure( 
-                            object,
-                            y=object at data[,nt,drop=FALSE],
-                            x=mu,
-                            times=times[nt+1],
-                            params=if (transform) {
-                              partrans(object,m,dir="forward")
-                            } else {
-                              m
-                            }
-                            )	
-              storeForEvidence1 <- log(sum(g))
-              ## sample indices -- From L&W AGM (2)
-##              k <- .Call(systematic_resampling,g)
-              k <- sample.int(n=Np,size=Np,replace=TRUE,prob=g)
-              params <- params[,k]
-              m <- m[,k]
-              g <- g[k]
+    ## sample new parameter vector as per L&W AGM (3) and Liu & West eq(3.2)
+    pvec <- try(
+                mvtnorm::rmvnorm(
+                                 n=Np,
+                                 mean=rep(0,npars.est),
+                                 sigma=hsq*params.var,
+                                 method="svd"
+                                 ),
+                silent=FALSE
+                )
+    if (inherits(pvec,"try-error"))
+      stop(error.prefix,"error in ",sQuote("rmvnorm"),call.=FALSE)
+    if (any(!is.finite(pvec)))
+      stop(error.prefix,"extreme particle depletion",call.=FALSE)
+    params[estind,] <- m[estind,]+t(pvec)
 
-              ## sample new parameter vector as per L&W AGM (3) and Liu & West eq(3.2)
-              pvec <- try(
-                          mvtnorm::rmvnorm(
-                                           n=Np,
-                                           mean=rep(0,npars.est),
-                                           sigma=hsq*params.var,
-                                           method="svd"
-                                           ),
-                          silent=FALSE
-                          )
-              if (inherits(pvec,"try-error"))
-                stop(error.prefix,"error in ",sQuote("rmvnorm"),call.=FALSE)
-              if (any(!is.finite(pvec)))
-                stop(error.prefix,"extreme particle depletion",call.=FALSE)
-              params[estind,] <- m[estind,]+t(pvec)
+    if (transform)
+      tparams <- partrans(object,params,dir="forward",
+                          .getnativesymbolinfo=ptsi.for)
+    
+    ## sample current state vector x^(g)_(t+1) as per L&W AGM (4)
+    X <- rprocess(
+                  object,
+                  xstart=x[,k,drop=FALSE],
+                  times=times[c(nt,nt+1)],
+                  params=if (transform) {
+                    tparams
+                  } else {
+                    params
+                  },
+                  offset=1,
+                  .getnativesymbolinfo=gnsi.rproc
+                  )
 
-              if (transform)
-                tparams <- partrans(object,params,dir="forward")
-              
-              ## sample current state vector x^(g)_(t+1) as per L&W AGM (4)
-              X <- rprocess(
-                            object,
-                            xstart=x[,k,drop=FALSE],
-                            times=times[c(nt,nt+1)],
-                            params=if (transform) {
-                              tparams
-                            } else {
-                              params
-                            },
-                            offset=1
-                            )
+    ## evaluate likelihood of observation given X (from L&W AGM (4))
+    numer <- dmeasure(
+                      object,
+                      y=object at data[,nt,drop=FALSE],
+                      x=X,
+                      times=times[nt+1],
+                      params=if (transform) {
+                        tparams
+                      } else {
+                        params
+                      },
+                      .getnativesymbolinfo=gnsi.dmeas
+                      )
+    ## evaluate weights as per L&W AGM (5)
 
-              ## evaluate likelihood of observation given X (from L&W AGM (4))
-              numer <- dmeasure(
-                                object,
-                                y=object at data[,nt,drop=FALSE],
-                                x=X,
-                                times=times[nt+1],
-                                params=if (transform) {
-                                  tparams
-                                } else {
-                                  params
-                                }
-                                )
-              ## evaluate weights as per L&W AGM (5)
+    weights <- numer/g
+    storeForEvidence2 <- log(mean(weights))
+    
+    ## apply box constraints as per the priors          
+    for (j in seq_len(Np)) {
+      ## the following seems problematic: will it tend to make the boundaries repellors
+      if (any((params[estind,j]<lower)|(params[estind,j]>upper))) {
+        weights[j] <- 0 
+      }
+      ## might this rejection method be preferable?
+      ## while (any((params[estind,j]<lower)|(params[estind,j]>upper))) {
+      ##   ## rejection method
+      ##   pvec <- try(
+      ##               mvtnorm::rmvnorm(
+      ##                                n=1,
+      ##                                mean=rep(0,npars.est),
+      ##                                sigma=hsq*params.var,
+      ##                                method="eigen"
+      ##                                ),
+      ##               silent=FALSE
+      ##               )
+      ##   if (inherits(pvec,"try-error"))
+      ##     stop(error.prefix,"error in ",sQuote("rmvnorm"),call.=FALSE)
+      ##   if (any(!is.finite(pvec)))
+      ##     stop(error.prefix,"extreme particle depletion",call.=FALSE)
+      ##   params[estind,j] <- m[estind,j]+pvec[1,]
+      ## }
+    }
 
-	      weights <- numer/g
-	      storeForEvidence2 <- log(mean(weights))
-              
-              ## apply box constraints as per the priors          
-              for (j in seq_len(Np)) {
-                ## the following seems problematic: will it tend to make the boundaries repellors
-                if (any((params[estind,j]<lower)|(params[estind,j]>upper))) {
-                  weights[j] <- 0 
-                }
-                ## might this rejection method be preferable?
-                ## while (any((params[estind,j]<lower)|(params[estind,j]>upper))) {
-                ##   ## rejection method
-                ##   pvec <- try(
-                ##               mvtnorm::rmvnorm(
-                ##                                n=1,
-                ##                                mean=rep(0,npars.est),
-                ##                                sigma=hsq*params.var,
-                ##                                method="eigen"
-                ##                                ),
-                ##               silent=FALSE
-                ##               )
-                ##   if (inherits(pvec,"try-error"))
-                ##     stop(error.prefix,"error in ",sQuote("rmvnorm"),call.=FALSE)
-                ##   if (any(!is.finite(pvec)))
-                ##     stop(error.prefix,"extreme particle depletion",call.=FALSE)
-                ##   params[estind,j] <- m[estind,j]+pvec[1,]
-                ## }
-              }
+    x[,] <- X                
+    
+    ## test for failure to filter
+    dim(weights) <- NULL
+    failures <- ((weights<tol)|(!is.finite(weights))) # test for NA weights
+    all.fail <- all(failures)
+    if (all.fail) {                     # all particles are lost
+      if (verbose) {
+        message("filtering failure at time t = ",times[nt+1])
+      }
+      nfail <- nfail+1
+      if (nfail > max.fail)
+        stop(error.prefix,"too many filtering failures",call.=FALSE)
+      evidence[nt] <- log(tol)          # worst log-likelihood
+      weights <- rep(1/Np,Np)
+      eff.sample.size[nt] <- 0
+    } else {                  # not all particles are lost
+      ## compute log-likelihood
+      evidence[nt] <- storeForEvidence1+storeForEvidence2
+      weights[failures] <- 0
+      weights <- weights/sum(weights)
+      ## compute effective sample-size
+      eff.sample.size[nt] <- 1/crossprod(weights)
+    }
 
-              x[,] <- X                
-              
-              ## test for failure to filter
-              dim(weights) <- NULL
-              failures <- ((weights<tol)|(!is.finite(weights))) # test for NA weights
-              all.fail <- all(failures)
-              if (all.fail) {                     # all particles are lost
-                if (verbose) {
-                  message("filtering failure at time t = ",times[nt+1])
-                }
-                nfail <- nfail+1
-                if (nfail > max.fail)
-                  stop(error.prefix,"too many filtering failures",call.=FALSE)
-                evidence[nt] <- log(tol)          # worst log-likelihood
-                weights <- rep(1/Np,Np)
-                eff.sample.size[nt] <- 0
-              } else {                  # not all particles are lost
-                ## compute log-likelihood
-                evidence[nt] <- storeForEvidence1+storeForEvidence2
-                weights[failures] <- 0
-                weights <- weights/sum(weights)
-                ## compute effective sample-size
-                eff.sample.size[nt] <- 1/crossprod(weights)
-              }
+    if (verbose) {
+      cat("effective sample size =",round(eff.sample.size[nt],1),"\n")
+    }
 
-              if (verbose) {
-                cat("effective sample size =",round(eff.sample.size[nt],1),"\n")
-              }
+    ## Matrix with samples (columns) from filtering distribution theta.t | Y.t
+    if (!all.fail) {
+      ## smp <- .Call(systematic_resampling,weights)
+      smp <- sample.int(n=Np,size=Np,replace=TRUE,prob=weights)
+      x <- x[,smp,drop=FALSE]
+      params[estind,] <- params[estind,smp,drop=FALSE]
+    }
 
-              ## Matrix with samples (columns) from filtering distribution theta.t | Y.t
-              if (!all.fail) {
-                ## smp <- .Call(systematic_resampling,weights)
-                smp <- sample.int(n=Np,size=Np,replace=TRUE,prob=weights)
-                x <- x[,smp,drop=FALSE]
-                params[estind,] <- params[estind,smp,drop=FALSE]
-              }
-              
-            }
-            
-            if (!is.null(seed)) {
-              assign(".Random.seed",save.seed,pos=.GlobalEnv)
-              seed <- save.seed
-            }
-            
-            ## if (transform) {
-            ##   params <- partrans(object,params,dir="forward")
-            ##   prior <- partrans(object,prior,dir="forward")
-            ## }
+    .getnativesymbolinfo <- FALSE
+    
+  }
 
-            ## replace parameters with point estimate (posterior median)
-            coef(object,transform=transform) <- apply(params,1,median)
+  if (!is.null(seed)) {
+    assign(".Random.seed",save.seed,pos=.GlobalEnv)
+    seed <- save.seed
+  }
+  
+  ## replace parameters with point estimate (posterior median)
+  coef(object,transform=transform) <- apply(params,1,median)
 
-            new(
-                "bsmcd.pomp",
-                object,
-                transform=transform,
-                post=params,
-                prior=prior,
-                est=as.character(est),
-                eff.sample.size=eff.sample.size,
-                smooth=smooth,
-                seed=as.integer(seed),
-                nfail=as.integer(nfail),
-                cond.log.evidence=evidence,
-                log.evidence=sum(evidence),
-                weights=weights
-                )
+  new(
+      "bsmcd.pomp",
+      object,
+      transform=transform,
+      post=params,
+      prior=prior,
+      est=as.character(est),
+      eff.sample.size=eff.sample.size,
+      smooth=smooth,
+      seed=as.integer(seed),
+      nfail=as.integer(nfail),
+      cond.log.evidence=evidence,
+      log.evidence=sum(evidence),
+      weights=weights
+      )
+}
+
+setMethod(
+          "bsmc",
+          "pomp",
+          function (object, params, Np, est,
+                    smooth = 0.1,
+                    ntries = 1,
+                    tol = 1e-17,
+                    lower = -Inf, upper = Inf,
+                    seed = NULL,
+                    verbose = getOption("verbose"),
+                    max.fail = 0,
+                    transform = FALSE,
+                    ...) {
+            bsmc.internal(
+                          object=object,
+                          params=params,
+                          Np=Np,
+                          est=est,
+                          smooth=smooth,
+                          ntries=ntries,
+                          tol=tol,
+                          lower=lower,
+                          upper=upper,
+                          seed=seed,
+                          verbose=verbose,
+                          max.fail=max.fail,
+                          transform=transform,
+                          ...
+                          )
           }
           )
 

Modified: pkg/pomp/R/dmeasure-pomp.R
===================================================================
--- pkg/pomp/R/dmeasure-pomp.R	2013-02-05 22:49:05 UTC (rev 825)
+++ pkg/pomp/R/dmeasure-pomp.R	2013-02-26 13:20:27 UTC (rev 826)
@@ -1,9 +1,11 @@
 ## evaluate the measurement model density function
 setGeneric("dmeasure",function(object,...)standardGeneric("dmeasure"))
 
-dmeasure.internal <- function (object, y, x, times, params, log = FALSE, ...) {
-  fun <- get.pomp.fun(object at dmeasure)
-  .Call(do_dmeasure,object,y,x,times,params,log,fun)
+dmeasure.internal <- function (object, y, x, times, params, log = FALSE, .getnativesymbolinfo = TRUE, ...) {
+  .Call(do_dmeasure,object,y,x,times,params,log,.getnativesymbolinfo)
 }
 
-setMethod("dmeasure","pomp",dmeasure.internal)
+setMethod("dmeasure","pomp",
+          function (object, y, x, times, params, log = FALSE, ...)
+          dmeasure.internal(object=object,y=y,x=x,times=times,params=params,log=log,...)
+          )

Modified: pkg/pomp/R/dprocess-pomp.R
===================================================================
--- pkg/pomp/R/dprocess-pomp.R	2013-02-05 22:49:05 UTC (rev 825)
+++ pkg/pomp/R/dprocess-pomp.R	2013-02-26 13:20:27 UTC (rev 826)
@@ -1,7 +1,10 @@
 ## evaluate the process model density function
 setGeneric("dprocess",function(object,...)standardGeneric("dprocess"))
 
-dprocess.internal <- function (object, x, times, params, log = FALSE, ...)
-  .Call(do_dprocess,object,x,times,params,log)
+dprocess.internal <- function (object, x, times, params, log = FALSE, .getnativesymbolinfo = TRUE, ...)
+  .Call(do_dprocess,object,x,times,params,log,.getnativesymbolinfo)
 
-setMethod("dprocess","pomp",dprocess.internal)
+setMethod("dprocess","pomp",
+          function (object, x, times, params, log = FALSE, ...)
+          dprocess.internal(object=object,x=x,times=times,params=params,log=log,...)
+          )

Modified: pkg/pomp/R/init-state-pomp.R
===================================================================
--- pkg/pomp/R/init-state-pomp.R	2013-02-05 22:49:05 UTC (rev 825)
+++ pkg/pomp/R/init-state-pomp.R	2013-02-26 13:20:27 UTC (rev 826)
@@ -1,12 +1,16 @@
+## initialize the state variables of the process model
+init.state.internal <- function (object, params, t0, ...) {
+  if (missing(t0)) t0 <- object at t0
+  if (missing(params)) params <- coef(object)
+  .Call(do_init_state,object,params,t0)
+}
+
 setGeneric("init.state",function(object,...)standardGeneric("init.state"))
 
-## initialize the process model
 setMethod(
           'init.state',
           'pomp',
-          function (object, params, t0, ...) { # the package algorithms will only use these arguments
-            if (missing(t0)) t0 <- object at t0
-            if (missing(params)) params <- coef(object)
-            .Call(do_init_state,object,params,t0)
+          function (object, params, t0, ...) {
+            init.state.internal(object=object,params=params,t0=t0,...)
           }
           )

Modified: pkg/pomp/R/mif-methods.R
===================================================================
--- pkg/pomp/R/mif-methods.R	2013-02-05 22:49:05 UTC (rev 825)
+++ pkg/pomp/R/mif-methods.R	2013-02-26 13:20:27 UTC (rev 826)
@@ -1,48 +1,48 @@
 ## this file contains short definitions of methods for the 'mif' class
 
-conv.rec <- function (object, ...)
-  stop("function ",sQuote("conv.rec")," is undefined for objects of class ",sQuote(class(object)))
-setGeneric('conv.rec')  
-
 ## extract the estimated log likelihood
 setMethod('logLik','mif',function(object,...)object at loglik)
 
 ## extract the convergence record
-setMethod(
-          'conv.rec',
-          'mif',
+conv.rec.internal <- function (object, pars, transform = FALSE, ...) {
+  if (transform) {
+    pars.improper <- c("loglik","nfail")
+    pars.proper <- setdiff(colnames(object at conv.rec),pars.improper)
+    retval <- cbind(
+                    t(
+                      partrans(
+                               object,
+                               params=t(object at conv.rec[,pars.proper]),
+                               dir="forward"
+                               )
+                      ),
+                    object at conv.rec[,pars.improper]
+                    )
+  } else {
+    retval <- object at conv.rec
+  }
+  if (missing(pars))
+    retval
+  else {
+    bad.pars <- setdiff(pars,colnames(retval))
+    if (length(bad.pars)>0)
+      stop(
+           "in ",sQuote("conv.rec"),": name(s) ",
+           paste(sQuote(bad.pars),collapse=","),
+           " correspond to no parameter(s) in ",
+           if (transform) sQuote("conv.rec(object,transform=TRUE)")
+           else sQuote("conv.rec(object,transform=FALSE)"),
+           call.=FALSE
+           )
+    retval[,pars]
+  }
+}
+
+setGeneric("conv.rec",function(object,...)standardGeneric("conv.rec"))
+
+setMethod('conv.rec','mif',
           function (object, pars, transform = FALSE, ...) {
-            if (transform) {
-              pars.improper <- c("loglik","nfail")
-              pars.proper <- setdiff(colnames(object at conv.rec),pars.improper)
-              retval <- cbind(
-                              t(
-                                partrans(
-                                         object,
-                                         params=t(object at conv.rec[,pars.proper]),
-                                         dir="forward"
-                                         )
-                                ),
-                              object at conv.rec[,pars.improper]
-                              )
-            } else {
-              retval <- object at conv.rec
-            }
-            if (missing(pars))
-              retval
-            else {
-              bad.pars <- setdiff(pars,colnames(retval))
-              if (length(bad.pars)>0)
-                stop(
-                     "in ",sQuote("conv.rec"),": name(s) ",
-                     paste(sQuote(bad.pars),collapse=","),
-                     " correspond to no parameter(s) in ",
-                     if (transform) sQuote("conv.rec(object,transform=TRUE)")
-                     else sQuote("conv.rec(object,transform=FALSE)"),
-                     call.=FALSE
-                     )
-              retval[,pars]
-            }
+            conv.rec.internal(object=object,pars=pars,transform=transform,...)
           }
           )
 

Modified: pkg/pomp/R/mif.R
===================================================================
--- pkg/pomp/R/mif.R	2013-02-05 22:49:05 UTC (rev 825)
+++ pkg/pomp/R/mif.R	2013-02-26 13:20:27 UTC (rev 826)
[TRUNCATED]

To get the complete diff run:
    svnlook diff /svnroot/pomp -r 826


More information about the pomp-commits mailing list