Calculating the z-score to arbitrary precision for a given percentile

I am trying to calculate the value for a given distribution that is greater than a certain percentile of other elements of that distribution. For example, for a distribution with mean X and standard deviation Y, what is the value T for which 95% of the distribution will be less than T and 5% will be greater than T. This can be done if one knows the corresponding z-score (T = X + z-score*Y). However, tables are incomplete(multiple z-scores correspond to the same percentile etc) and I'd like to be able to do this to higher accuracy than the table allows.

Is it the case that in order to calculate z for an arbitrary distribution I have to solve $$0.95 = \frac{1}{\sigma\sqrt{2\pi}}\int_0^z e^{(-(x-\mu)/\sigma)^2}dx$$ for z? This would involve some messy error function results, correct? Does anyone have any tips about solving this numerically?

• This has probably already been done in some arbitrary precision computing library (although probably through erf and its inverse, so you would need a few other arbitrary precision operations as well). For instance, upon a quick Google search I find sciencedirect.com/science/article/pii/S0890540112000697 which implements erf and erfc (the latter actually being what you want here). Imitating this from scratch will be more difficult than you expect. Also, your lower bound for the problem you stated should be $-\infty$. – Ian Sep 15 '17 at 21:13
• I wonder what the application is? I'm sure there are good and interesting reasons to compute a normal integral to arbitrary precision. Statistical application is not one of them. – Mark McClure Sep 21 '17 at 16:48

You say you are interested in a 'given' distribution, but as your question progresses, you talk about z-scores and normal distributions. The CDF $\Phi$ of a standard normal distribution cannot be expressed in closed form. However, various software packages have inverse CDF functions, called 'quantile' functions, that provide the answers you're looking for.

For example, if $X \sim \mathsf{Norm}(\mu = 100,\, \sigma = 15),$ then you can find $q_{.95} = 124.6728.$ such that $P(X \le q_{.95}) = 0.95000.$ In R statistical software this is computed as shown below. In R, the quantile function is qnorm and the CDF is pnorm; arguments $\mu$ and $\sigma$ can be supplied, so that it is not necessary to standardize.

q.95 = qnorm(.95, 100, 15);  q.95;  pnorm(q.95, 100, 15)
## 124.6728   # 95th quantile of NORM(100, 15)
## 0.95       # check


The 95th quantile of the standard normal distribution is 1.644854:

qnorm(.95)
## 1.644854


The quantiles shown above are sufficiently accurate for any practical purpose I can imagine, but if you have a theoretical interest in more accuracy, it should be possible to find software that gives more digits of accuracy.

The quantile function in R uses a rational approximation to $\Phi^{-1}$ (by Wichura) that is accurate to about as many places as the floating point arithmetic in R can express. More decimal places than given above could be found by iterative numerical integration. Although the quantile function cannot be expressed in closed form, it is known to be continuous, so there is no theoretical limit to the number of digits of accuracy obtainable.

Note: Most statistical computer packages implement CDFs and quantile computations for some commonly-used distributions other than normal. Below are 95th percentiles for $\mathsf{Exp}(1)$, which has CDF $F(x) = 1 - e^{-x},$ for $x > 0;$ and for $\mathsf{Unif}(0,1),$ which has CDF $F(x) = x,$ for $0 < x < 1.$ These follow the same kind of syntax in R as for the normal distribution. However, in these cases the CDF is known and invertible, so computation is straightforward.

q.e = qexp(.95);  q.e;  1 - exp(-q.e)
## 2.995732
## 0.95
qunif(.95)
## 0.95


I ended up using the inverse error function as specified in the boost special function library (I'm working in C++). That was able to get me what I wanted without having to numerically integrate.

Thanks for the responses.