[R] how R implement qnorm()

[peter dalgaard]

On Oct 19, 2012, at 06:05 , Thomas Lumley wrote:

> On Fri, Oct 19, 2012 at 12:21 PM, Sheng Liu <sheng.liu at live.ca> wrote:
>> Thanks a lot. It's very helpful.
>> I've read through the c code. Just FYI and for my completion of the
>> question, I post some of my thought on it:
>> To me it looks like the algorithm is actually inquiring through an
>> approximation table (the approximations, at least for pnom, is "derived
>> from those in "Rational Chebyshev approximations for the error function" by
>> W. J. Cody, Math. Comp., 1969, 631-637."), it is not using any function
>> such as intergration or inverse erf to compute the value based on an
>> equation as I thought, though lack a bit subtlety but the up side is the
>> speed.
> 
> It's not an approximation table, it's a rational Chebyshev
> approximation, a ratio of polynomials.  qnorm also uses a rational
> polynomial approximation.

I inferred that Sheng probably knew that and meant that there's a table of coefficients of a rational polynomial approximation.

> At some level this *has to be* what is going on: computers don't
> implement pnorm/qnorm or erfc in hardware, and there is no closed-form
> expression for them in terms of quantities that are implemented in
> hardware, so the functions must be some sort of approximate expression
> using arithmetic and other hardware computations.

Up to a few years ago, or maybe a decade, that was also the way special functions and even square roots and division were implemented in hardware using lookup tables and add/multiply circuits. All hell broke loose when one of the coefficients got entered incorrectly -- some may remember the FDIV Pentium bug of 1994. I used to have a really nice paper from Byte magazine c. 1990 which detailed the process, but I think it got lost in space. (L. Brett Glass: "Math Coprocessors" could be the one. http://dl.acm.org/citation.cfm?id=86680). 

Nowadays, we have single CPU cycle transcendentals, so I suppose that has all been reduced to hard-core optimized electronic circuitry. Fringe-market customers like statisticians still need to implement their functions in software, of course. 

>> From the point of view of R, the only relevant issues are precision,
> speed, and portability, and these approximations are portable,
> accurate, and fast.
> 
>    -thomas
> 
>> Hope this can help others who had similar questions as well.
>> 
>> Sheng
>> 
>> 
>> On Thu, Oct 18, 2012 at 2:32 AM, peter dalgaard <pdalgd at gmail.com> wrote:
>> 
>>> 
>>> On Oct 18, 2012, at 09:55 , Prof Brian Ripley wrote:
>>> 
>>>>> R is a bit confusing as it requires inverse error function (X =
>>>>> - sqrt(2)* erf-1 (2*P)), while R doesn't have a build in one. The InvErf
>>>>> function most people use is through qnorm( InvErf=function(x)
>>>> 
>>>> I think you are wrong about 'most people': this is the notation used by
>>> a small group of non-statisticians (mainly physicists, I think).
>>> 
>>> Well, he's right in the sense that it is erf/erfc that are commonly
>>> documented in collections of special functions (like Abramowitz & Stegun),
>>> and also those that appear in common C math libraries. In both cases of
>>> course because physicists have dominated in their development.
>>> 
>>> Of course "most people" use Excel these days (which has the inverse normal
>>> distribution but AFAICS not the inverse error function).
>>> 
>>> --
>>> Peter Dalgaard, Professor,
>>> Center for Statistics, Copenhagen Business School
>>> Solbjerg Plads 3, 2000 Frederiksberg, Denmark
>>> Phone: (+45)38153501
>>> Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com
>>> 
>>> ______________________________________________
>>> R-help at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide
>>> http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>>> 
>> 
>>        [[alternative HTML version deleted]]
>> 
>> ______________________________________________
>> R-help at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
> 
> 
> 
> -- 
> Thomas Lumley
> Professor of Biostatistics
> University of Auckland
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com