[R] how R implement qnorm()

[Prof Brian Ripley]

On 18/10/2012 00:16, Sheng Liu wrote:
> how R implement qnorm()
>
> I wonder anyone knows the mathematical process that R calculated the
> quantile?

It's on the help page!

'For qnorm, the code is a C translation of

Wichura, M. J. (1988) Algorithm AS 241: The Percentage Points of the 
Normal Distribution. Applied Statistics, 37, 477–484.'

and the exact code is in the R sources.  E.g. online at 
https://svn.r-project.org/R/trunk/src/nmath/qnorm.c

>
> The reason I asked is soly by curiosity. I know the probability of a normal
> distribution is calculated through integrate the Gaussian function, which
> can be implemented easily (see code),  while the calculation of quantile
> (or Zα) in 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).

> qnorm((1+x)/2)/sqrt(2) ). When you type qnorm in the console, it doesn't
> show it as it is an internal function, I searched around can't found too
> much information, my hunch is R might be using some internal library that's
> in the chipset which can calculate erf-1(x), but it is not accessible to
> user.

No, FPUs do not commonly have inverse erf functions.

>
> Any information is welcomed. thanks.
>
> Sheng
>
>
>
> code for implementation of pnorm()
> --------------------------------------------------
>
> p.Gaussian=function (z, mean=0,sd=1) {
>
> Gaussian=function(x) {1/(sqrt(2*pi)*sd)*exp(-(x-mean)^2/(2*sd^2))}
> per=integrate(Gaussian,lower=-Inf,upper=z)
>
> return (per$value)
> }
>
>
> code for implementation of qnorm()
> --------------------------------------------------
> # I've figured out one that uses the uniroot function to get x, it
> approximate qnorm() well but not exactly. I would be very happy to see the
> implementation through a mathematical formula such as using the X = -
> sqrt(2)* erf-1 (2*P), P is the probability).
>
> q.Gaussian=function(p,mean=0,sd=1) {
>     variable = function(x) p.Gaussian(x)-p
>     z = uniroot(variable, interval=c(-4,4))
>     v = z$root*sd+mean
>     return(v)
>
> }
>
> 	[[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.
>


-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595