Normal distribution on 32E and 21S

01292021, 05:46 PM
Post: #1




Normal distribution on 32E and 21S
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.


01292021, 09:51 PM
(This post was last modified: 01302021 03:30 AM by Albert Chan.)
Post: #2




RE: Normal distribution on 32E and 21S
(01292021 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.musllibc.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 FX115MS 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 Qfunction = Casio's R function. Sorry about the confusion. 

01292021, 11:27 PM
Post: #3




RE: Normal distribution on 32E and 21S
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 

01302021, 12:30 AM
Post: #4




RE: Normal distribution on 32E and 21S
Search at TOS for
“Q function of HP21S unveiled” and “Q function of HP32E unveiled” by former forum member Mike (Stgt). 

01302021, 02:31 AM
Post: #5




RE: Normal distribution on 32E and 21S
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 (HP21S) http://forum.ho41.org/viewtopic.php?f=10&t=427 (HP32E) Copy each address to your browser's entry window, replace "ho" with "hp" and hit enter. Bob Prosperi 

01302021, 01:13 PM
Post: #6




RE: Normal distribution on 32E and 21S
(01292021 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. 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 * icdf(p) = sqrt(2) * erf_inv(2*(p.5)) >>> p = 0.7 >>> erf_inv(2*(p.5), verbal=1) mpf('0.37083921155017374') 3.20529566034e05 mpf('0.37080715859357027') 1.2418308831e14 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.95709478826e09 mpf('1.6449763571331861') >>> sqrt(2) * _ # icdf(0.99) mpf('2.3263478740408399') 

01302021, 01:50 PM
Post: #7




RE: Normal distribution on 32E and 21S
(01302021 02:31 AM)rprosperi Wrote: You beat me Gerson!! 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. 

02042021, 04:54 PM
(This post was last modified: 02042021 05:05 PM by trojdor.)
Post: #8




RE: Normal distribution on 32E and 21S
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. 424427 (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 ENTER > = 

02042021, 05:40 PM
Post: #9




RE: Normal distribution on 32E and 21S
(02042021 04:54 PM)trojdor Wrote: Dave, 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. 

02082021, 10:33 AM
Post: #10




RE: Normal distribution on 32E and 21S
(02042021 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 tdistribution 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 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)