A follow-up to 'The meaning of an analyt. result concerning the… mean of the square of a reciprocal of a norm. distrib. rand. variable'

Question

This question concerns the same subject matter as this previous question of mine. However, a moderator felt that the questions I posed there are significantly different from the question I am about to pose here. (The comment in which the moderator makes this observation is the next-to-last comment to this answer.)

Background information

This question should be self-contained, so I will repeat the description of the setting.

Let $X$ be normally distributed with a mean $\mu$ and a 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 is true:

$$\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]\hspace{8em}$$ $$\hspace{4em}=\frac{1}{\sigma^{2}}\left(\mu\,\mathbb{E}[Y]-1\right)$$ $$\hspace{15em}=\frac{1}{\sigma^{2}}\left(\frac{\sqrt{2}\, \mu}{\sigma}F\left(\frac{\mu}{\sqrt{2}\,\sigma}\right)-1\right).\hspace{4em} (1)$$

My question, basically, will be this: what is the meaning (or, interpretation) of $\mathbf{\frac{1}{\boldsymbol{\sigma}^{2}}\left(\frac{\sqrt{2}\, \boldsymbol{\mu}}{\boldsymbol{\sigma}}\boldsymbol{F}\left(\frac{\boldsymbol{\mu}}{\sqrt{2}\,\boldsymbol{\sigma}}\right)-1\right)}$?

The left-hand side of Eq. (1) 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 a derivation of this equation, see the end of this answer of mine. For a probabilistic interpretation of the counterterm, see this answer by whuber.

Comparison to random sampling

Numerical experiments show that for sample sizes as big as $10^{8}$, for $\sigma<0.1$ and $\mu>0.5$, the value of $\frac{1}{\sigma^{2}}\left(\frac{\sqrt{2}\, \mu}{\sigma}F\left(\frac{\mu}{\sqrt{2}\,\sigma}\right)-1\right)$ agrees remarkably well with the mean of $1/x^{2}$-values from the random sample.

The following table presents the results for 100 samples of $10^{8}$ points each, drawn from a normal distribution with $\mu=0.9$ and $\sigma=0.10,\,0.15$, and $0.16$. For each sample, we compute the mean of $1/x^{2}$-values. In the table, we list the smallest and the largest mean obtained among the 100 samples. We see that the analytic prediction $\frac{1}{\sigma^{2}}\left(\frac{\sqrt{2}\, \mu}{\sigma}F\left(\frac{\mu}{\sqrt{2}\,\sigma}\right)-1\right)$ begins to fail rather significantly once $\sigma$ reaches about 0.15.

                                                                        $\mathbf{\boldsymbol\mu=0.9}$
                                         $\sigma$           0.10              0.15             0.16
$\rule{25em}{0.5pt}$
$\frac{1}{\sigma^{2}}\left(\frac{\sqrt{2}\, \mu}{\sigma}F\left(\frac{\mu}{\sqrt{2}\,\sigma}\right)-1\right)$        1.283 39         1.356          1.376
               min sample mean         1.283 32         1.355          1.376
              max sample mean         1.283 47         1.364          4.011

If we make $\mu$ closer to zero, failure will happen for a smaller value of $\sigma$. If $\mu$ is reduced to 0.6 (analytic prediction 3.050), with $\sigma=0.10$, the minimum and maximum mean in one run of 100 samples ($10^{8}$ points each) were 3.050 and 3.446, while in another run of 100 samples, they were 3.050 and 5.636.

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 analytic prediction $\frac{1}{\sigma^{2}}\left(\frac{\sqrt{2}\, \mu}{\sigma}F\left(\frac{\mu}{\sqrt{2}\,\sigma}\right)-1\right)$ must fail. Eventually, as we keep increasing the size of the sample and recomputing the mean of the $1/x^{2}$-values, 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.

A corresponding phenomenon in the context of numerical integration

Whatever causes this stability can also fool numerical integrators.

Using standard methods, we can show that the probability density function for $1/Y^{2}$ is

$$p(z)=\frac{1}{z^{3/2}}\frac{1}{2\, \sqrt{2 \pi }\, \sigma}\Big(e^{-\frac{1}{2 \sigma ^2}(\mu +1/\sqrt{z})^{2}}+e^{-\frac{1}{2 \sigma ^2}(\mu -1/\sqrt{z})^{2}}\Big).$$

