more testing for girf, igirf, iubf

R/girf.R

@@ -240,7 +240,8 @@ momgirf.internal <- function (object,
     pStop_(paste(sQuote(c("rprocess","dmeasure")),collapse=", ")," are needed basic components.")
 
   if (length(params)==0)
-    stop(ep,sQuote("params")," must be specified",call.=FALSE)
+    pStop_(ep,sQuote("params"),
+      " must be specified if not provided in the spatPomp model object")
 
   tol <- as.numeric(tol)
   gnsi <- as.logical(.gnsi)
@@ -286,9 +287,7 @@ momgirf.internal <- function (object,
         times=times[(nt+2):(nt+1+lookahead_steps)],
         params=params,
         gnsi=TRUE),
-      error = function (e) {
-        pStop_(conditionMessage(e))
-      }
+      error = function (e) pStop_("rprocess error in girf : ", conditionMessage(e))
     )
     fcst_samp_var <- xx
     dim(fcst_samp_var) <- c(length(unit_names(object)), lookahead_steps, Np)
@@ -515,14 +514,7 @@ bootgirf.internal <- function (object,
         params=params.matrix,
         times = times[(nt+2):(nt+1+lookahead_steps)],
         ...),
-      error = function (e) {
-        pomp::flow(object,
-          x0=x,
-          t0=times[nt+1],
-          params=params.matrix,
-          times = times[(nt+2):(nt+1+lookahead_steps)],
-          method = 'adams')
-      }
+      error = function (e) pStop_("girf error in flow : ", conditionMessage(e))
     )
     resids <- Xg - Xskel[,rep(1:Np, each=Nguide),,drop=FALSE] # residuals
     rm(Xg, Xskel, x_with_guides)
@@ -543,14 +535,7 @@ bootgirf.internal <- function (object,
             params=params.matrix,
             times = times[(nt + 1 + 1):(nt + 1 + lookahead_steps)],
             ...),
-          error = function (e) {
-            pomp::flow(object,
-              x0=X.start,
-              t0=tt[s+1],
-              params=params.matrix,
-              times = times[(nt + 1 + 1):(nt + 1 + lookahead_steps)],
-              method = 'adams')
-          }
+          error = function (e) pStop_("girf error in flow : ", conditionMessage(e))
         )
         if(length(znames) > 0){
           skel.start <- skel[,,1]
@@ -586,10 +571,7 @@ bootgirf.internal <- function (object,
           log=TRUE,
           .gnsi=gnsi
         )),
-        error = function (e) {
-          stop(ep,"error in calculation of log_dmeas_weights: ",
-            conditionMessage(e),call.=FALSE)
-        }
+        error = function (e) pStop_("girf error in calculation of log_dmeas_weights: ", conditionMessage(e))
         )
         ldw <- array(log_dmeas_weights, c(U,Nguide,Np)) # log_dmeas_weights is an array with dim U*(Np*Nguide)*1. Reorder it as U*Nguide*Np
         log_fcst_lik <- colSums(log(apply(exp(ldw),c(1,3),sum)/Nguide)) # average dmeas (natural scale) over Nguide sims, then take log, and then sum over 1:U (for each particle)

R/igirf.R

@@ -281,9 +281,10 @@ igirf.momgirf <- function (object, params, Ninter, lookahead, Nguide,
   if (length(tol) != 1 || !is.finite(tol) || tol < 0)
     pStop_(ep, sQuote("tol")," should be a small positive number.")
 
-  do_ta <- length(.indices)>0L
-  if (do_ta && length(.indices)!=Np)
-    pStop_(ep, sQuote(".indices")," has improper length.")
+  ## ancestor tracking is not supported
+  ## do_ta <- length(.indices)>0L
+  ## if (do_ta && length(.indices)!=Np)
+  ##   pStop_(ep, sQuote(".indices")," has improper length.")
 
   times <- time(object,t0=TRUE)
   ntimes <- length(times)-1
@@ -332,9 +333,7 @@ igirf.momgirf <- function (object, params, Ninter, lookahead, Nguide,
         times=times[(nt+2):(nt+1+lookahead_steps)],
         params=tp_with_guides,
         gnsi=TRUE),
-      error = function (e) {
-        pStop_("in ", sQuote("igirf"), " : ", conditionMessage(e)) 
-      }
+      error = function (e) pStop_("in ", sQuote("igirf"), " : ", conditionMessage(e)) 
     )
     fcst_samp_var <- xx
     dim(fcst_samp_var) <- c(length(unit_names(object)), lookahead_steps, Np)
