Thanks I’m currently testing S1 with discrete success (S2 never completed with the following code)
set.seed(1) # for reproducible example
s = 50
r = 3
x <- rgamma(10000,shape=s,rate=r)
(stan(
model_code=“
functions{
real ggamma_lpdf(real x, real k, real mu, real sigma) {
real y;
real w;
real d;
y = log(x);
w = (y-mu)/sigma;
d = (k-.5)*log(k) - log(sigma) - lgamma(k) +
(sqrt(k)w - kexp(1/sqrt(k)*w)) - y;
return d;
}
}
data {
int<lower=1> G;
vector[G] y;
}
parameters { real<lower=0> h[3]; }
model { for(g in 1:G) y[g] ~ ggamma(h[1], h[2], h[3]); }
”,
data = list(G=length(x), y=x)
))
Here the more general picture of what I’m trying to achieve. Maybe there is a fundamental issue with the way I am thinking of my hyerarchical negative binomial process --> Hyerarchical negative binomial process