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

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.

• Since the integral can be expressed in closed form in terms of confluent hypergeometric functions as $$2^{\frac{p}{2}-1}\Gamma(\tfrac{p}{2})\mathrm{F}(\tfrac{1-p}{2}; \tfrac{1}{2}; -\tfrac{1}{2}z^2) + 2^{\frac{p-1}{2}}\Gamma(\tfrac{p+1}{2})\mathrm{F}(1-\tfrac{p}{2}; \tfrac{3}{2}; -\tfrac{1}{2}z^2),$$ is there an existing library that can evaluate $\log({}_1\mathrm{F}_1)$? Failing that, check the references at the DLMF. Also, the question is likely difficult for very large $p$ or $z$, so can you say more precisely what the bounds on them are? – Kirill Oct 6 '16 at 23:53
• For example: arxiv.org/abs/1407.7786 – Kirill Oct 6 '16 at 23:59
• As a start, we can take $p = 32$ and $z \in [-4, 4]$. – Pine Tree Oct 7 '16 at 0:42
• That's pretty small; have you tried just evaluating the recurrence directly, or are you saying that you did and it failed (how?)? – Kirill Oct 7 '16 at 19:39
• @PineTree How is "accurately" defined, and is it more important that $f_p(z)$ is accurate or that $\log f_p(z)$ is accurate? For example, if "accurate" is defined as relative error, then $\log f_2(z)$ will lose accuracy for $\lvert z \rvert \lt \frac{1}{4}$ as $f_2(z)$ is near unity, even if $f_2(z)$ itself is computed accurately. Conversely, since $\log$ is a contracting operation, we could tolerate some amount of error in $f_p(z)$, say 30 ulp, while still getting an accurate $\log f_p(z)$ with only about 5 ulp errors. – njuffa Oct 25 '16 at 4:16

## 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):
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...

• (+1) I just want to say, in my comment above the hypergeometric functions are both positive so there wouldn't be any cancellation issues. – Kirill Oct 11 '16 at 18:28