I’m continuing the model from Using integration to fit the difference of gamma distributed random variable here, because this is not about integration anymore. I’m trying to follow @nhuurre’s suggestion to fit the difference between random variables with a gamma distribution, and I was trying to do SBC on that. But I’m having problems to just fit the model:

This is the model:

```
functions {
real gammadiffs_lpdf(real y, real alpha1, real alpha2, real beta, real g2) {
real g1 = g2 + y; // Y = g1 - g2
return gamma_lpdf(g1| alpha1, beta) + gamma_lpdf(g2| alpha2, beta);
}
}
data {
int<lower = 0> J;
}
transformed data {
real<lower = 0> alpha_1_ = lognormal_rng(log(4), .1);
real<lower = 0> alpha_2_ = lognormal_rng(log(.5), .1);
real<lower = 0> beta_ = lognormal_rng(log(.5), .1);
real y[J];
for (j in 1:J){
real X1 = 0;
real X2 = 0;
X1 = gamma_rng(alpha_1_, beta_);
X2 = gamma_rng(alpha_2_, beta_);
y[j] = X1 - X2;
}
}
parameters {
real<lower=0> alpha_1;
real<lower=0> alpha_2;
real<lower=0> beta;
vector<lower=0>[J] g2;
}
transformed parameters {
}
model {
target += lognormal_lpdf(alpha_1 | log(4), .1);
target += lognormal_lpdf(alpha_2 | log(.5), .1);
target += lognormal_lpdf(beta | log(.5), .1);
for(j in 1:J)
target += gammadiffs_lpdf(y[j]| alpha_1, alpha_2, beta, g2[j]);
}
generated quantities {
real y_[J] = y;
vector[3] pars_;
int ranks_[3] = {alpha_1 > alpha_1_,
alpha_2 > alpha_2_,
beta > beta_};
vector[J] log_lik;
pars_[1] = alpha_1_;
pars_[2] = alpha_2_;
pars_[3] = beta_;
for(j in 1:J)
log_lik[j] = gammadiffs_lpdf(y[j]| alpha_1, alpha_2, beta, g2[j]);
}
```

I’m having problems with the initial values, 3 out of 10 fits they get rejected with

```
Chain 1: Rejecting initial value:
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can't start sampling from this initial value.
```

But sometimes they do work.

What might be going on here? I think I didn’t forget any `lower`

, did I?

This is not even the model that I want to fit, this was the very simple one that was supposed to work :/