[R] how R implement qnorm()

[Thomas Lumley]

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.

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.

>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