Search code examples
rnumericgamma-distribution

R calculated gamma distribution density wrong?

I am trying to numerically calculate marginal likelihood (marginalize over a positive parameter). I am using Gamma distribution as prior for that parameter. Here I looked at the behavior of Gamma distribution for two specific parameter settings:

s = 28.4; r = 17000
plot(x, dgamma(x, shape=s, rate = r), type = 'l', ylab = 'density')
abline(v = s/r, col = 'red')

I got the following results:

density plot1

Then I tried the following to get a tighter Gamma distribution: 

lines(x, dgamma(x, shape=s*1000, rate = r*1000), col = 'blue')

and the result:

density plot2

I am confused. As distribution gets tighter, the height should grow taller, otherwise the area won't integrate to 1. Did I miss anything? Or is there any numerally problems? Thanks!

Solution

  • Your x variable needs to have more samples to capture the narrow peak in the second density function:

    x = seq(0, .01, .000001)
s = 28.4; r = 17000
plot(x, dgamma(x, shape=s, rate = r), type = 'l', ylab = 'density')
abline(v = s/r, col = 'red')
lines(x, dgamma(x, shape=s*1000, rate = r*1000), col = 'blue')

    enter image description here