[R] how R implement qnorm()
Prof Brian Ripley
ripley at stats.ox.ac.uk
Thu Oct 18 09:55:47 CEST 2012
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
More information about the R-help
mailing list