Expectation of the reciprocal of the Gaussian tail

Question

Let $N(\mu, \sigma^2)$ be the Gaussian distribution. For a real number $l\in\mathbb{R}$, we further denote by $N(l; \mu,\sigma^2)$ the Gaussian tail distribution, namely, $N(l;\mu,\sigma^2)$ has the density $$ p(x) = \begin{cases} 0 & \text{if } x \leq l,\\ C e^{-\frac{(x-\mu)^2}{2\sigma^2}} & \text{otherwise}. \end{cases} $$ Here, $C$ is a positive number which makes the integral equal to $1$. Note that if $l>0$ and $X\sim N(l; \mu,\sigma^2)$, then $\frac{1}{X}$ is well-defined.

Suppose $X\sim N(l; \mu,\sigma^2)$ is distributed with respect to the Gaussian tail distribution where $l>0$. What is $\mathbb{E}[\frac{1}{X}]$?

I calculated $C$ to be $$ C = \frac{\sqrt{2}}{\mathrm{erfc}(\frac{l-\mu}{\sqrt{2}\sigma})\sqrt{\pi}\sigma}, $$ where $\mathrm{erfc}$ denotes the complementary error function. With this, I have $$ \mathbb{E}[\frac{1}{X}] = C \int_{l}^{\infty} \frac{1}{x} e^{-\frac{(x-\mu)^2}{2\sigma^2}}\mathrm{d} x. $$ In the case that $\mu=0$, this integral can be expressed in terms of the incomplete $\Gamma$ function. If I am not mistaken, the integral equals $\Gamma(0;\frac{l^2}{2})/2$ in this case.

However, I was not able to solve the integral for $\mu\neq 0$.

Edit: I want to use the closest integer to $\mathbb{E}[\frac{1}{X}]$ in a code. Hence, I am also happy with a good estimate or a fast algorithm that computes it (preferably, the algorithm does not actually compute the integral).

Edit2 : I would also like to add the following diagrams, that explain the behaviour of $\mathbb{E}[\frac{1}{X}]$ in several regimes.

In each picture, $\mu$ and $\sigma$ is given in the title, the $x$-axis represents different values of $l$ and the $y$-axis is the value of $\mathbb{E}[\frac{1}{X}]$. The best one is the last one: Here we have $\mu=100$ and $\sigma=2$. In this case, for small values of $l$ we always have $\mathbb{E}[\frac{1}{X}]\approx 0.01$, which makes sense.

Answer 1

We want the following integral $$ I = \int_{\ell}^\infty f(x) \;dx $$ with $$ f(x)=\frac{e^{-\frac{(x-\mu)^2}{2\sigma^2}}}{x} $$ Use the following change of variables $$ x = \sigma\tan\left(\frac{\pi z}{2}\right)+\ell $$ $$ dx = \frac{\pi\sigma}{2}\sec^2\left(\frac{\pi z}{2}\right) \;dz $$ $$ x=\ell\Rightarrow z=0,\;\;x\rightarrow\infty\Rightarrow z\rightarrow 1 $$

Which gives $$ I = \frac{\pi\sigma}{2}\int_{0}^1 f\left(\sigma\tan\left(\frac{\pi z}{2}\right)+\ell\right) \sec^2\left(\frac{\pi z}{2}\right) \;dz $$ That integrand is a well-behaved enough for any numerical integration scheme of your choosing (the function is strictly positive, continuous, has finite first derivatives on the interval, and the tail decays faster than the secant blows up so the integrand is zero at $z=1$). I like the adaptive approaches where you approximate on $N$ points, then on $2N$ points, then $4N$ points, etc., and look at the differences between successive approximations; when the differences are smaller than some small tolerance value, e.g., $10^{-6}$ or so, then you can stop and error is less than the tolerance.

I'll finally note that I do not believe this would have any other closed form, and if it has some expression in terms of more exotic special functions, you'd still have to find a code that implements that special function, and under-the-hood you'd see it does various approximation techniques based on polynomial fits of the function (unlike the adaptive quadrature approach I suggested).