>>>>> 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.