Sampling issues in my model: Reparamterization does not seem to work

Modeling
1

I wrote the below Stan model, but upon running it via RStan and CmdStanR, I get many errors/warnings.

data {
  int<lower=1> N; 
  int<lower=1> Hstar;
  real<lower=1,upper=Hstar> mu_hat;
  real<lower=0> sigma_hat;
  int<lower=1> K; // posterior predictive check
}
parameters {
  real<lower=1,upper=Hstar> H;
  real<lower=0,upper=1> R;
  real<lower=N> Nstar;
}
model {
  H ~ normal(mu_hat, sigma_hat);
  R ~ beta(H, Hstar - H);
  Nstar ~ gamma(N, R); // alternatively, Nstar ~ gamma(N * Hstar, H)
}
generated quantities {
  real<lower=1,upper=Hstar> H_rep[K];
  real<lower=0,upper=1> R_rep[K];
  real<lower=N> Nstar_rep[K];
  for (i in 1:K) {
    H_rep[i] = normal_rng(mu_hat, sigma_hat);
    R_rep[i] = beta_rng(H, Hstar - H);
    Nstar_rep[i] = gamma_rng(N, R); // alternatively, Nstar_rep = gamma_rng(N * Hstar, H) 
  }
}

My data is as follows:

N <- 76
Hstar <- 4
mu_hat <- 3.6475
sigma_hat <- 0.4777725
K <- 5

Specifically, in the generated quantities block, I require H_rep and Nstar_rep to each have the same upper and lower bounds as given in the parameters block. However, Stan doesn’t like this.

Running the model via RStan produces errors and posterior sampling fit is not generated when running via CmdStanR. CmdStanR complains that draws are essentially outside allowable bounds.

I’ve tried increasing both the number of iterations and the adapt_delta parameter, but issues still persist. I really do think a reparameterization of the model is needed.

Any ideas on how to fix the problem?

UPDATE

The following is a reparameterized model, but the issue still persists:

data {
  int<lower=1> N; 
  int<lower=1> Hstar;
  real<lower=1,upper=Hstar> mu_hat;
  real<lower=0> sigma_hat;
  int<lower=1> K; // posterior predictive check
}
parameters {
  real H_tilde; // reparameterized
  real<lower=0,upper=1> R;
  real<lower=N> Nstar;
}
transformed parameters {
  real<lower=1,upper=Hstar> H;
  H = mu_hat + sigma_hat * H_tilde; // reparameterized
}
model {
  H_tilde ~ std_normal(); // reparameterized
  R ~ beta(H, Hstar - H);
  Nstar ~ gamma(N, R); // alternatively, Nstar ~ gamma(N * Hstar, H)
}
generated quantities {
  real H_rep[K];
  real<lower=0,upper=1> R_rep[K];
  real Nstar_rep[K];
  for (i in 1:K) {
    H_rep[i] = normal_rng(mu_hat, sigma_hat);
    R_rep[i] = beta_rng(H, Hstar - H);
    Nstar_rep[i] = gamma_rng(N, R); // alternatively, Nstar_rep = gamma_rng(N * Hstar, H) 
  }
}

I am running on:

R version 3.6.3 (2020-02-29)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.6

2

Can you describe the context of the data you’re trying to model? I find it a little odd to see an int data variable is used as a bound on real parameter.

Oh, and while you’re doing the “non-centered” parameterization by hand, see here for how to use the new-ish <offset=...,multiplier=...> syntax plus a data toggle to achieve the non-centered parameterization more easily.

3

Oh! I’m not sure you need (or even are permitted to express) the bounds on the variables declared in the generated quantities section.

4

Interesting… I’ll try not declaring bounds. Thanks.