Gamma density plot: dgamma working while own function returning error

Question

I am trying to plot a Gamma density function for the parameters (shape) alfa=3457 and (rate) beta=84. If I do this using:

curve(dgamma(x,shape=3457,rate=84),from=35,to=50,xlab="posterior theta",ylab="density")

everything works out just fine and I get the density centred on ~41. While if I first construct the density as:

para_density = function(x,s=3457,r=84,N) ((r^s)/(gamma(s)))*x^(s-1)*exp(-r*x)

and then plot it with:

curve(para_density(x,s=3457,r=84,N=N),from=0,to=100)

R tells me gamma(s) returns Inf, which is sensible given its meaning.

How come then I can plot it with identical parametrisation using the plot and dgamma functions?

Thank y'all in advance.

Answer 1

Well, just use logarithm then

Code (untested!)

mygamma <- function(x, s, r) {
    l <- s*log(r) - lgamma(s) + (s-1.0)*log(x) - x*r
    exp(l)
}

UPDATE

Yes, it works

library(ggplot2)

p <- ggplot(data = data.frame(x = 0), mapping = aes(x = x))
p <- p + stat_function(fun = function(x) dgamma(x,shape=3457,rate=84))
p <- p + stat_function(fun = function(x) mygamma(x,s=3457,r=84))
p <- p + xlim(35.0, 50.0)
print(p)