Conversation
|
Important Review skippedAuto incremental reviews are disabled on this repository. Please check the settings in the CodeRabbit UI or the You can disable this status message by setting the Use the checkbox below for a quick retry:
WalkthroughFixed bounds constraints on data variables in the Stan model file. Refactored parameter declarations by introducing a raw parameter with [0,1] bounds and transforming it to the required range. Restructured model branches to use per-item likelihood contributions via log-difference of cumulative distribution functions. Changes
🚥 Pre-merge checks | ✅ 5✅ Passed checks (5 passed)
✏️ Tip: You can configure your own custom pre-merge checks in the settings. ✨ Finishing touches🧪 Generate unit tests (beta)
Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out. Comment |
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
There was a problem hiding this comment.
Actionable comments posted: 2
🤖 Fix all issues with AI agents
In `@inst/stan/dist_fit.stan`:
- Around line 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.
- Around line 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.
🧹 Nitpick comments (2)
NEWS.md (1)
61-61: Prefer a user-facing description over an internal file path.The NEWS entry references the internal Stan file
inst/stan/dist_fit.stan, which is an implementation detail. Per the coding guidelines, NEWS items should describe user-facing changes. Consider rephrasing to focus on what users ofdist_fit()will experience, e.g.:Suggested rewording
-- Fixed bounds error in parameters for `inst/stan/dist_fit.stan` and added necessary data bounds. +- Fixed parameter bounds and added necessary data constraints in the `dist_fit()` model, ensuring valid prior and likelihood computations.As per coding guidelines,
NEWS.md: "Describe user-facing changes, never reference issues or PRs."inst/stan/dist_fit.stan (1)
56-63: Clarify the comment about truncation removal.Line 58's comment
// T[0, ] only needed if loc or scale arguments are paramsis technically correct — the normalizing constant for truncation at 0 depends only on data (prior_sd[1]and10), making it a constant that doesn't affect sampling. However, readers could misread "loc or scale" as referring to the lognormal'smu/sigmarather than the prior normal's arguments.A slightly clearer comment might be:
Suggested wording
- // T[0, ] only needed if loc or scale arguments are params + // T[0, ] omitted: truncation normaliser is constant (args are all data)The interval-censored lognormal likelihood via
log_diff_expon lines 60–62 is correctly implemented.
| 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) | ||
| } |
There was a problem hiding this comment.
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]) | ||
| // implies: beta[1] ~ normal(prior_beta[1], par_sigma[1]) | ||
| } |
There was a problem hiding this comment.
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.
|
This is how benchmark results would change (along with a 95% confidence interval in relative change) if e1294f8 is merged into main:
|
|
A test is now failing because the Stan program is being given a negative standard deviation as data, and I have now included a This is really impressive test coverage! |
|
This is how benchmark results would change (along with a 95% confidence interval in relative change) if 01a0ea1 is merged into main:
|
|
[edit: finishing---didn't mean to submit earlier.] Here's the error that's coming up from the unit tests after I enforced a lower bound of 0 for a variable called Error : Exception: model_dist_fit_namespace::model_dist_fit: prior_sd[1] is -0.0196371, but must be greater than or equal to 0.000000 (in 'dist_fit', line 8, column 2 to column 41)This Let me try to summarize the statistical problem with the model as it currently exists: https://github.com/epiforecasts/EpiNow2/blob/main/inst/stan/dist_fit.stan There are data declarations here that do not constrain data {
array[dist > 0] real prior_mean;
array[dist > 0] real prior_sd;
array[dist == 2] real<lower = 0> par_sigma;
...
transformed data {
array[dist == 2] real prior_beta;
if (dist == 2) {
prior_beta[1] = prior_mean[1] / prior_sd[1]^2;
}
...
parameters {
array[dist == 2] real<lower = 0> beta_raw;
...
transformed parameters {
array[dist == 2] real<lower = 0> beta;
if (dist == 2) {
beta[1] = prior_beta[1] + par_sigma[1] * beta_raw[1];
}
...
model {
beta_raw[1] ~ normal(0, 1);
...
gamma_lcdf(up[i] | alpha, beta),
...Working backwards, we need prior_beta[1] + par_sigma[1] * beta_raw[1] > 0
par_sigma[1] * beta_raw[1] > -prior_beta[1]
beta_raw[1] > -prior_beta[1] / par_sigma[1]Because The usual way I see something like this done is to not constrain beta[1] = exp(prior_beta[1] + par_sigma[1] * beta_raw[1]);As CodeRabbit pointed out, it's not really clear what distribution is being implied for |
|
Thanks - moving this to be modelled with a log link makes sense to me. Just noting that we're planning to remove this model in the not-too-distant future anyway (will put a comment on the linked Issue). |
2 participants
Description
This PR closes #1313
Initial submission checklist
usingusing Stan).devtools::test()anddevtools::check()devtools::document()).and checked usingfor Stan).lintr::lint_package()I made up some data and made sure all three model choices ran without errors and that the data fixes really caught the errors. If there are places to add tests for Stan code like this, I can do that.
After the initial Pull Request
Summary by CodeRabbit
Release Notes