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@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, rip...@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
______________________________________________
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.
