>>>>> Marc Girondot <marc.giron...@u-psud.fr>
>>>>> on Sat, 10 Sep 2011 11:43:20 +0200 writes:
>>>>> Marc Girondot <marc.giron...@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.giron...@u-psud.fr
> Web: http://www.ese.u-psud.fr/epc/conservation/Marc.html
> Skype: girondot
> ______________________________________________
> R-help@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.
______________________________________________
R-help@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.
