Merge 5b02799 into dee8c4f

DESCRIPTION

@@ -2,7 +2,7 @@ Type: Package
 Package: EpiNow2
 Title: Estimate Real-Time Case Counts and Time-Varying
     Epidemiological Parameters
-Version: 1.3.6.4000
+Version: 1.3.6.5000
 Authors@R: 
     c(person(given = "Sam",
              family = "Abbott",

NEWS.md

@@ -14,6 +14,7 @@ This release is in development. For a stable release install 1.3.5 from CRAN.
   code base with `data.table::fcase()`. By @jamesmbaazam in #383 and reviewed by @seabbs.
 * Reviewed the example in `calc_backcalc_data()` to call `calc_backcalc_data()`
   instead of `create_gp_data()`. By @jamesmbaazam in #388 and reviewed by @seabbs.
+* Improved compilation times by reducing the number of distinct stan models and deprecated `tune_inv_gamma()`. By @sbfnk in #394 and reviewed by @seabbs.
 
 # EpiNow2 1.3.5
 

R/dist.R

@@ -216,26 +216,26 @@ dist_fit <- function(values = NULL, samples = NULL, cores = 1,
     N = length(values),
     low = lows,
     up = ups,
-    iter = samples + 1000,
-    warmup = 1000
+    lam_mean = numeric(0),
+    prior_mean = numeric(0),
+    prior_sd = numeric(0),
+    par_sigma = numeric(0)
   )
 
+  model <- stanmodels$dist_fit
+
   if (dist %in% "exp") {
-    model <- stanmodels$exp
-    data <- c(data, lam_mean = mean(values))
-  } else if (dist %in% "gamma") {
-    model <- stanmodels$gamma
-    data <- c(data,
-      prior_mean = mean(values),
-      prior_sd = sd(values),
-      par_sigma = 1.0
-    )
+    data$dist <- 0
+    data$lam_mean <- array(mean(values))
   } else if (dist %in% "lognormal") {
-    model <- stanmodels$lnorm
-    data <- c(data,
-      prior_mean = log(mean(values)),
-      prior_sd = log(sd(values))
-    )
+    data$dist <- 1
+    data$prior_mean <- array(log(mean(values)))
+    data$prior_sd <- array(log(sd(values)))
+  } else if (dist %in% "gamma") {
+    data$dist <- 2
+    data$prior_mean <- array(mean(values))
+    data$prior_sd <- array(sd(values))
+    data$par_sigma <- array(1.0)
   }
 
   # set adapt delta based on the sample size
@@ -249,6 +249,8 @@ dist_fit <- function(values = NULL, samples = NULL, cores = 1,
   fit <- rstan::sampling(
     model,
     data = data,
+    iter = samples + 1000,
+    warmup = 1000,
     control = list(adapt_delta = adapt_delta),
     chains = chains,
     cores = cores,
@@ -484,9 +486,15 @@ bootstrapped_dist_fit <- function(values, dist = "lognormal",
 
 
     out <- list()
-    out$mean_samples <- sample(rstan::extract(fit)$mu, samples)
-    out$sd_samples <- sample(rstan::extract(fit)$sigma, samples)
-
+    if (dist == "lognormal") {
+      out$mean_samples <- sample(rstan::extract(fit)$mu, samples)
+      out$sd_samples <- sample(rstan::extract(fit)$sigma, samples)
+    } else if (dist == "gamma") {
+      alpha_samples <- sample(rstan::extract(fit)$alpha, samples)
+      beta_samples <- sample(rstan::extract(fit)$beta, samples)
+      out$mean_samples <- alpha_samples / beta_samples
+      out$sd_samples <- sqrt(alpha_samples) / beta_samples
+    }
     return(out)
   }
 
@@ -741,12 +749,12 @@ sample_approx_dist <- function(cases = NULL,
 
 #' Tune an Inverse Gamma to Achieve the Target Truncation
 #'
-#' @description `r lifecycle::badge("questioning")`
+#' @description `r lifecycle::badge("deprecated")`
 #' Allows an inverse gamma distribution to be. tuned so that less than 0.01 of
 #' its probability mass function falls outside of the specified bounds. This is
 #' required when using an inverse gamma prior, for example for a Gaussian
 #' process. As no inverse gamma priors are currently in use and this function
-#' has some stability issues it may be deprecated at a later date.
+#' has some stability issues it has been deprecated.
 #'
 #' @param lower Numeric, defaults to 2. Lower truncation bound.
 #'
@@ -756,32 +764,20 @@ sample_approx_dist <- function(cases = NULL,
 #' distribution that achieves the target truncation.
 #' @export
 #'
-#' @examples
+#' @keywords internal
 #'
-#' tune_inv_gamma(lower = 2, upper = 21)
 tune_inv_gamma <- function(lower = 2, upper = 21) {
-  model <- stanmodels$tune_inv_gamma
-  # optimise for correct upper and lower probabilites
-  fit <- rstan::sampling(model,
-    data = list(
-      u = upper,
-      l = lower
-    ),
-    iter = 1,
-    warmup = 0,
-    chains = 1,
-    algorithm = "Fixed_param",
-    refresh = 0
-  )
-
-  alpha <- rstan::extract(fit, "alpha")
-  beta <- rstan::extract(fit, "beta")
-
-  out <- list(
-    alpha = round(unlist(unname(alpha)), 1),
-    beta = round(unlist(unname(beta)), 1)
+  lifecycle::deprecate_stop(
+    "1.3.6", "tune_inv_gamma()",
+    details = paste0(
+      "As no inverse gamma priors are currently in use and this function has ",
+      "some stability issues it has been deprecated. Please let the package ",
+      "authors know if deprecating this function has caused any issues. ",
+      "For the last active version of the function see the one contained ",
+      "in version 1.3.5 at ",
+      "https://github.com/epiforecasts/EpiNow2/blob/bad836ebd650ace73ad1ead887fd0eae98c52dd6/R/dist.R#L739" # nolint
+    )
   )
-  return(out)
 }
 
 #' Specify a distribution.

R/stanmodels.R

@@ -1,18 +1,15 @@
 # Generated by rstantools.  Do not edit by hand.
 
 # names of stan models
-stanmodels <- c("estimate_infections", "estimate_secondary", "estimate_truncation", "exp", "gamma", "lnorm", "simulate_infections", "simulate_secondary", "tune_inv_gamma")
+stanmodels <- c("dist_fit", "estimate_infections", "estimate_secondary", "estimate_truncation", "simulate_infections", "simulate_secondary")
 
 # load each stan module
+Rcpp::loadModule("stan_fit4dist_fit_mod", what = TRUE)
 Rcpp::loadModule("stan_fit4estimate_infections_mod", what = TRUE)
 Rcpp::loadModule("stan_fit4estimate_secondary_mod", what = TRUE)
 Rcpp::loadModule("stan_fit4estimate_truncation_mod", what = TRUE)
-Rcpp::loadModule("stan_fit4exp_mod", what = TRUE)
-Rcpp::loadModule("stan_fit4gamma_mod", what = TRUE)
-Rcpp::loadModule("stan_fit4lnorm_mod", what = TRUE)
 Rcpp::loadModule("stan_fit4simulate_infections_mod", what = TRUE)
 Rcpp::loadModule("stan_fit4simulate_secondary_mod", what = TRUE)
-Rcpp::loadModule("stan_fit4tune_inv_gamma_mod", what = TRUE)
 
 # instantiate each stanmodel object
 stanmodels <- sapply(stanmodels, function(model_name) {

_pkgdown.yml

@@ -114,7 +114,6 @@ reference:
     desc: Functions to define, fit and parameterise distributions
     contents:
       - contains("dist")
-      - tune_inv_gamma
   - title: Simulate
     desc: Functions to help with simulating data or mapping to reported cases
     contents:

inst/stan/dist_fit.stan

@@ -0,0 +1,69 @@
+data {
+  int dist; // 0: exp; 1: lnorm; 2: gamma
+  int N;
+  vector[N] low;
+  vector[N] up;
+  real lam_mean[dist == 0];
+  real prior_mean[dist > 0];
+  real prior_sd[dist > 0];
+  real<lower = 0> par_sigma[dist == 2];
+}
+
+transformed data {
+  real prior_alpha[dist == 2];
+  real prior_beta[dist == 2];
+
+  if (dist == 2) {
+    prior_alpha[1] = (prior_mean[1] / prior_sd[1])^2;
+    prior_beta[1] = prior_mean[1] / prior_sd[1]^2;
+  }
+}
+
+parameters {
+  real<lower = 0> lambda[dist == 0];
+  real mu[dist == 1];
+  real<lower = 0> sigma[dist == 1];
+  real<lower = 0> alpha_raw[dist == 2];
+  real<lower = 0> beta_raw[dist == 2];
+}
+
+transformed parameters{
+  real<lower = 0> alpha[dist == 2];
+  real<lower = 0> beta[dist == 2];
+
+  if (dist == 2) {
+    alpha[1] = prior_alpha[1] + par_sigma[1] * alpha_raw[1];
+    beta[1] = prior_beta[1] + par_sigma[1] * beta_raw[1];
+  }
+}
+
+model {
+  if (dist == 0) {
+    lambda[1] ~ uniform(1 / (5. * lam_mean[1]), 1 / (0.2 * lam_mean[1]));
+  } else if (dist == 1) {
+    mu[1] ~ normal(prior_mean[1], 10);
+    sigma[1] ~ normal(prior_sd[1], 10) T[0,];
+  } else if (dist == 2) {
+    alpha_raw[1] ~ normal(0, 1);
+    beta_raw[1] ~ normal(0, 1);
+  }
+
+  for(i in 1:N){
+    if (dist == 0) {
+      target += log(
+        exponential_cdf(up[i] , lambda) -
+        exponential_cdf(low[i], lambda)
+      );
+    } else if (dist == 1) {
+      target += log(
+        lognormal_cdf(up[i], mu, sigma) -
+        lognormal_cdf(low[i], mu, sigma)
+      );
+    } else if (dist == 2) {
+      target += log(
+        gamma_cdf(up[i], alpha, beta) -
+        gamma_cdf(low[i], alpha, beta)
+      );
+    }
+  }
+}