Normal distribution on 32E and 21S

[Dave Britten]

Does anybody happen to know what algorithms are used internally on the 32E and 21S to calculate cumulative normal and inverse normal? The 32E manual only says "The HP 32E computes Q^-1 using a sophisticated iterative algorithm", and I couldn't really find any details in the 21S manual. It sort of leaves to the imagination how many good digits to expect from the result. The various Application Pacs and books have typically used polynomial approximations, but I'm wondering if either of these models integrate numerically, and if so, to what accuracy. Or perhaps it's something else entirely.

[Albert Chan]

(01-29-2021 05:46 PM)Dave Britten Wrote:  The various Application Pacs and books have typically used polynomial approximations, but I'm wondering if either of these models integrate numerically, and if so, to what accuracy. Or perhaps it's something else entirely.

Q(x) is likely implemented using rational polynomial approximation of erf(x).

see http://git.musl-libc.org/cgit/musl/tree/src/math/erf.c

If not, we still do not need to numerically integrate.
erf integral can be turn into summation (or continued fraction form, if |x| is big)

see Abramowitz & Stegun, Error Function and Fresnel Integrals, 7.1.5 , 7.1.15
see XCas implementation, translated to Python

Quote:The 32E manual only says "The HP 32E computes Q^-1 using a sophisticated iterative algorithm"

This may be Newton's method: Q^-1(p) = z       → Q(z) - p = 0

>>> from mpmath import *
>>> Q = lambda x: erf(x/sqrt(2)) / 2
>>> pdf = lambda x: exp(-x*x/2) / sqrt(2*pi)
>>>
>>> p = 0.2
>>> z = 2.6*p                   # guess for Q^-1(p)
>>> z -= (Q(z) - p)/pdf(z); print z
0.524395467595098
>>> z -= (Q(z) - p)/pdf(z); print z
0.524400512701367
>>> z -= (Q(z) - p)/pdf(z); print z
0.524400512708041
>>> print Q(z)                   # should recover p
0.2

Comment: I used Casio FX-115MS definition: Q(z) = ∫(pdf(x), x = 0.5 .. z)
Casio also named P(z) = cdf(z), R(z) = 1 - cdf(z) = erfc(z/sqrt(2)) / 2

With Casio definition: P^-1(0.7) = Q^-1(0.2) = R^-1(0.3) = 0.524400512708041

Just noticed statistical Q-function = Casio's R function. Sorry about the confusion.

[Paul Dale]

I'd expect something like Newton's method for Q^-1. Calculating the CDF gives the PDF at the same time, so the first derivative is free. Likely, there will be an attempt to get a decent initial estimate and possibly some bounding or correction applied to avoid shooting off somewhere where it shouldn't be.

A higher order method is also viable. Second and higher order derivatives also come for free. Thus, Halley's method is also possible. This quickly becomes a balancing game between code size and speed of convergence.


Pauli

[Gerson W. Barbosa]

Search at TOS for

“Q function of HP-21S unveiled”

and

“Q function of HP-32E unveiled”

by former forum member Mike (Stgt).

[rprosperi]

You beat me Gerson!!

Mike has written quite a bit on these functions, with lots of theory and examples.

They can be found here:

http://forum.ho41.org/viewtopic.php?f=10&t=429 (HP-21S)

http://forum.ho41.org/viewtopic.php?f=10&t=427 (HP-32E)

Copy each address to your browser's entry window, replace "ho" with "hp" and hit enter.

[Albert Chan]

(01-29-2021 11:27 PM)Paul Dale Wrote:  I'd expect something like Newton's method for Q^-1. Calculating the CDF gives the PDF at the same time, so the first derivative is free. Likely, there will be an attempt to get a decent initial estimate and possibly some bounding or correction applied to avoid shooting off somewhere where it shouldn't be.

A higher order method is also viable. Second and higher order derivatives also come for free. Thus, Halley's method is also possible. This quickly becomes a balancing game between code size and speed of convergence.

Taking your advice, I tried a nice erf^-1(x) approximation, from Sergei Winitzki, 2008
A handy approximation for the error function and its inverse

As expected, Halley's method is not that complicated, because erf(x)'' / erf(x)' = -2x

Code:

