The meaning of an analytical result concerning the (formally nonexistent) mean of the square of a reciprocal of a normally distributed random variable

Question

This question arose as I was writing this answer to this question.

Let $X$ be normally distributed with mean $\mu$ and standard deviation $\sigma$, and let $Y=1/X$.

First, note that the integral representing $\mathbb{E}[Y]$ (the mean value of $Y$) exists in the principal value sense. It may be related to the Hilbert transform of a Gaussian, which, in turn, is expressible in terms of the Dawson function $F(x)=e^{-x^{2}}\int_{0}^{x}e^{t^{2}}dt$. One obtains $\mathbb{E}[Y]=\frac{\sqrt{2}}{\sigma}F\left(\frac{\mu}{\sqrt{2}\,\sigma}\right)$. For references, see e.g. here, here, and this paper.

At the same time, the integral representing the mean value of $Y^{2}$ does not exist, not even in the principal value sense.

However, the following seems to be true (for a derivation, see the end of the answer of mine I already mentioned above):

$$\lim_{\epsilon\to 0^{+}} \left[\int_{-\infty}^{-\epsilon}\frac{1}{x^{2}}\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}}dx+\int_{\epsilon}^{\infty}\frac{1}{x^{2}}\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}}dx-\frac{1}{\epsilon}\frac{2}{\sqrt{2\pi}\sigma}e^{-\frac{\mu^{2}}{2\sigma^{2}}}\right]=\frac{1}{\sigma^{2}}\left(\mu\,\mathbb{E}[Y]-1\right).\hspace{4em} (1)$$

This is almost the principal value of $\mathbb{E}[Y^{2}]$, except that it also includes the 'counterterm' $\frac{1}{\epsilon}\frac{2}{\sqrt{2\pi}\sigma}e^{-\frac{\mu^{2}}{2\sigma^{2}}}$, for which I, at least, have no interpretation in terms of probability theory.

Comparison to random sampling

Nevertheless, numerical experiments show that for sample sizes as big as $10^{8}$, for $\sigma=1/10$ and $\mu>0.5$, the 'prediction' $\frac{1}{\sigma^{2}}\left(\mu\,\mathbb{E}[Y]-1\right)$ agrees remarkably well with the mean of $1/x^{2}$-values from the random sample. In one typical run, for $\mu=0.9$, I got 1.283 39 from the formula vs 1.283 41 from the random sample. In 100 runs with the same parameters, the minimum mean was 1.283 32 and the maximum 1.283 47. If $\mu$ is reduced to 0.6 (analytic prediction 3.050), the minimum and maximum mean were 3.050 and 3.446.

To illustrate the dependence on $\sigma$, let's keep $\mu$ at 0.9. For $\sigma=0.15$, the analytic prediction is 1.356, while the min and max mean (again, 100 samples of $10^{8}$ points each) were 1.355 and 1.364. If $\sigma=0.16$, the analytic prediction is 1.376, while the min and max mean were 1.376 and 4.011; we see that the prediction is starting to fail here.

I realize that eventually, for any fixed values of $\mu$ and $\sigma$, if we keep increasing the size of the sample (or repeat sufficiently often the comparison with samples of the same size) the prediction must fail. Eventually, as we keep increasing the size of the sample and recomputing the mean of the squares of the inverses, we will find that this mean keeps increasing with the size of the sample.

But there also seems to be a remarkable stability for samples of 'intermediate' sizes.

This situation reminds me a little bit of the 'conflict' between the Second Law of Thermodynamics and the Poincaré recurrence theorem: yes, the theorem says that if we wait long enough, the entropy of an isolated system must come as close as we wish to its initial value. But this happens on enormous time scales; on any 'practical' time scale, the entropy of the system cannot decrease. In fact, as far as anything observable and macroscopic in our universe, only the result for the 'intermediate' time scales matters.

My questions are:

1. Is there a probability-theory interpretation of Eq. (1), and if yes, what is it?

1.(a) In particular, is there a probability-theory interpretation of the right-hand side of Eq. (1)? In terms of probability theory, what is it actually computing? What is its probabilistic interpretation?

