I found the generalised gamma function in BUGS
these are the attempts
#1
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 a, real b, real c) {
return( (c*a * (log (c) + log(b)) ) + ((c*a - 1) * log (x) ) + (-(b*x)^c ) - log( tgamma(a) ));
}
}
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)
)
#2
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 a, real b, real c) {
return(log(c*b^(c*a) * x^(c*a - 1) * exp(-(b*x)^c) / tgamma(a)));
}
}
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)
)
Thanks