from mpmath import *
def erf_inv(y, k=2/sqrt(pi), b=4.33, verbal=False):
    x = log1p(-y*y)/2               # max_rel_err = 0.002
    x = sign(y) * sqrt(sqrt(b*b+x*((2-pi)*b+x)) - (b+x))
    while 1:
        f = erf(x) - y
        h = -f/(k*exp(-x*x)+x*f)    # Halley correction
        if verbal: print repr(x), float(h)
        if x+h*h==x: return x+h
        x += h

icdf(p) = sqrt(2) * erf_inv(2*(p-.5))

>>> p = 0.7
>>> erf_inv(2*(p-.5), verbal=1)
mpf('0.37083921155017374') -3.20529566034e-05
mpf('0.37080715859357027') -1.2418308831e-14
mpf('0.37080715859355784')
>>> sqrt(2) * _ # icdf(0.7)
mpf('0.52440051270804067')

>>> p = 0.99
>>> erf_inv(2*(p-.5), verbal=1)
mpf('1.6433866997434172') 0.00158965243267
mpf('1.6449763521760914') 4.95709478826e-09
mpf('1.6449763571331861')
>>> sqrt(2) * _ # icdf(0.99)
mpf('2.3263478740408399')

[Dave Britten]

(01-30-2021 02:31 AM)rprosperi Wrote:  You beat me Gerson!!

Mike has written quite a bit on these functions, with lots of theory and examples.

They can be found here:

http://forum.ho41.org/viewtopic.php?f=10&t=429 (HP-21S)

http://forum.ho41.org/viewtopic.php?f=10&t=427 (HP-32E)

Copy each address to your browser's entry window, replace "ho" with "hp" and hit enter.

Perfect, thanks guys! Presumably then it just uses some iterative solver to compute the inverses. And Mike was even kind enough to do a 42S version already.

[trojdor]

Dave,
If you want to see some of the 'behind the curtain' notes on the algorithm adapted for use on the 32e, it's:
Algorithm AS 66: The Normal Integral
Author(s): I. D. Hill
Source: Journal of the Royal Statistical Society. Series C (Applied Statistics), Vol. 22, No. 3
(1973), pp. 424-427
(which in turn is based on Adams (1969))

Here's a link to the Wiley/RSS paper:
http://csg.sph.umich.edu/abecasis/gas_po...tegral.pdf

[Dave Britten]

(02-04-2021 04:54 PM)trojdor Wrote:  Dave,
If you want to see some of the 'behind the curtain' notes on the algorithm adapted for use on the 32e, it's:
Algorithm AS 66: The Normal Integral
Author(s): I. D. Hill
Source: Journal of the Royal Statistical Society. Series C (Applied Statistics), Vol. 22, No. 3
(1973), pp. 424-427
(which in turn is based on Adams (1969))

Here's a link to the Wiley/RSS paper:
http://csg.sph.umich.edu/abecasis/gas_po...tegral.pdf

Thanks, I'm no export on Fortran, nor calculus for that matter, so we'll see how much sense I can make of it. I see a bunch of coefficients, so I'm assuming it's some kind of series expansion.

As an aside, anybody happen to know how the 21S handles the F distribution and its inverse? All of the various HP application pacs seem to require that one of the degrees of freedom be even in their implementations of it, but the 21S imposes no such limitation.

[deffunc]

(02-04-2021 05:40 PM)Dave Britten Wrote:  As an aside, anybody happen to know how the 21S handles the F distribution and its inverse? All of the various HP application pacs seem to require that one of the degrees of freedom be even in their implementations of it, but the 21S imposes no such limitation.

I don't know which algorithm is used by the 21S, but a while back I had need for the t-distribution on my 27S. I implemented the continued fraction algorithm for the incomplete beta function recommended in chapter 6 of Numerical Recipes, 2nd Ed., for Q(t) and used Solve when I needed the inverse. Both the t and F distributions are special cases of the incomplete beta function. The method is slow to converge when the degrees of freedom are large, but mine were under 100 so that was never a problem. The 21S was restricted to dof < 200, if I remember correctly.

Numerical Recipes can be read online for free, and the 3rd Ed. explains how to use Halley's method for fast inverse. There is also Abramowitz and Stegun on their site, which lists the series and continued fraction formulae, and also polynomial and rational function approximations.

http://numerical.recipes/oldverswitcher.html