diff --git a/NEWS.md b/NEWS.md
index 7264cfa51..dfc3367a3 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -58,6 +58,7 @@ The function interface remains unchanged.
 
 ## Model changes
 
+- Fixed bounds error in parameters for `dist_fit()` to ensure valid prior and likelihood bounds.
 - MCMC runs are now initialised with parameter values drawn from a distribution that approximates their prior distributions.
 - Added an option to compute growth rates using an estimator by Parag et al. (2022) based on total infectiousness rather than new infections, see `growth_method` argument in rt_opts().
 - Added support for fitting the susceptible population size.
diff --git a/inst/stan/dist_fit.stan b/inst/stan/dist_fit.stan
index 3285ff59a..041ea9234 100644
--- a/inst/stan/dist_fit.stan
+++ b/inst/stan/dist_fit.stan
@@ -1,15 +1,18 @@
 data {
-  int dist; // 0: exp; 1: lnorm; 2: gamma
+  int dist; // 0: exponential; 1: lognormal; 2: gamma
   int N;
   vector[N] low;
   vector[N] up;
-  array[dist == 0] real lam_mean;
-  array[dist > 0] real prior_mean;
-  array[dist > 0] real prior_sd;
+  array[dist == 0] real<lower=0> lam_mean;  // should not be exactly 0
+  array[dist > 0] real<lower=0> prior_mean;
+  array[dist > 0] real<lower=0> prior_sd;
   array[dist == 2] real<lower = 0> par_sigma;
 }
 
 transformed data {
+  if (dist == 0 && lam_mean[1] == 0) {
+    reject("lam_mean must be strictly positive, but found 0");
+  }
   array[dist == 2] real prior_alpha;
   array[dist == 2] real prior_beta;
 
@@ -20,7 +23,7 @@ transformed data {
 }
 
 parameters {
-  array[dist == 0] real<lower = 0> lambda;
+  array[dist == 0] real<lower = 0, upper = 1> lambda_raw;
   array[dist == 1] real mu;
   array[dist == 1] real<lower = 0> sigma;
   array[dist == 2] real<lower = 0> alpha_raw;
@@ -28,42 +31,46 @@ parameters {
 }
 
 transformed parameters{
+  array[dist == 0] real lambda;
   array[dist == 2] real<lower = 0> alpha;
   array[dist == 2] real<lower = 0> beta;
 
+  if (dist == 0) {
+    real lb = 0.2 / lam_mean[1];
+    real ub = 5 / lam_mean[1];
+    lambda[1] = lb + (ub - lb) * lambda_raw[1];
+    // implies: lambda[1] ~ uniform(lb, ub)
+  }
+  
   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];
+    // implies:
+    // alpha[1] ~ normal(prior_alpha[1], par_sigma[1]) T[prior_alpha[1],]
+    // beta[1] ~ normal(prior_beta[1], par_sigma[1]) T[prior_beta[2],]
   }
 }
 
 model {
   if (dist == 0) {
-    lambda[1] ~ uniform(1 / (5. * lam_mean[1]), 1 / (0.2 * lam_mean[1]));
+    for (i in 1:N) {
+      target += log_diff_exp(exponential_lcdf(up[i] | lambda),
+			     exponential_lcdf(low[i] | lambda));
+    }
   } else if (dist == 1) {
     mu[1] ~ normal(prior_mean[1], 10);
-    sigma[1] ~ normal(prior_sd[1], 10) T[0,];
+    // T[0, ] omitted in following: it's constant with data args
+    sigma[1] ~ normal(prior_sd[1], 10);
+    for (i in 1:N) {
+      target += log_diff_exp(lognormal_lcdf(up[i] | mu, sigma),
+			     lognormal_lcdf(low[i] | mu, sigma));
+    }
   } 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_diff_exp(
-        exponential_lcdf(up[i] | lambda),
-        exponential_lcdf(low[i] | lambda)
-      );
-    } else if (dist == 1) {
-      target += log_diff_exp(
-        lognormal_lcdf(up[i] | mu, sigma),
-        lognormal_lcdf(low[i] | mu, sigma)
-      );
-    } else if (dist == 2) {
-      target += log_diff_exp(
-        gamma_lcdf(up[i] | alpha, beta),
-        gamma_lcdf(low[i] | alpha, beta)
-      );
+    alpha_raw[1] ~ std_normal();
+    beta_raw[1] ~ std_normal();
+    for (i in 1:N) {
+      target += log_diff_exp(gamma_lcdf(up[i] | alpha, beta),
+			     gamma_lcdf(low[i] | alpha, beta));
     }
   }
 }