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 <[email protected]> wrote: > > On Oct 19, 2012, at 06:05 , Thomas Lumley wrote: > > > On Fri, Oct 19, 2012 at 12:21 PM, Sheng Liu <[email protected]> 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 <[email protected]> > 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: [email protected] Priv: [email protected] > >>> > >>> ______________________________________________ > >>> [email protected] 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]] > >> > >> ______________________________________________ > >> [email protected] 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 > > > > ______________________________________________ > > [email protected] 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: [email protected] Priv: [email protected] > > > > > > > > > [[alternative HTML version deleted]] ______________________________________________ [email protected] 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.
