estimating poisson_scale

reconhub · Jun 27, 2019 · 83ddf91 · 83ddf91
1 parent 61af77e
commit 83ddf91
Show file tree

Hide file tree

Showing 4 changed files with 35 additions and 2 deletions.
diff --git a/R/bind_moves.R b/R/bind_moves.R
@@ -135,6 +135,17 @@ bind_moves <- function(moves = custom_moves(), config, data,
                                                 )
   }
 
+  ## remove move$sigma if disabled
+  if (!any(config$move_poisson_scale)) {
+    out$poisson_scale <- NULL
+  } else {
+    out$poisson_scale <- bind_to_function(out$poisson_scale,
+                                          data = data,
+                                          config = config,
+                                          custom_ll = likelihoods$potential_colonised,
+                                          custom_prior = priors$poisson_scale
+                                          )
+  }
 
   ## perform binding for new unknown movements
   known_moves <- names(custom_moves())

diff --git a/R/create_config.R b/R/create_config.R
@@ -180,10 +180,12 @@ create_config <- function(..., data = NULL) {
                    move_t_inf = TRUE,
                    move_mu = TRUE, move_kappa = TRUE, move_pi = TRUE,
                    move_eps = TRUE, move_lambda = TRUE, move_sigma = TRUE,
+                   move_poisson_scale = TRUE,
                    move_potential_colonised = TRUE,
                    n_iter = 1e4, sample_every = 50,
                    sd_mu = 0.0001, sd_pi = 0.1,
                    sd_eps = 0.1, sd_lambda = 0.05, sd_sigma = 0.05,
+                   sd_poisson_scale = 0.1,
                    sd_potential_colonised = 1,
                    prop_alpha_move = 1/4,
                    prop_t_inf_move = 0.2,
@@ -423,7 +425,15 @@ create_config <- function(..., data = NULL) {
     stop("move_potential_colonised is NA")
   }
 
+  ## check move_poisson_scale
+  if (!all(is.logical(config$move_poisson_scale))) {
+    stop("move_poisson_scale is not a logical")
+  }
+  if (any(is.na(config$move_poisson_scale))) {
+    stop("move_poisson_scale is NA")
+  }
 
+
   ## check n_iter
   if (!is.numeric(config$n_iter)) {
     stop("n_iter is not a numeric value")
@@ -503,6 +513,17 @@ create_config <- function(..., data = NULL) {
     stop("sd_sigma is infinite or NA")
   }
 
+  ## check sd_poisson_scale
+  if (!is.numeric(config$sd_poisson_scale)) {
+    stop("sd_poisson_scale is not a numeric value")
+  }
+  if (config$sd_poisson_scale < 1e-10) {
+    stop("sd_poisson_scale is close to zero or negative")
+  }
+  if (!is.finite(config$sd_poisson_scale)) {
+    stop("sd_poisson_scale is infinite or NA")
+  }
+
   ## check sd_potential_colonised
   if (!is.numeric(config$sd_potential_colonised)) {
     stop("sd_potential_colonised is not a numeric value")

diff --git a/R/custom_moves.R b/R/custom_moves.R
@@ -90,7 +90,8 @@ custom_moves <- function(...) {
                      swap_cases = cpp_move_swap_cases,
                      t_inf = cpp_move_t_inf,
                      kappa = cpp_move_kappa,
-                     potential_colonised = cpp_move_potential_colonised
+                     potential_colonised = cpp_move_potential_colonised,
+                     poisson_scale = cpp_move_poisson_scale
                      )
 
 

diff --git a/R/outbreaker_move.R b/R/outbreaker_move.R
@@ -27,7 +27,7 @@ outbreaker_move <- function(moves, data, param_current,
   if(config$pb) {
     pb <- utils::txtProgressBar(min = 1, max = config$n_iter, style = 3)
   }
-  
+
   ## RUN MCMC ##
   for (i in seq.int(2, config$n_iter, 1)) {
     ## move parameters / augmented data