Gamma density function

Question

I'm trying different ways to plot for my data. I have a vector ld <- seq(0.001,0.4,0.005) , and parameters for the gamma distribution are 25 (shape) and 240 (rate). My first way for plotting was:

pr <- dgamma(ld, 25,240)
plot(ld,pr,type="b")

I also tried:

pr <- ld^{25-1}*exp(-240*ld)
plot(ld,pr,type="b")

The two plots should be the same, however I found the y-axis is different in scale. I'm wondering how to achieve the same plot using the second method? Thank you!

Answer 1

The plots aren't the same because your formula for the density of the gamma distribution isn't correct.

It should be:

pr <- 240^25 * ld^(25-1)*exp(-240 * ld)/factorial(25 - 1)
plot(ld,pr,type="b")

Or, more generally:

my_dgamma <- function(x, alpha, beta)
{
  beta^alpha * x^(alpha - 1) * exp(-beta * x) / factorial(alpha - 1)
}

(Of course, the denominator could be written as gamma(alpha), but that feels a bit like cheating here.)