Numerical evaluation of gaussian-like integral expressible as a recurrence relation

Question

I'm looking to numerically evaluate $\log f_p(z)$ and its derivative $f^\prime_p(z)/f_p(z)$ accurately and efficiently in floating-point, where

$$ f_p(z)=\int_0^\infty r^{p-1} \exp\left(-\tfrac{1}{2} (r-z)^2\right)\,dr \ . $$

Let $\phi(x)$ and $\Phi(x)$ be the PDF and CDF of the standard normal distribution $x \sim \mathcal{N}(0,1)$ respectively. Then the integral $f_p(z)$ can be expressed as the recurrence relation

\begin{align} f_1(z) &= \sqrt{2\pi}\,\Phi(z) \\ f_2(z) &= \sqrt{2\pi}\,\phi(z) + z\,f_1(z) \\ f_p(z) &= (p-2)\,f_{p-2}(z) + z\,f_{p-1}(z)\ . \end{align}

How can I evaluate $\log f_p(z)$ and its derivative $f^\prime_p(z)/f_p(z)$ in floating point efficiently and accurately? The parameter $p$ can be assumed to be fixed such that extra pre-computation can be performed for each $p$. If needed, I can also restrict $z$ to an interval $z \in [c_1,c_2]$.

A direct evaluation by precomputing the polynomial coefficients of $\phi(z)$ and $\Phi(z)$ is feasible, but the coefficients involved quickly become unmanageably huge and would require extended-precision intermediate values even though the integrand and result are both well-behaved.

Answer 1

Step 1: analyze the recurrence

I implemented the recurrence in Python, using basic numpy and scipy.special for the erf function. Why? Because it is simple and cheap (for reasonable $p$).

import numpy as np
import scipy.special as sp
def F1(z):
    return np.sqrt(np.pi/2)*(1.0+sp.erf(z/np.sqrt(2)))

def F2(z):
    return F1(z)*z + np.exp(-z**2/2.0)

def F(p,z):
    if (p==1):
        return F1(z)
    elif (p==2):
        return F2(z)
    else:
        fpm2 = F1(z)
        fpm1 = F2(z)

        for pi in range(3,p+1):
            fp = (pi-2)*fpm2 + z* fpm1
            fpm2 = fpm1
            fpm1 = fp
        return fp

I compared results from this to fixed precision results using Maple 2016 (Digits set to 500). The following graph shows the relative error between the recurrence results and Maple results.

You can clearly see that the recurrence is forward-stable (stable in increasing $p$) for positive $z$. For negative $z$, things are different. For $z<0$, $f_p$ as a function of the parameter $p$ for a fixed $z$ has a bump (it is not monotonous). See next graph.

Plotting the relative error gives:

One observes a loss of accuracy for increasing $p$. For $p$ above 27, you can see that the log plot of $f_p$ does not continue, simply because the recurrence gives negative values.

Step 2: can we do black-box numerical integration?

Again, for $z>0$, this would be meaningless since the recurrence works so well. We could try it for negative values of $z$. So let's have a look at the integrand:

Positive point: the integrand decreases very rapidly (after a peak). So one could approximate the infinite range integration to a finite range based on the value of the integrand.

Problematic point: the function is enormously peaked. Let's see if QUADPACK (used as back-end in scipy.integrate.quad) can handle this. Note that I put the right point of integration arbitrarily at $r=200$. This could be optimized (if need be).

from scipy.integrate import quad
def F_neg(p,z):
    Qres = np.zeros_like(z)
    for i,zz in enumerate(z):
        Q,err = quad(integrand,0,200,args=(float(p),float(zz)),epsrel=1e-12,epsabs=1e-12)
        Qres[i] = Q
    return Qres

Again, a plot of the relative error comparing this Python code to Maple results:

Perfect!

Preliminary conclusion: use the recurrence for $z>0$, use the numerical quadrature for $z<0$.

I don't think that using closed form solutions using hypergeometric functions (or LaguerreL functions) will help much since you will have linear combinations of these, possibly leading to cancellation problems.

Work to be done: check range of $z$ and $p$ for which the above holds. But that is for the poster of the question...