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:

 A: 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. 
A: 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")