2. Why does the right-hand side of Eq. (1) agree so well with the results from the random samples of large-but-not-truly-enormously-large size?

3. Suppose we will be taking samples of size $N$. Suppose that $\sigma$ is fixed. How small can $\mu$ get before we have to start to worry that the mean of $1/x^{2}$-values will start to become unpredictable? Or, equivalently, for a given $\mu$, how large can $N$ get before we have to start to worry?

Answer 1

Because (2) and (3) have easy, clear answers, let's focus on (1).

It is revealing to eliminate the distracting details of this situation. Let us, then, consider any distribution that is absolutely continuous in a neighborhood $\mathcal U$ of $0$ with a continuous density function $f$ in that neighborhood. (To simplify the notation in the following, let's extend $f$ if necessary to a function of $\mathbb R$ by setting it to zero at all points not in $\mathcal U.$) Because $f$ is nonnegative, the problem comes down to analyzing

$$\int_{\mathcal U} \frac{1}{x^2} f(x)\,\mathrm{d}x = \lim_{\epsilon\to 0^-}\int_{-\infty}^{-\epsilon} \frac{f(x)}{x^2} \, \mathrm{d}x + \lim_{\epsilon\to 0^-}\int_\epsilon^\infty \frac{ f(x)}{x^2}\, \mathrm{d}x = \lim_{\epsilon\to 0^-}\int_\epsilon^\infty \frac{g(x)}{x^2} \, \mathrm{d}x\tag{*}$$

where $g(x) = f(x) + f(-x).$ (The limits define the left hand side, at least for Riemann integrals.)

When $f(0)\ne 0$ the limits diverge (as I showed in detail at https://stats.stackexchange.com/a/299765/919). But at what rate? We can study this by defining a family of (Pareto) random variables $X_\epsilon$ supported on $[\epsilon,\infty)$ for $\epsilon\gt 0$ with cumulative distribution function

$$F(x;\epsilon) = \Pr(X_\epsilon \le x) = 1 - \frac{\epsilon}{x},\ x \ge \epsilon.$$

Because $f$ is integrable, $g$ is a (non-negative) measurable function so $g(X_\epsilon)$ has an expectation (possibly infinite) given by

$$E[g(X_\epsilon)] = \int_{\mathbb{R}} g(x)\,\mathrm{d}F(x;\epsilon) = \int_\epsilon^\infty g(x)\,\frac{\epsilon}{x^2}\,\mathrm{d}x = \epsilon\int_\epsilon^\infty \frac{g(x)}{x^2}\,\mathrm{d}x.$$

Now as $\epsilon\to 0$ (from above), the probability distribution of $X_\epsilon$ becomes concentrated near $0.$ (Proof: pick any quantile $\alpha$ and observe that the $\alpha$ quantile of $X_\epsilon,$ equal to $\epsilon/(1-\alpha),$ converges to zero.) Since we have assumed $g$ is continuous at $0,$ this implies the expectation of $g(X_\epsilon)$ converges to $g(0).$ Subtracting and taking limits yields

$$\lim_{\epsilon\to0^+} \epsilon \int_\epsilon^\infty \frac{g(x)}{x^2}\,\mathrm{d}x - g(0) = 0.$$

This implies the difference between the integral in $(*)$ and $g(0)/\epsilon$ is bounded for sufficiently small $\epsilon.$ In this sense, the integral missing from the principal part formula, $$\int_{-\epsilon}^\epsilon \frac{f(x)}{x^2}\,\mathrm{d}x,$$ is diverging like $2f(0)/\epsilon.$

---

When $f$ is the density of a Normal$(\mu,\sigma)$ distribution,

$$g(0) = 2f(0) = \frac{2}{\sigma\sqrt{2\pi}}\exp\left(-\frac{\mu^2}{2\sigma^2}\right).$$

Because this is $1/(\sigma\sqrt{2\pi})$ times the value cited in the question, I performed a quick numerical check.

The curve on the left (for a standard Normal distribution) stays bounded as $\epsilon$ approaches zero, as claimed here. The curve on the right, using the result in the question, diverges.