3
$\begingroup$

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.

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ 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? $\endgroup$
    – Kirill
    Oct 6, 2016 at 23:53
  • $\begingroup$ For example: arxiv.org/abs/1407.7786 $\endgroup$
    – Kirill
    Oct 6, 2016 at 23:59
  • $\begingroup$ As a start, we can take $p = 32$ and $z \in [-4, 4]$. $\endgroup$
    – Pine Tree
    Oct 7, 2016 at 0:42
  • $\begingroup$ That's pretty small; have you tried just evaluating the recurrence directly, or are you saying that you did and it failed (how?)? $\endgroup$
    – Kirill
    Oct 7, 2016 at 19:39
  • $\begingroup$ @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. $\endgroup$
    – njuffa
    Oct 25, 2016 at 4:16

1 Answer 1

Reset to default
5
$\begingroup$

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.

Simple implementation of recurrence

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. Behaviour of $f_p$ as a function of the parameter $p$ for fixed $z=-4$

Plotting the relative error gives: Relative error for $f_p$ using the recurrence for fixed $z=-4$

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:

Integrand defining $f_p$ for $p=11$

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: enter image description here

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...

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.