One can check that it is properly normalized, $\int_{0}^{\infty}p(z)\,dz=1$, and that the mean $\int_{0}^{\infty}z\,p(z)\,dz$ diverges, since $z\,p(z)$ goes to zero only as $1/\sqrt{z}$ for large $z$.

However, for example Mathematica's numerical integrator ${\scriptstyle\texttt{NIntegrate[z*p[z], \{z, 0, Infinity\}]}}$ will, for $\mu=0.9$ and $\sigma=0.1$, be fooled into thinking that $\int_{0}^{\infty}z\,p(z)\,dz$ converges. In fact, it will produce the output that is numerically identical to $\frac{1}{\sigma^{2}}\left(\frac{\sqrt{2}\, \mu}{\sigma}F\left(\frac{\mu}{\sqrt{2}\,\sigma}\right)-1\right)=$ 1.283.

It is not difficult to understand why the integrator might be fooled. The function $f(z)=z\,p(z)$ is monotonic and concave up for all positive $z$. We may ask the question, how big must $L$ get before it becomes evident that $\int_{0}^{L}f(z)\,dz$ will diverge as $L$ grows without bound? One answer is: this happens when $f(z)$ starts to decay slower than $1/z$. In other words, we are looking for $L^{*}$ such that the slope of $f(z)$ at $L^{*}$ become less negative than the slope of $f(L^{*})(L^{*}/z)$. When $\mu=0.9$ and $\sigma=0.1$, we get $L^{*}\approx 5500$. Now note that $f(L^{*})=2.5\times10^{-19}$, and that $\int_{L^{*}}^{L}\frac{f(L^{*})}{\sqrt{z}}\,dz \sim 5\times 10^{-19}\sqrt{L}$. Since $\int_{0}^{L^{*}}f(z)\,dz \sim 1$, a naive integrator would have to extend the upper limit of integration to $L \sim 10^{36}$ to notice the divergence. (To its credit, Mathematica's ${\scriptstyle\texttt{NIntegrate}}$ does detect the divergence once the lower limit of integration is raised to at least 15.)

An analogy from physics

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 question is:

What is it that $\mathbf{\frac{1}{\boldsymbol{\sigma}^{2}}\left(\frac{\sqrt{2}\, \boldsymbol{\mu}}{\boldsymbol{\sigma}}\boldsymbol{F}\left(\frac{\boldsymbol{\mu}}{\sqrt{2}\,\boldsymbol{\sigma}}\right)-1\right)}$ actually computes?

Here is an example, with $\mathbb{E}[Y]$ as computed through $\frac{\sqrt{2}}{\sigma}F\left(\frac{\mu}{\sqrt{2}\,\sigma}\right)$, of what sort of statement would count as an answer:

The quantity $\frac{\sqrt{2}}{\sigma}F\left(\frac{\mu}{\sqrt{2}\,\sigma}\right)$ will be exactly equal to the mean of an infinitely large sample of $1/x$-values (where the $x$-values themselves are drawn from a normal distribution with a mean $\mu$ and standard deviation $\sigma$).

As we know, the corresponding statement for $\frac{1}{\sigma^{2}}\left(\frac{\sqrt{2}\, \mu}{\sigma}F\left(\frac{\mu}{\sqrt{2}\,\sigma}\right)-1\right)$ does not work, because the mean of $1/x^{2}$-values diverges for infinitely large samples.

So what do we say instead?

The quantity $\mathbf{\frac{1}{\boldsymbol{\sigma}^{2}}\left(\frac{\sqrt{2}\, \boldsymbol{\mu}}{\boldsymbol{\sigma}}\boldsymbol{F}\left(\frac{\boldsymbol{\mu}}{\sqrt{2}\,\boldsymbol{\sigma}}\right)-1\right)}$ will be exactly equal to the mean of… what?

(Here the 'what' will preferably involve random samples. Also, preferably, the statement will assume that $\mu$ and $\sigma$ are fixed finite numbers.)

Answer 1

This thread really concerns an answer in search of a question, which makes it especially difficult to formulate a response. I can, however, provide one derivation leading to statistical interpretations of (1).

---

Suppose $X$ is a random variable and you want to analyze something like the expectation of $|X|^{-p}$ for $p \ge 1.$ As we know, this can be problematic when $X$ has some probability near $0,$ in the sense I described at https://stats.stackexchange.com/a/299765/919, because values of $|X|^{-p}$ "blow up" near $0.$

Suppose, though, that (like Normally distributed random variables), $X$ is absolutely continuous in a neighborhood of $0$ with a cumulative distribution function $F.$ Let's explore the singularity at $0$ by "subtracting off infinity."

Specifically, writing $\mathrm{d}F(x) = f(x)\,\mathrm{d}x$ for the density, consider the expectation-like integral

$$H_f(\epsilon;p) = \int_{-\epsilon}^\epsilon |x|^{-p}\left(f(x) - f(0)\right)\,\mathrm{d}x.$$

This probes the behavior of the expectation within an arbitrarily small interval around its singularity.

Ordinarily this integral would have to be evaluated by splitting it into a negative and positive part and taking limits at the singularity $x=0,$ because the factor $|x|^{-p}$ makes it look like the integrand blows up there. However, if it is the case that the integrand is bounded for sufficiently small $\epsilon$ (which will happen when $f(x)$ approaches $f(0)$ sufficiently rapidly), then no actual singularity appears and this integral makes sense and yields a definite value. In this case, we may fold the positive and negative parts together in a single integral to produce

$$H_f(\epsilon;p) = \int_0^\epsilon x^{-2}\left(f(x) - 2f(0) + f(-x)\right)\,\mathrm{d}x.$$

In the question, $p=2.$ Assuming the density $f = F^\prime$ is differentiable at $0,$ one application of L'Hopital's Rule gives

$$\lim_{x\to 0} |x|^{-2} (f(x) - 2f(0) + f(-x)) = \lim_{x\to 0} \frac{f(x)-2f(0)+f(-x)}{x^2} = \lim_{x\to 0} \frac{f^\prime(x)-f^\prime(-x)}{2x}.$$

When $f^\prime$ is differentiable at $0,$ the numerator approaches $0,$ permitting another application of L'Hopital's Rule,

$$\lim_{x\to 0} \frac{f^\prime(x)-f^\prime(-x)}{2x} = \lim_{x\to 0} \frac{2f^{\prime\prime}(x)}{2}.$$

(Equivalently, you may use Taylor's Theorem with Remainder to obtain the same result.)

The existence of this limit is thus tantamount to supposing $f^{\prime\prime}$ is continuous at $0$ and gives, finally,

$$\lim_{x\to 0} |x|^{-2} (f(x) - 2f(0) + f(-x)) = f^{\prime\prime}(0).$$

This shows that the integrand does not blow up near zero: subtracting $f(0)$ from $f(x)$ in the integration has removed the singularity.

Thus, when $f$ is thrice differentiable in a neighborhood of $0,$

$$H_f(\epsilon, 2) = \int_0^\epsilon \left(f^{\prime\prime}(0) + o(x)\right)\,\mathrm{d}x = \epsilon f^{\prime\prime}(0) + o(\epsilon^2)$$

where the error term $o(\epsilon^2)$ is proportional to the third derivative of $f$ somewhere near $0.$

(The generalization to other integral values of $p$ should now be clear.)

Having removed the singularity, let's consider a version of the original (infinite) expectation where in the formula $E[X^{-2}]=\int x^{-2}\mathrm{d}F(x)$ the probability element $\mathrm{d}F(x)$ is replaced by $\mathrm{d}(F(x) - f(0)x)$ $= (f(x)-f(0))\mathrm{d}x:$

$$\begin{aligned} E^{*}[X^{-2}] &= \int_\mathbb{R} x^{-2}(f(x) - f(0))\,\mathrm{d}x \\ &= \int_{\infty}^{-\epsilon} + \int_{-\epsilon}^\epsilon + \int_{\epsilon}^\infty x^{-2}(f(x) - f(0))\,\mathrm{d}x \\ &= \int_{\infty}^{-\epsilon} + \int_{\epsilon}^\infty x^{-2}f(x)\,\mathrm{d}x\ -\ f(0)\left(\int_{\infty}^{-\epsilon} + \int_{\epsilon}^\infty x^{-2}\,\mathrm{d}x\right) + H_f(\epsilon, 2)\\ &= \int_{\infty}^{-\epsilon} + \int_{\epsilon}^\infty x^{-2}f(x)\,\mathrm{d}x - \frac{2f(0)}{\epsilon} + \left(\epsilon f^{\prime\prime}(0) + o(\epsilon^2)\right) \\ &= \int_{\infty}^{-\epsilon} + \int_{\epsilon}^\infty x^{-2}f(x)\,\mathrm{d}x - \frac{2f(0)}{\epsilon} + o(\epsilon). \end{aligned}$$

When taking the limit, by definition the last $o(\epsilon)$ term will vanish. Thus we may ignore it, permitting the first two re-interpretations of the expression in the question:

Formula $(1)$ in the question is the integral of $(f(x) - f(0))/x^2.$ If you like, when $F$ is absolutely continuous everywhere with nonvanishing density (such as any Normal distribution), this can be written as an expectation: $$\int_{\infty}^{-\epsilon} + \int_{\epsilon}^\infty x^{-2}f(x)\,\mathrm{d}x - \frac{2f(0)}{\epsilon} = E\left[\frac{1 - f(0)/f(X)}{X^2}\right] + o(\epsilon) = E^{*}\left[X^{-2}\right].$$

The latter statement follows upon substituting $f(x)-f(0) = (1 - f(0)/f(x)) f(x)$ in the definition of $E^{*}.$

---

As a practical demonstration of this second interpretation, I performed $100$ simulations with $\mu=2,$ $\sigma=1$ of $10000$ trials each. I generated $x_i$ according to a Normal$(\mu,\sigma)$ distribution (with density $f_{\mu,\sigma}$) and averaged the values of $(1-f_{\mu,\sigma}(0)/f_{\mu,\sigma}(x_i))/x_i^2.$ This is a standard Monte-Carlo estimate of $E^{*}\left[X^{-2}\right].$

The figure is a "funnel plot" of those $100$ Monte-Carlo estimates and their standard errors.

On it I have drawn the value of $E^{*}$ in red, obtained via numerical integration, and an approximate value based on formula $(1)$ of the question with $\epsilon = 1/10$ as a dashed black line: those lines coincide and they coincide with the funnel's tip, which is the best of the Monte-Carlo estimates. This example shows that the two re-interpretations of $(1)$ are practical and statistically meaningful.

# For the expectation:
f <- function(x, mu=0, sigma=1) {
  ifelse(x==0, dnorm(0, mu, sigma) * (mu^2 - sigma^2)/sigma^4/2, 
         (dnorm(x, mu, sigma) - dnorm(0, mu, sigma)) / x^2)
}

# Expression (1):
es <- Vectorize(function(eps, mu, sigma) {
  e.1 <- integrate(function(x) dnorm(x,mu,sigma)/x^2, -Inf, -abs(eps),
                   stop.on.error=FALSE)$value
  e.2 <- integrate(function(x) dnorm(x,mu,sigma)/x^2, abs(eps), Inf, 
               stop.on.error=FALSE)$value
  e.1 + e.2 - 2 * dnorm(0, mu, sigma) / abs(eps)
}, "eps")

# The simulation:
n <- 1e4   # Size of each simulation
mu <- 2    # Normal mean
sigma <- 1 # Normal SD
sim <- replicate(1e2, {
  x <- rnorm(n, mu, sigma)
  y <- f(x, mu, sigma) / dnorm(x, mu, sigma)
  e.star <- mean(y)
  se <- sd(y) / sqrt(length(y))
  c(`E*`=e.star, `Standard error`=se)
})

# Numerical integral:
e.star <- integrate(function(x) f(x,mu,sigma), -Inf, Inf, stop.on.error=FALSE)$value

# Summary of results:
plot(t(sim), log="y", main="Monte-Carlo Estimates", sub="Funnel Plot",
     pch=21, bg=gray(.5,alpha=.25), cex=0.75)
abline(v=e.star, col="Red", lwd=2)
abline(v=es(0.1, mu, sigma), col="Black", lty=2, lwd=2)