[tdyuocul] expressing low probabilities as sigmas

for a given probability p, where (what z score) on the bell curve (standard normal distribution, area of the left tail) corresponds to that probability?

using Octave or Matlab, the answer is z = norminv(p) .  norminv is also known as the probit function.

consider using z-score as a general means of expressing a small probabilities, even if the probability did not come from the normal distribution.  (though a great many things are normally distributed, or approximately so.)  this was inspired by particle physics which denotes significance of an experimental result in terms of "sigmas".  this is the same as z score, the number of standard deviations away from the mean.

we describe how to calculate z if p is smaller than that which can be represented by floating point.  because p is very small, we assume the input probability will be given as a natural logarithm: logp = log(p).  For example, if p = 10^-1000, then logp = -1000*log(10) ~= -2302.6.

we will use the scaled complementary error function erfcx (available in Matlab and Octave, but not Mathematica) to further avoid floating point underflow.  erfcx(x) "factors out" the rapidly shrinking exp(-x^2) portion of the complementary error function erfc(x) = 1-erf(x).

we need to solve z^2/2 - log(erfcx(-z/sqrt(2))/2) + logp = 0.

to make an initial guess, we ignore the log(erfcx) term, assuming it is small compared to z^2/2.  a quick and dirty initial guess is then
z = -sqrt(-2*logp)

digression: erfcx(x) can be approximated as erfcx(x) ~= 1/(sqrt(pi)*x) for large x, but unfortunately that does not get us anywhere: z^2/2 + log(-z) + log(sqrt(2*pi)) + logp = 0 .  if we could quickly solve that, then we would have a better initial guess for z.  however, that equation doesn't seem solvable in closed form.  update: Mathematica gives z = (+/-) Sqrt[ ProductLog[1 / (2 * E^(2*logp) * Pi)]] .  ProductLog is Mathematica's name for the Lambert W function.

given the initial guess of z, we can improve the solution with Newton's method.  the derivative of erfcx(x) is 2*(x*erfcx(x) - 1/sqrt(pi)) .  the Newton's iteration is

z = z + (z^2/2-log(erfcx(-z/sqrt(2))/2)+logp)*erfcx(-z/sqrt(2))/sqrt(2/pi)

here are some computational results.  the less extreme values can be verified with norminv.

p=0.1, z=-1.2816
p=0.01, z=-2.3263
p=0.001, z=-3.0902
p=10^-10, z=-6.3613
p=10^-100, z=-21.2735
p=10^-1000, z=-67.7857

p=2^-32, z=-6.2303
p=2^-64, z=-9.0802
p=2^-128, z=-13.0559
p=2^-256, z=-18.6332
p=2^-512, z=-26.4837
p=2^-1024, z=-37.5563
p=2^-82589933, z=-10700.18

quadruple the exponent, double the z score.

probabilities that are large negative powers of 2 are seen in cryptanalysis, describing things like the probability of a successful attack.

the final line is p=1/(m+1), where m is the currently largest known prime.  even that astronomically low probability corresponds to a reasonably describable location under the bell curve, only ten thousand standard deviations below the mean.