[R] dgamma in jags within r

Marc Girondot marc.girondot at u-psud.fr
Sat Sep 10 11:43:20 CEST 2011


I define priors in jags within r using a gamma distribution. I would 
like to control the shape but I have problems. Any help will be usefull.

 From help of dgamma
___________________
The Gamma distribution with parameters shape = a and scale = s has density
f(x)= 1/(s^a Gamma(a)) x^(a-1) e^-(x/s)
and rate=1/scale

 From jags user manual
____________________
dgamma(r, mu) has a density of
μ^r*x^(r−1)*exp(−μx)
--------------------
Gamma(r)

So I conclude that

µ=1/s then µ in jags is the rate=1/s parameter of dgamma
and r in jags is the shape=a parameter of dgamma

Then

dgamma(r, mu) in jags syntax is dgamma(x, shape=r, rate=mu) in r syntax

# lets try:
jagsgamma <- function(x, r, mu) {(mu^r*x^(r-1)*exp(-mu*r))/gamma(r)}
x <- seq(0,40,by=0.1)
# parameters of the gamma
jagsr=0.001
jagsmu=0.001
# plot the density using the formula of jags
plot(x, jagsgamma(x, jagsr, jagsmu), type="l", xlab="x", ylab="Density")
# plot the density using the dgamma of r
par(new=TRUE)
plot(x, dgamma(x, shape=jagsr, rate=jagsmu), type="l", axes=FALSE, 
col="red", xlab="", ylab="")

It seems to work. Both curves are superimposed.

But it is not at all with these parameters:

# parameters of the gamma
jagsr=10
jagsmu=1
# plot the density using the formula of jags
plot(x, jagsgamma(x, jagsr, jagsmu), type="l", xlab="x", ylab="Density")
# plot the density using the dgamma of r
par(new=TRUE)
plot(x, dgamma(x, shape=jagsr, rate=jagsmu), type="l", axes=FALSE, 
col="red", xlab="", ylab="")

Probably something trivial is wrong but I do not see what.

-- 
__________________________________________________________
Marc Girondot, Pr

Laboratoire Ecologie, Systématique et Evolution
Equipe de Conservation des Populations et des Communautés
CNRS, AgroParisTech et Université Paris-Sud 11 , UMR 8079
Bâtiment 362
91405 Orsay Cedex, France

Tel:  33 1 (0)1.69.15.72.30   Fax: 33 1 (0)1.69.15.73.53
e-mail: marc.girondot at u-psud.fr
Web: http://www.ese.u-psud.fr/epc/conservation/Marc.html
Skype: girondot



More information about the R-help mailing list