Skip to content

Issue1313 param bounds #1314

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
bob-carpenter wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Issue1313 param bounds #1314

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
6fab81f
corrected bounds; fixes #1313
bob-carpenter Feb 11, 2026
e1294f8
news for issue #1313
bob-carpenter Feb 11, 2026
01a0ea1
code review doc patches
bob-carpenter Feb 11, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions NEWS.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
61 changes: 34 additions & 27 deletions inst/stan/dist_fit.stan
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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;

Expand All @@ -20,50 +23,54 @@ 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;
array[dist == 2] real<lower = 0> beta_raw;
}

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)
}
Comment on lines +38 to +43
Copy link
Copy Markdown
Contributor

@coderabbitai coderabbitai Bot Feb 11, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟡 Minor

Division by zero if lam_mean[1] is exactly zero.

The lower=0 bound on lam_mean (line 6) permits a value of exactly 0, which would cause division by zero here when computing lb and ub. Consider using lower=0 with a runtime check or tightening to a small positive lower bound (e.g. real<lower=1e-10>).

This is likely only a theoretical concern since passing a zero mean for an exponential distribution is degenerate, but it leaves an unguarded code path.

🤖 Prompt for AI Agents 
In `@inst/stan/dist_fit.stan` around lines 35 - 40, The code computes lb and ub by
dividing by lam_mean[1], which can be zero due to lam_mean's lower=0 bound and
causes division by zero; update the model to guard against zero by either
enforcing a strictly positive lower bound on lam_mean (e.g., change its
declaration to real<lower=1e-10> lam_mean[...] ) or add a runtime check before
computing lb/ub (e.g., when dist == 0 verify lam_mean[1] > 0 and handle or
throw) and ensure any use of lambda_raw/lambda (symbols lambda_raw, lambda,
lam_mean, and the dist branch) uses the safe nonzero value.


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],]
}
Comment on lines 45 to 51
Copy link
Copy Markdown
Contributor

@coderabbitai coderabbitai Bot Feb 11, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟡 Minor

Comments claim normal but effective prior is shifted half-normal.

alpha_raw and beta_raw are declared with lower=0 (lines 26–27), and the sampling statements use std_normal() (lines 65–66). This gives a half-normal(0, 1) prior on the raw parameters, not a standard normal. After the transform alpha = prior_alpha + par_sigma * alpha_raw, the implied prior on alpha is a shifted, scaled half-normal with support [prior_alpha, ∞), not a symmetric normal(prior_alpha, par_sigma) as the comments on lines 45–46 state.

If the asymmetric (half-normal) prior is intentional, the comments should be corrected:

Suggested fix 
-    // implies: alpha[1] ~ normal(prior_alpha[1], par_sigma[1])
-    // implies: beta[1] ~ normal(prior_beta[1], par_sigma[1])
+    // implies: alpha[1] ~ prior_alpha[1] + par_sigma[1] * half_normal(0, 1)
+    // implies: beta[1] ~ prior_beta[1] + par_sigma[1] * half_normal(0, 1)

If a symmetric normal was intended instead, the lower=0 constraints on alpha_raw/beta_raw should be removed (and positivity of alpha/beta ensured via rejection or a different parameterisation).

Also applies to: 64-66

🤖 Prompt for AI Agents 
In `@inst/stan/dist_fit.stan` around lines 42 - 47, The comments claiming alpha
and beta "imply: ... ~ normal(prior_*, par_sigma)" are wrong because alpha_raw
and beta_raw are declared lower=0 and drawn with std_normal(), producing a
half-normal(0,1) on the raw scale; after transform alpha = prior_alpha +
par_sigma * alpha_raw (and similarly for beta) the implied prior is a shifted,
scaled half-normal with support [prior_alpha, ∞), not a symmetric normal. Fix by
either (A) updating the comments next to alpha, beta (and the analogous block
using std_normal()) to state the correct prior: "implies: alpha[1] ~
shifted_scaled_half_normal(prior_alpha[1], par_sigma[1]) with support
[prior_alpha[1], ∞)" (and same for beta), or (B) if a symmetric normal was
intended, remove the lower=0 constraints on alpha_raw and beta_raw so
std_normal() yields a full normal on the raw scale (ensuring any positivity
constraints on alpha/beta are handled differently); adjust the corresponding
sampling statements and comments accordingly.

}

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));
}
}
}
Loading