[R] dgamma in jags within r

[Martin Maechler]

>>>>> Marc Girondot <marc.girondot at u-psud.fr>
>>>>>     on Sat, 10 Sep 2011 11:43:20 +0200 writes:
>>>>> Marc Girondot <marc.girondot at u-psud.fr>
>>>>>     on Sat, 10 Sep 2011 11:43:20 +0200 writes:

    > 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

all correct..

    > # lets try:
    > jagsgamma <- function(x, r, mu) {(mu^r*x^(r-1)*exp(-mu*r))/gamma(r)}

well: Here (above) is the error....
I hope you were always just asking "where did I go wrong?", 
as you just can be sure that R is right ...
Hint: It's one letter inside 'exp(*)'..
      

As I did not see your typo immediately,
I've improved the following and "donate" the code here:

p.both.gamma <- function(x, r.jags, mu.jags, ylab = "Density", ...) {
    ## plot the density using the formula of jags
    matplot(x, cbind(jagsgamma(x, r.jags, mu.jags),
                     dgamma(x, shape=r.jags, rate=mu.jags)),
            type="l", lty=1, ylab=ylab, ...)
    mtext(substitute(list(r[jags] == R, mu[jags] == M),
                     as.list(formatC(c(R=r.jags, M=mu.jags)))))
    legend("topright", c("jagsgamma", "dgamma"), lty=1, col=1:2, bty = "n")
}

x <- seq(0,40, by=0.1)
# parameters of the gamma
p.both.gamma(x, r.jags = 0.001, mu.jags = 0.001)
p.both.gamma(x, r.jags = 0.001, mu.jags = 0.001,
             log = "xy")
## It seems to work. Both curves are superimposed.

## MM: something in between:
p.both.gamma(x, r.jags = 0.1, mu.jags = 0.5)
p.both.gamma(x, r.jags = 0.1, mu.jags = 0.5, log = "xy")

## But it is not at all with these parameters:

p.both.gamma(x, r.jags = 10, mu.jags = 1)




    > 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

    > ______________________________________________
    > R-help at r-project.org mailing list
    > https://stat.ethz.ch/mailman/listinfo/r-help
    > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
    > and provide commented, minimal, self-contained, reproducible code.