gamma(.., ..)^2 -> generalized_gamma(?, ?, ?)

59 views
Skip to first unread message

S

unread,
May 12, 2017, 8:43:27 PM5/12/17
to Stan users mailing list
Hello,

I have shape and rate parameters of a gamma distribution that seem to be related as some form of rate=shape^2

Now I want to observe if rate followes a generalized_gamma distribution. Howeve my dummy attempts for developing such function failed in test

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

Bob Carpenter

unread,
May 13, 2017, 7:08:32 PM5/13/17
to stan-...@googlegroups.com
First thing to do is test your gamma implementations before
trying to fit.

It helps if you say how things failed.

Rather than log(tgamma(x)), you want lgamma(x).

It should look like this:

return log(c)
+ c * a * log(b)
+ (c * a - 1) * log(x)
- (b * x)^c
- lgamma(a);

then either use the function exposer in RStan or test within the
Stan program to make sure you get the right values using some kind
of independent implementation.

- Bob

> On May 12, 2017, at 8:43 PM, S <mangiol...@gmail.com> wrote:
>
> Hello,
>
> I have shape and rate parameters of a gamma distribution that seem to be related as some form of rate=shape^2
>
> Now I want to observe if rate followes a generalized_gamma distribution. Howeve my dummy attempts for developing such function failed in test
>
> I found the generalised gamma function in BUGS
>
>
>
>
>
>
>
> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

Reply all
Reply to author
Forward
0 new messages