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