Expectation of the inverse of a shifted squared of a normal random variable!

Question

I am trying to calculate $\int_{-\infty}^{\infty} \frac{p(x)}{x^2+\epsilon}dx$ where $p(x)$ is the density function of a normal random variable. Numerical experiments easily show that it is nothing but $1/\epsilon$ [1]. Can somebody guide me how to solve this?

[1] More precisely I calculate $$\frac{1}{N}\sum_{i=1}^{N}\frac{1}{X_i^2+\epsilon},$$ where $X_i\sim \mathcal{N}(0,1)$.

EDIT: Had a bug in my code! As it is pointed out, the correct convergence is to some point less than $1/\epsilon$.

EDIT 2: To the future reader: Help yourself by reading both of the answers, since they both are constructive.

Answer 1

Use Feynman trick by introducing auxiliary variable $\alpha$ such that the integral is simplified if we differentiate with the respect to $\alpha$: $$I(\alpha)=\frac {e^{\epsilon/2}} {\sqrt {2\pi}}\int_{-\infty}^{\infty} \frac {e^{-\frac 1 2(x^2+\epsilon)\alpha}}{x^2+\epsilon}dx$$ Then we are interested in $$I(1)=I(0)+\int_0^1\frac{\partial I}{\partial \alpha}d\alpha$$ Differentiating inside the integral and integrating over $x$ yields

$$\frac{\partial I}{\partial \alpha}=-\frac 1 2 e^{\epsilon/2} e^{-\epsilon\alpha/2}\frac 1{\sqrt \alpha}$$ Making substitution $\alpha=u^2$ yields

$$\int_0^1\frac{\partial I}{\partial \alpha}d\alpha=-e^{\epsilon/2}\int_0^1e^{-\frac 1 2 \epsilon u^2}du=-1+O(\epsilon)$$

and separately

$$I(0)=e^{\epsilon/2}\sqrt{\frac \pi {2\epsilon}}=\sqrt{\frac \pi {2\epsilon}}+O(\epsilon^{1/2})$$ Combining, we get

$$E\Big(\frac 1 {X^2+\epsilon}\Big)=\sqrt{\frac \pi {2\epsilon}}-1+O(\epsilon^{1/2})$$

This contradicts your numerical results for small $\epsilon$ and it's clear why. For small $\epsilon$, the main contribution to the integral comes from small $x$, so we can approximate the density by a constant $p(0)$, yielding the main term of $\sqrt{\frac \pi {2\epsilon}}$. The only way to get $1/\epsilon$ is to consider large $\epsilon$ so that in the region where $p(x)$ is not small, we can approximate the $\frac 1 {x^2+\epsilon}$ by $\frac 1 \epsilon$.

If you want to confirm this numerically for small $\epsilon$, it's better to do importance sampling and generate $Y_i$ from some distribution concentrated on $x^2$ upto order $\epsilon$ and then reweight function values by ratio of densities. Notice that $\frac {\sqrt \epsilon}\pi\frac 1 {x^2+\epsilon}$ is pdf of a standard Cauchy distributed variable $C$ scaled by $\sqrt \epsilon$ which is easy to generate. Hence expectation of interest can be rewritten as

$$\frac \pi {\sqrt \epsilon} E[p(Y)]$$

where $Y\sim \sqrt\epsilon C$ and $p$ is pdf of a $N(0,1)$. Doing computations this way will give us a much smaller variance than estimation [1].

Answer 2

Of course $\int p(x)/(x^2+\epsilon)dx < \int p(x)/\epsilon\,dx = \frac{1}{\epsilon}$. Thus the answer isn't exactly $1/\epsilon$, let's prove that we can bound the order order though.

We note that for $0<\delta < 1/\sqrt{2\pi}$ we have $$ R_\delta:= \{x:|x| \leq \sqrt{-\log(2\pi \delta^2)}\} $$ is the region on which $p(x) \geq \delta$. Thus we have for any $\delta$ in that range $$ \int \frac{p(x)}{x^2+\epsilon}dx \geq \delta \int_{R_\delta} \frac{1}{x^2+\epsilon}\,dx = \cdots = \frac{2\delta}{\sqrt{\epsilon}}\arctan\left( \frac{\sqrt{-\log(2\pi \delta^2)}}{\sqrt{\epsilon}} \right) $$ and taking $\delta = \sqrt{e^{-\epsilon}}/\sqrt{2\pi}$ we find $$ \int \frac{p(x)}{x^2+\epsilon}dx \geq \frac{\sqrt{\pi}}{2\sqrt{2}} \frac{e^{-\epsilon/2}}{\sqrt{\epsilon}} \geq \frac{\sqrt{\pi}}{2\sqrt{2}} \frac{e^{-1/2}}{\sqrt{\epsilon}} $$ for all $\epsilon < 1$. Thus your integral is bounded between $1/\epsilon$ and $C/\epsilon^{1/2}$ for $\epsilon < 1$.

Moreover, $\epsilon/(x^2+\epsilon) \leq 1$ so by dominated convergence $$ \frac{1}{\frac{1}{\epsilon}}\int\frac{ p(x)}{x^2+\epsilon}\,dx = \int p(x) \frac{ \epsilon}{x^2+\epsilon}\,dx \to \int 0\,dx = 0 $$ so that in fact $1/\epsilon$ is the wrong rate.