# Why doesn't my gamma density plot match my histogram of samples?

I'm confused as to why plot is showing such a steep curve for the gamma PDF. I thought that plotting the draws from an inverse-gamma should approximate the PDF, however it doesn't match the plotted PDF given the same parameters. I take the inverse of the rgamma(), which makes more sense if you can see the variables I'm using. However, I've replaced variables with hard-coded values to keep my MWE brief. I am using R and my code is below:

grid <- seq(0,100,by=0.1) # hard-coded for MWE
sig.post.shape <- 91 # hard-coded for MWE
sig.post.rate <- 1247.52 # hard-coded for MWE
set.seed(1)
hist(1/rgamma(grid, shape = sig.post.shape, rate = sig.post.rate), breaks=10)
plot(grid,1/dgamma(grid, shape = sig.post.shape, rate = sig.post.rate), type = "l")


The resulting graph is:

• In all forums explaining what software you use is courteous and informative and will be until everyone uses the same software. Even people using the same software (presumably R) would like to see your graphs to start thinking about the question. People using different software can't be expected to write their own code to do this. – Nick Cox Mar 25 '15 at 9:27
• @NickCox Fixed. – user2205916 Mar 25 '15 at 10:25

If $X$ has density (pdf) $f$ then $1/X$ does not have density $1/f$ (which is not even a density). Indeed, the change-of-variables formula teaches us that $Y:=h(X)$ has density $$\bigl|{(h^{-1})}'(y)\bigr|f(h^{-1}(y))$$ when $h$ is a "nice" invertible transformation.

I gave the formula for a general $h$, because $h^{-1}(y)=1/y$ when $h(x)=1/x$, and that could cause some confusion.

Then denoting by $f$ the density of your simulated Gamma distribution (rgamma), you have to compare your histogram with the density $$\frac{1}{y^2} f\left(\frac{1}{y}\right)$$ :

grid <- seq(0,100,by=0.1)
sig.post.shape <- 91
sig.post.rate <- 1247.52
set.seed(1);
hist(1/rgamma(grid, shape = sig.post.shape, rate = sig.post.rate), breaks=10, prob=TRUE)
lines(grid,1/grid^2*dgamma(1/grid, shape = sig.post.shape, rate = sig.post.rate), type = "l")


The problem here isn't at root anything to do with graphics.

The density of an inverse gamma distribution is not the reciprocal of the density of a gamma distribution, which is what your last syntax line implies. To see that, it's sufficient to note that such a relation would map near zero densities of the gamma to near infinite densities of the inverse gamma. In effect, this is what R is trying to draw, with the bizarre results you note. The function you draw isn't even a density function.

UPDATE. Stéphane Laurent followed this quickly, with a much fuller, definitive version. I am letting this stand as a Mickey Mouse "executive summary" answer.