@@ -372,9 +371,7 @@ igirf.momgirf <- function (object, params, Ninter, lookahead, Nguide,
           times=times[(nt+2):(nt+1+lookahead_steps)],
           params=tparams,
           gnsi=TRUE),
-        error = function (e) {
-          pStop_("in ", sQuote("igirf"), " : ",conditionMessage(e)) 
-        }
+        error = function (e) pStop_("in ", sQuote("igirf"), " : ",conditionMessage(e)) 
       )
       dim(meas_var_skel) <- c(length(unit_names(object)), lookahead_steps, Np)
       fcst_var_upd <- array(0, dim = c(length(unit_names(object)), lookahead_steps, Np))
@@ -522,9 +519,10 @@ igirf.bootgirf <- function (object, params, Ninter, lookahead, Nguide,
   if (length(tol) != 1 || !is.finite(tol) || tol < 0)
     stop(ep,sQuote("tol")," should be a small positive number.",call.=FALSE)
 
-  do_ta <- length(.indices)>0L
-  if (do_ta && length(.indices)!=Np)
-    stop(ep,sQuote(".indices")," has improper length.",call.=FALSE)
+  ## ancestor tracking is not supported
+  ## do_ta <- length(.indices)>0L
+  ## if (do_ta && length(.indices)!=Np)
+  ##   pStop_(ep,sQuote(".indices")," has improper length.")
 
   times <- time(object,t0=TRUE)
   ntimes <- length(times)-1
@@ -597,49 +595,45 @@ igirf.bootgirf <- function (object, params, Ninter, lookahead, Nguide,
       }
       log_guide_fun <- vector(mode = "numeric", length = Np) + 1
 
-      for(l in 1:lookahead_steps){
-        if(nt+1+l-lookahead <= 0) discount_denom_init <- object@t0
-        else discount_denom_init <- times[nt+1+l - lookahead]
-        discount_factor <- 1 - (times[nt+1+l] - tt[s+1])/(times[nt+1+l] - discount_denom_init)/ifelse(lookahead==1,2,1)
+      for(lag in 1:lookahead_steps){
+        if(nt+1+lag-lookahead <= 0) discount_denom_init <- object@t0
+        else discount_denom_init <- times[nt+1+lag - lookahead]
+        discount_factor <- 1 - (times[nt+1+lag] - tt[s+1])/(times[nt+1+lag] - discount_denom_init)/ifelse(lookahead==1,2,1)
 	## to ensure that the discount factor does not become too small for L=1
 	## and small s (which can lead to very uninformative guide function), increase
 	## the discount factor to at least 1/2 when L=1.
 
         ## construct pseudo-simulations by adding simulated noise terms (residuals)
 	## to the skeletons
-        pseudosims <- skel[,rep(1:Np, each=Nguide),l,drop=FALSE] +
+        pseudosims <- skel[,rep(1:Np, each=Nguide),lag,drop=FALSE] +
           resids[,rep(guidesim_index-1, each=Nguide)*
-	    Nguide+rep(1:Nguide, Np),l,drop=FALSE] -
+	    Nguide+rep(1:Nguide, Np),lag,drop=FALSE] -
           resids[,rep(guidesim_index-1, each=Nguide)*
 	    Nguide+rep(1:Nguide, Np),1,drop=FALSE] +
           resids[,rep(guidesim_index-1, each=Nguide)*
 	    Nguide+rep(1:Nguide, Np),1,drop=FALSE] *
 	      sqrt((times[nt+2]-tt[s+1])/(times[nt+2]-times[nt+1]))
-        rm(skel)
         log_dmeas_weights <- tryCatch(
         (vec_dmeasure(
           object,
-          y=object@data[,nt+l,drop=FALSE],
+          y=object@data[,nt+lag,drop=FALSE],
           x=pseudosims,
-          times=times[nt+1+l],
+          times=times[nt+1+lag],
           params=tp_with_guides,
           log=TRUE,
           .gnsi=gnsi
         )),
-        error = function (e) {
-          stop(ep,"error in calculation of log_dmeas_weights: ",
-            conditionMessage(e),call.=FALSE)
-        }
+        error = function (e) pStop_("error in calculation of log_dmeas_weights: ", conditionMessage(e))
         )
 
         ## uncomment for debugging: tracking down particles with NA weight
         ## if(any(is.na(log_dmeas_weights))){
         ##   na_ix <- which(is.na(log_dmeas_weights[,,1]))[1]
         ##   na_ix_col <- (na_ix %/% nunits) + (na_ix %% nunits > 0)
         ##   illegal_dunit_measure_error(
-        ##     time=times[nt+1+l],
+        ##     time=times[nt+1+lag],
         ##     lik=log_dmeas_weights[,na_ix_col,1,drop=FALSE],
-        ##     datvals=object@data[,nt+l],
+        ##     datvals=object@data[,nt+lag],
         ##     states=pseudosims[,na_ix_col,1L],
         ##     params=tp_with_guides[,na_ix_col]
         ##   )
@@ -655,6 +649,7 @@ igirf.bootgirf <- function (object, params, Ninter, lookahead, Nguide,
         log_resamp_weights <- log_fcst_lik*discount_factor
         log_guide_fun = log_guide_fun + log_resamp_weights
       }
+      rm(skel)
       log_s_not_1_weights <- log_guide_fun - log_filter_guide_fun
       if (!(s==1 & nt!=0)){
         log_weights <- log_s_not_1_weights
@@ -673,10 +668,7 @@ igirf.bootgirf <- function (object, params, Ninter, lookahead, Nguide,
             log=TRUE,
             .gnsi=gnsi
           ),
-          error = function (e) {
-            stop(ep,"error in calculation of log_meas_weights: ",
-              conditionMessage(e),call.=FALSE)
-          }
+          error = function (e) pStop_("error in calculation of log_meas_weights: ", conditionMessage(e))
         )
         gnsi <- FALSE
         log_weights <- as.numeric(log_meas_weights) + log_s_not_1_weights
@@ -707,9 +699,7 @@ igirf.bootgirf <- function (object, params, Ninter, lookahead, Nguide,
             fsv=array(0,dim=c(length(unit_names(object)), lookahead_steps, Np)), 
             tol=tol
           ),
-          error = function (e) {
-            stop(ep,conditionMessage(e),call.=FALSE) 
-          }
+          error = function (e) pStop_("error in igirf computations : ", conditionMessage(e)) 
         )
         guidesim_index = guidesim_index[xx$ancestry] 
         cond.loglik[nt+1, s] <- xx$loglik + max_log_weights

R/iubf.R

@@ -104,6 +104,8 @@ setMethod(
       pStop_(ep,sQuote("nbhd")," must be specified")
     if (missing(Nparam))
       pStop_(ep,sQuote("Nparam")," must be specified")
+    if (Nparam<2.5)
+      pStop_(ep,sQuote("Nparam")," must be at least 3")
     if (missing(prop))
       pStop_(ep,sQuote("prop")," must be specified")
     cooling.type <- match.arg(cooling.type)

tests/bm.R

@@ -2,7 +2,7 @@ library(spatPomp)
 i <- 1
 U <- switch(i,2,10)
 N <- switch(i,2,10)
-Np <- switch(i,10,100)
+Np <- switch(i,5,100)
 
 # For CRAN tests, need to limit to two cores
 # https://stackoverflow.com/questions/50571325/r-cran-check-fail-when-using-parallel-functions
@@ -45,6 +45,8 @@ b_model_no_runit_measure <- spatPomp(b_model,runit_measure=NULL)
 b_model_no_dunit_measure <- spatPomp(b_model,dunit_measure=NULL)
 b_model_t0_equal_t1 <- spatPomp(b_model,t0=1)
 b_model5 <- bm(U=U,N=5) 
+b_model_with_accumvars <- b_model
+b_model_with_accumvars@accumvars <- rownames(states(b_model))
 
 b_model_zero_dmeasure <- spatPomp(b_model,
   dmeasure = spatPomp_Csnippet(
@@ -212,14 +214,14 @@ paste("check girf on girfd_spatPomp: ",
 
 ## check girf for zero measurement density situations
 b_girf_mom_inf <- girf(b_model_zero_dmeasure,Np = floor(Np/2),lookahead = 1,
-  Nguide = floor(Np/2),
+  Nguide = 3,
   kind = 'moment',Ninter=2)
 paste("bm moment girf loglik, zero measurement: ",
   round(logLik(b_girf_mom_inf),10))
 
 set.seed(0)
 b_girf_boot_inf <- girf(b_model_zero_dmeasure,Np = floor(Np/2),lookahead = 1,
-  Nguide = floor(Np/2),
+  Nguide = 3,
   kind = 'bootstrap',Ninter=2)
 paste("bm bootstrap girf loglik, zero measurement: ",
   round(logLik(b_girf_boot_inf),10))
@@ -240,6 +242,22 @@ try(girf(b_model,kind='boot'))
 try(girf(b_model,kind='boot',Np=5))
 try(girf(b_model,kind='boot',Np=5,Nguide=3,tol=1:1000))
 
+try(girf(b_model_no_params,Np = 3,lookahead = 1, Nguide = 3,
+  kind = 'moment',Ninter=2))
+try(girf(b_model_no_params,Np = 3,lookahead = 1, Nguide = 3,
+  kind = 'boot',Ninter=2))
+
+try(girf(b_model,Np = 1:10,lookahead = 1, Nguide = 3,
+  kind = 'moment',Ninter=2))
+try(girf(b_model,Np = "JUNK",lookahead = 1, Nguide = 3,
+  kind = 'moment',Ninter=2))
+
+girf(b_model_with_accumvars,Np = 3,lookahead = 2, Nguide = 3,
+  kind = 'moment',Ninter=2)
+girf(b_model_with_accumvars,Np = 3,lookahead = 2, Nguide = 3,
+  kind = 'boot',Ninter=2)
+
+
 ## ------------------------------------------------------------
 ## Now, we test the inference methods
 ## ____________________________________________________________
@@ -264,12 +282,12 @@ b_igirf_geom <- igirf(b_model,
   cooling.fraction.50 = 0.5,
   Np=Np,
   Ninter = 2,
-  lookahead = 1,
-  Nguide = 5,
+  lookahead = 2,
+  Nguide = 3,
   kind = 'moment',
   verbose = FALSE
 )
-paste("bm igirf loglik, geometric cooling, verbose=F: ",round(logLik(b_igirf_geom),10))
+paste("bm igirf loglik, geometric cooling, verbose=F: ",round(logLik(b_igirf_geom),5))
 
 set.seed(1)
 b_igirf_geom_repeat <- igirf(b_igirf_geom,params=coef(b_model))
@@ -281,9 +299,9 @@ b_igirf_hyp <- igirf(b_model,
   rw.sd = b_rw.sd,
   cooling.type = "hyperbolic",
   cooling.fraction.50 = 0.5,
-  Np=floor(Np/2),
+  Np=3,
   Ninter = 2,
-  lookahead = 1,
+  lookahead = 2,
   Nguide = floor(Np/2),
   kind = 'moment',
   verbose = TRUE
@@ -313,10 +331,10 @@ b_igirf_boot_hyp <- igirf(b_model,
   rw.sd = b_rw.sd,
   cooling.type = "hyperbolic",
   cooling.fraction.50 = 0.5,
-  Np=floor(Np/2),
+  Np=3,
   Ninter = 2,
-  lookahead = 1,
-  Nguide = floor(Np/2),
+  lookahead = 2,
+  Nguide = 3,
   kind = 'bootstrap',
   verbose = TRUE
 )
@@ -349,6 +367,13 @@ try(igirf(b_model,Ngirf=2,rw.sd=b_rw.sd,cooling.fraction.50=0.5,
   cooling.type="geometric",Np=4))
 try(igirf(b_model,Ngirf=2,rw.sd=b_rw.sd,cooling.fraction.50=0.5,
   cooling.type="geometric",Np="JUNK",Nguide=4))
+try(igirf(b_model,Ngirf="JUNK",rw.sd=b_rw.sd,cooling.fraction.50=0.5,
+  cooling.type="geometric",Np=3,Nguide=4))
+try(igirf(b_model,Ngirf=2,rw.sd=b_rw.sd,cooling.fraction.50=1000,
+  cooling.type="geometric",Np=3,Nguide=4))
+
+try(igirf(b_model,kind="moment",Ngirf=2,rw.sd=b_rw.sd,cooling.fraction.50=0.5,cooling.type="geometric",Np=3,Nguide=4,tol=-1))
+
 
 igirf(b_igirf_boot_geom,Np=3,
   .paramMatrix=cbind(coef(b_model),coef(b_model),coef(b_model)))
@@ -357,20 +382,19 @@ try(igirf(b_igirf_boot_geom,Np=function(x) "JUNK"))
 try(igirf(b_igirf_boot_geom,Np=5:15))
 
 
-try(igirf(b_model,kind='moment', Ngirf = 2, Nguide=2,
+igirf(b_model_with_accumvars,kind='moment', Ngirf = 2, Nguide=2,
   rw.sd = b_rw.sd, cooling.type = "hyperbolic", cooling.fraction.50 = 0.5,
-  Np=floor(Np/2), Ninter = 2))
-
+  Np=3, Ninter = 1)
+igirf(b_model_with_accumvars,kind='boot', Ngirf = 2, Nguide=2,
+  rw.sd = b_rw.sd, cooling.type = "hyperbolic", cooling.fraction.50 = 0.5,
+  Np=3, Ninter = 1)
 
-# try(igirf(b_model_no_vunit_measure,kind='moment'))
-# try(igirf(b_model_no_rprocess,kind='moment'))
-# try(igirf(b_model,kind='moment'))
-# try(igirf(b_model,kind='moment',Np=5))
-# try(igirf(b_model,kind='moment',Np=5,Nguide=3,tol=1:1000))
-# try(igirf(b_model_no_rprocess,kind='boot'))
-# try(igirf(b_model,kind='boot'))
-# try(igirf(b_model,kind='boot',Np=5))
-# try(igirf(b_model,kind='boot',Np=5,Nguide=3,tol=1:1000))
+igirf(b_model_zero_dmeasure,kind='moment', Ngirf = 2, Nguide=2,
+  rw.sd = b_rw.sd, cooling.type = "hyperbolic", cooling.fraction.50 = 0.5,
+  Np=3, Ninter = 1)
+igirf(b_model_zero_dmeasure,kind='boot', Ngirf = 2, Nguide=2,
+  rw.sd = b_rw.sd, cooling.type = "hyperbolic", cooling.fraction.50 = 0.5,
+  Np=3, Ninter = 1)
 
 
 ## ----------------------------------------------------------
@@ -436,8 +460,8 @@ paste("bm iubf loglik, geometric cooling, verbose=F: ",round(logLik(b_iubf_geom)
 
 b_iubf_hyp <- iubf(b_model,
   Nubf = 2,
-  Nrep_per_param = floor(Np/2),
-  Nparam = floor(Np/2),
+  Nrep_per_param = 3,
+  Nparam = 3,
   nbhd = b_bag_nbhd,
   prop = 0.8,
   rw.sd =b_rw.sd,
@@ -454,13 +478,14 @@ try(iubf(b_model,Nubf=2,Nrep_per_param=3))
 try(iubf(b_model,Nubf=2,Nrep_per_param=3,rw.sd=b_rw.sd))
 try(iubf(b_model,Nubf=2,Nrep_per_param=3,rw.sd=b_rw.sd,cooling.fraction.50=1000))
 try(iubf(b_model,Nubf=2,Nrep_per_param=3,rw.sd=b_rw.sd,cooling.fraction.50=0.5))
-try(iubf(b_model,Nubf=2,Nrep_per_param=3,rw.sd=b_rw.sd,cooling.fraction.50=0.5,Nparam=2))
-try(iubf(b_model,Nubf=2,Nrep_per_param=3,rw.sd=b_rw.sd,cooling.fraction.50=0.5,Nparam=2,nbhd=b_bag_nbhd))
-try(iubf(b_model,Nubf=2,Nrep_per_param=1,rw.sd=b_rw.sd,cooling.fraction.50=0.5,Nparam=2,nbhd=b_bag_nbhd))
+try(iubf(b_model,Nubf=2,Nrep_per_param=3,rw.sd=b_rw.sd,cooling.fraction.50=0.5,Nparam=3))
+try(iubf(b_model,Nubf=2,Nrep_per_param=3,rw.sd=b_rw.sd,cooling.fraction.50=0.5,Nparam=3,nbhd=b_bag_nbhd))
+try(iubf(b_model,Nubf=2,Nrep_per_param=1,rw.sd=b_rw.sd,cooling.fraction.50=0.5,Nparam=3,nbhd=b_bag_nbhd))
 
 ## max_lookback is not triggered for b_model with N=2
 iubf(b_model5,Nubf=2, Nparam = 3,Nrep_per_param=3,nbhd=b_bag_nbhd,rw.sd=b_rw.sd,cooling.fraction.50=0.5,prop=0.8) 
 
+try(iubf(b_model5,Nubf=2, Nparam = 2,Nrep_per_param=3,nbhd=b_bag_nbhd,rw.sd=b_rw.sd,cooling.fraction.50=0.5,prop=0.8))
 
 ## --------------------------------------------
 ## using bm to test simulate and plot