Φ(−600) = 6.546588 205807 692852 105927 713888 108782 119412 831853 177211 16943e − 78177

Specific Domains Statistics
1

The paper Approximations to the Normal Distribution Function and An Extended Table for the Mean Range of the Normal Variables by Kiani, M and Panaretos, J and Psarakis, S and Saleem, M (2008) reviews algorithms for computing the CDF of the univariate Normal distribution. Then they develop some new approximations.

From use of their new approximations, we learn that:

Φ(−600) = 6.54658820580769285210592771388810878211941283185317721116943e − 78177

This required them to compute 1080 terms of their series approximation.

Now Julia:

julia> sn=Normal(BigFloat("0.0"),BigFloat("1.0"))
Normal{BigFloat}(μ=0.0, σ=1.0)

julia> cdf(sn,-600)
6.546588205807692852105927713888108782119412831853177211169430430310177386453905e-78177

How did Julia do it?

Best I can figure out, the univariate Normal CDF is computed using the error function from the GNU MFPR library. Apparently it just keeps going until all the bits of BigFloat accurate.

2

In Base Julia:

function my_erfc(x::BigFloat)
   z = BigFloat()
   ccall((:mpfr_erfc, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), z, x, Base.MPFR.ROUNDING_MODE[])
   return z
end

my_erfc(big(600.0)*inv(sqrt(big(2))))/2

(adapting code from IrrationalConstants, StatsFuns and SpecialFunctions)

3

Yup, and MPFR uses an asymptotic series for the complementary error function of large arguments (whenever x^2 \ge p \log2, where p is the precision), which converges fairly quickly in this regime.

image
image1120×356 32.6 KB

e.g. you get a decent approximation with only two terms:

julia> my_erfc_asympt(z) = (1 - 1/2z^2) * exp(-z^2) / (√π*z)
my_erfc_asympt (generic function with 1 method)

julia> setprecision(52)
52

julia> my_erfc_asympt(big(600.0)*inv(sqrt(big(2))))/2
6.5465882048940283e-78177
4

Thank you for posting the expression used. I had tried to find it, but never found the right path to drill down to it.