diff --git a/DESCRIPTION b/DESCRIPTION
index a4e077603..9779ecdde 100644
--- a/DESCRIPTION
+++ b/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",
diff --git a/NEWS.md b/NEWS.md
index 70b3a6105..e90d8d73f 100644
--- a/NEWS.md
+++ b/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
 
diff --git a/R/dist.R b/R/dist.R
index 35e5dd0b2..dd249df27 100644
--- a/R/dist.R
+++ b/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.
diff --git a/R/stanmodels.R b/R/stanmodels.R
index 120f51678..8c6dc0118 100644
--- a/R/stanmodels.R
+++ b/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) {
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 0520d7330..e83c873fb 100644
--- a/_pkgdown.yml
+++ b/_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:
diff --git a/inst/stan/dist_fit.stan b/inst/stan/dist_fit.stan
new file mode 100644
index 000000000..b247d14f1
--- /dev/null
+++ b/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)
+      );
+    }
+  }
+}
diff --git a/inst/stan/exp.stan b/inst/stan/exp.stan
deleted file mode 100644
index a067a39f6..000000000
--- a/inst/stan/exp.stan
+++ /dev/null
@@ -1,18 +0,0 @@
-data {
-  int N;
-  vector[N] low;
-  vector[N] up;
-  real lam_mean;
-}
-parameters {
-  real<lower = 0> lambda;
-}
-model {
-  lambda ~ uniform(1/(5*lam_mean),1/(0.2*lam_mean));
-
-  for(i in 1:N){
-    target += log(exponential_cdf(up[i] , lambda) - exponential_cdf(low[i] , lambda));
-  }
-
-}
-
diff --git a/inst/stan/gamma.stan b/inst/stan/gamma.stan
deleted file mode 100644
index 2ae54fa12..000000000
--- a/inst/stan/gamma.stan
+++ /dev/null
@@ -1,40 +0,0 @@
-
-data {
-  int N;
-  vector[N] low;
-  vector[N] up;
-  real prior_mean;
-  real prior_sd;
-  real <lower = 0> par_sigma;
-}
-
-transformed data {
-  real prior_alpha = ((prior_mean) / prior_sd)^2;
-  real prior_beta = (prior_mean) / (prior_sd^2);
-}
-
-parameters {
-  real <lower = 0> alpha_raw;
-  real <lower = 0> beta_raw;
-}
-
-transformed parameters{
-  real <lower = 0> alpha;
-  real <lower = 0> beta;
-
-  alpha = prior_alpha + par_sigma * alpha_raw;
-  beta = prior_beta + par_sigma * beta_raw;
-}
-model {
-  alpha_raw ~ normal(0, 1);
-  beta_raw ~ normal(0, 1);
-
-  for(i in 1:N){
-    target += log(gamma_cdf(up[i] , alpha, beta) - gamma_cdf(low[i] , alpha, beta));
-  }
-}
-
-generated quantities {
-  real mu = alpha / beta;
-  real sigma = sqrt(alpha / (beta^2));
-}
diff --git a/inst/stan/lnorm.stan b/inst/stan/lnorm.stan
deleted file mode 100644
index 24de645a9..000000000
--- a/inst/stan/lnorm.stan
+++ /dev/null
@@ -1,23 +0,0 @@
-
-data {
-  int N;
-  vector[N] low;
-  vector[N] up;
-  real prior_mean;
-  real prior_sd;
-}
-
-parameters {
-  real mu;
-  real<lower = 0> sigma;
-}
-
-model {
-  mu ~ normal(prior_mean, 10);
-  sigma ~ normal(prior_sd, 10) T[0,];
-
-  for(i in 1:N){
-   target += log(lognormal_cdf(up[i] , mu, sigma) - lognormal_cdf(low[i] , mu, sigma));
-  }
-}
-
diff --git a/inst/stan/tune_inv_gamma.stan b/inst/stan/tune_inv_gamma.stan
deleted file mode 100644
index 899fb4881..000000000
--- a/inst/stan/tune_inv_gamma.stan
+++ /dev/null
@@ -1,41 +0,0 @@
-// Adapted from here:
-// https://betanalpha.github.io/assets/case_studies/gaussian_processes.html
-
-functions {
-  vector tail_delta(vector y, vector theta, real[] x_r, int[] x_i) {
-    vector[2] deltas;
-    real alpha = exp(y[1]);
-    real beta = exp(y[2]);
-    alpha = alpha <= 0 ? 1e-6 : alpha;
-    beta = beta <= 0 ? 1e-6 : beta;
-    alpha = alpha >= 1e6 ? 1e6 : alpha;
-    beta = beta >= 1e6 ? 1e6 : beta;
-    deltas[1] = inv_gamma_cdf(theta[1], alpha, beta) - 0.01;
-    deltas[2] = 1 - inv_gamma_cdf(theta[2], alpha, beta) - 0.01;
-    return deltas;
-  }
-}
-
-data {
- real u;
- real l;
-}
-
-transformed data {
-  vector[2] theta = [l, u]';
-
-  real delta = 1;
-  real a = square(delta * (u + l) / (u - l)) + 2;
-  real b =  ((u + l) / 2) * (square(delta * (u + l) / (u - l)) + 1);
-  vector[2] y_guess = [log(a), log(b)]';
-
-  real x_r[0];
-  int x_i[0];
-
-  vector[2] y = algebra_solver(tail_delta, y_guess, theta, x_r, x_i);
-}
-
-generated quantities {
-  real alpha = exp(y[1]);
-  real beta = exp(y[2]);
-}
diff --git a/man/EpiNow2-package.Rd b/man/EpiNow2-package.Rd
index 3f75b2a6c..cb1ba1935 100644
--- a/man/EpiNow2-package.Rd
+++ b/man/EpiNow2-package.Rd
@@ -30,7 +30,7 @@ Authors:
   \item Hamada S. Badr \email{badr@jhu.edu} (\href{https://orcid.org/0000-0002-9808-2344}{ORCID})
   \item Michael DeWitt \email{me.dewitt.jr@gmail.com} (\href{https://orcid.org/0000-0001-8940-1967}{ORCID})
   \item EpiForecasts
-  \item Sebastian Funk \email{sebastian.funk@lshtm.ac.uk}
+  \item Sebastian Funk \email{sebastian.funk@lshtm.ac.uk} (\href{https://orcid.org/0000-0002-2842-3406}{ORCID})
 }
 
 Other contributors:
diff --git a/man/tune_inv_gamma.Rd b/man/tune_inv_gamma.Rd
index 1db2c38f3..c1c907d99 100644
--- a/man/tune_inv_gamma.Rd
+++ b/man/tune_inv_gamma.Rd
@@ -16,14 +16,11 @@ A list of alpha and beta values that describe a inverse gamma
 distribution that achieves the target truncation.
 }
 \description{
-\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#questioning}{\figure{lifecycle-questioning.svg}{options: alt='[Questioning]'}}}{\strong{[Questioning]}}
+\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#deprecated}{\figure{lifecycle-deprecated.svg}{options: alt='[Deprecated]'}}}{\strong{[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.
-}
-\examples{
-
-tune_inv_gamma(lower = 2, upper = 21)
+has some stability issues it has been deprecated.
 }
+\keyword{internal}