Thanks, again! Sorry for my misleading expression, I only knew the value is
inquired from a polynomial approximation, but I have no idea how it is done
in such a great detail. It's a great lesson for people like me who want a
deep understanding of the basics.
Sheng
On Thu, Oct 18, 2012 at 11:46 PM, peter dalgaard wrote:
>
> On Oct 19, 2012, at 06:05 , Thomas Lumley wrote:
>
> > On Fri, Oct 19, 2012 at 12:21 PM, Sheng Liu 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
> 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@cbs.dk Priv: PDalgd@gmail.com
> >>>
> >>> ______________________________________________
> >>> R-help@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@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@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@cbs.dk Priv: PDalgd@gmail.com
>
>
>
>
>
>
>
>
>
[[alternative HTML version deleted]]