# Why doesn't the method of moments work when calculating the variance of the inverse gamma distribution?

I'm trying to calculate the variance of the inverse gamma distribution using the method of movements. According to wikipedia the variance should be:

$$\sigma^2 =\frac{\beta^2}{(\alpha-1)^2(\alpha-2)}$$

Where $\alpha$ is the shape and $\beta$ is the scale of the inverse gamma distribution. When trying this out in R it works reasonably well when $\alpha$ is not to close to 2. For example:

library(MCMCpack) # for the rinvgamma function
a <- 10
b 100

# The variance according to the method of movements
b^2/((a-1)^2*(a-2))
## 15.4321

# The variance by generating inverse gamma distributed random numbers and
# calculating the sample variance
var(rinvgamma(n=9999, shape=a, scale=b)) #
## 15.84388


But when $\alpha$ gets close to 2 the method of movements doesn't seem to work anymore. In the following example the sample variance is much smaller than the method of movements variance:

a <- 2.2
b <- 100

# The variance according to the method of movements
b^2/((a-1)^2*(a-2))
## 34722.22

# The variance by generating inverse gamma distributed random numbers and
# calculating the sample variance
var(rinvgamma(n=9999, shape=a, scale=b)) #
##14479.56


Why doesn't the method of movements work? Am I doing something wrong that can be fixed or is there some other way that I can calculate the variance of an inverse gamma distribution?

• With such a huge theoretical variance, you need a quite large sample to estimate it accurately. Even with the normal distribution you will observe discrepancies var(rnorm(10000,100,35000)) (less severe, of course).
– user10525
Nov 6, 2012 at 11:14
• When $\alpha \rightarrow 2$, the variance goes to infinity. You cannot hope to approach "infinity" accurately Nov 6, 2012 at 11:15
• @Procrastinator: sorry, we have posted at the same time! Nov 6, 2012 at 11:15
• @RasmusBååth Well, you are not using the method of moments precisely. The expression b^2/((a-1)^2*(a-2)) is actually the theoretical variance. The method of moments consists of solving the system $\dfrac{\beta}{\alpha-1}=\bar{x}$ and $\dfrac{\beta^2}{(\alpha-1)^2(\alpha-2)}=\dfrac{1}{n}\sum_{j=1}^n(x_j-\bar{x})^2$ in terms of $(\alpha,\beta)$. After this, you could plug the solutions $(\hat{\alpha},\hat{\beta})$ into the expression for the variance $\dfrac{\beta^2}{(\alpha-1)^2(\alpha-2)}$ in order to obtain an estimator of this quantity.
– user10525
Nov 6, 2012 at 12:30
• @Procrastinator. Ah, ok... Thank you so much for clearing that up! Nov 6, 2012 at 12:37

I've made this community wiki so as not to take any credit for @Procrastinator or anyone else's solution, but this question does have an answer.

The first element to consider is that if the true distribution's $$\alpha$$ is approaching 2 from above, the true variance approaches infinity. Therefore, even if the variance is finite, as $$\alpha$$ approaches 2 from above, the "true" variance grows similar to $$\lim_{x\rightarrow0}\frac{1}{x}$$ (hyperbolically?) and the number of simulations needed to capture that extraordinary growth may well be intractable within any feasible time scale (like a human life). The Central Limit Theorem only holds in cases of finite variance, and as the variance gets out of control, its power weakens. For all we know, the third moment is already infinite and the Berry-Esseen theorem no longer holds, and even though the variance is finite, the rate of convergence is completely unknown and almost certainly glacial. The 9,999 observations in the question simply weren't enough to capture the tail.

A second important point is that $$\sigma^2=\frac{\beta^2}{(\alpha−1)^2(\alpha−2)}$$ is the true distribution given know $$\alpha$$ and $$\beta$$. When using method of moments, we are making the inference in the opposite direction. We have known empirical $$\bar{x}$$ and $$\bar{s}^2$$ from which we want to estimate an $$\hat{\alpha}$$ and $$\hat{\beta}$$. Running through the algebra, if I have not made an error, what is actually being solved is: \begin{aligned} \hat{\alpha} &= \frac{\bar{x}^2}{\bar{s}^2}+2\\ &=\frac{1}{CV^2} + 2\\ \hat{\beta} &= \left(\frac{\bar{x}^2}{\bar{s}^2}+1\right)\bar{x}\\ &=\left(\frac{1}{CV^2} + 1\right)\bar{x} \end{aligned}

While the Law of Large Numbers says these estimates will asymptotically approach their true values, once again, the rate of convergence may be on the order of the heat-death of the universe, and this is more likely as the parameters approach their hard limits. For an inverse gamma with a true $$\alpha \leq 2$$ there is no finite number of simulations whose variance would be captured by MoM. Practically, as the function has a $$+2$$ in it, but more importantly, as the true variance is infinite, the sum of squared errors does not converge. Similarly, For an inverse gamma with a true $$\alpha \leq 1$$ there is no finite number of simulations whose mean would be captured by MoM. Practically, note the $$+1$$. More importantly, in this case, as more observations are generated, larger and larger values are seen and the mean diverges.

• In short the variance is only defined for $\alpha>2$ and the mean for $\alpha>1$. One should use alternative measurements like harmonic mean rather than the mean and full width half maximum (FWHM) for a measure of dispersion rather than the variance.
– Carl
Dec 16, 2021 at 16:21
• The meaning of "huge" variance in your second paragraph is unclear. After all, one could take a beautiful textbook example where MoM works nicely (a Normal distribution, for instance) and just rescale all the data to create a "huge" variance according to any quantitative criterion. Something deeper than that must underlie your sense of "huge," but what is it specifically?
– whuber
Dec 16, 2021 at 18:33
• I was not as precise as I should be, the main issue isn't the variance per se, as you correctly point out, it is $\alpha$ approaching 2, which manifests a variance whose exponent can be measured with five or ore digits. I'll think about how to re-write it better. Also, this is community wiki, please feel free to correct and enhance; I greatly value your insight and input, @whuber. Dec 16, 2021 at 19:22
• This is a thorny issue. Since there isn't anything inherently problematic with $\alpha \gt 2$--the variance is still finite--the question must revolve around the estimation procedure. The solution would seem to located in the chance that with $\alpha \approx 2,$ the sampling variability of the variance will be so great as to preclude precise estimation of $\alpha$ without huge amounts of data. The problem doesn't seem to be inherent in MoM. This seems to be what you're trying to articulate in the second half of this answer.
– whuber
Dec 16, 2021 at 19:56
• The difficulties with convergence have nothing to do with the size of the variance, though. With a simple rescaling of the data, you can make any nonzero variance equal to unity. Convergence becomes a problem when the tails are heavy, as you allude to later in your comment. But that is not measured by variance.
– whuber
Dec 17, 2021 at 16:23