[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