how R implement qnorm()
I wonder anyone knows the mathematical process that R calculated the
quantile?
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)
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.
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.
