Solved – Why doesn’t the gamma density plot match the histogram of samples

data visualizationhistogram

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:
enter image description here

Best Answer

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")

enter image description here