# How to get more terms with the Series[] expansion of InverseErf[x] around x=1?

The inverse error function, which is given by InverseErf[x], is quite important in statistics as it gives the confidence levels around a 1D Gaussian, for example for the $$1$$ and $$2\sigma$$ ($$\sim 68.3\%$$ and $$\sim95.5\%$$ respectively):

Sqrt[2] InverseErf[0.6827]
=1

Sqrt[2] InverseErf[0.9545]
=2


So, while playing with it and wanting to find some asymptotics close to probability 1 (i.e. many $$\sigma$$'s), I did a series expansion around $$x=1$$ (note in hindsight I'm doing the expansion of InverseErf[x]^2 to avoid a square-root when comparing with the result of a paper, see below) and no matter what order I asked for

Series[InverseErf[x]^2, {x, 1, -30}, Assumptions -> x < 1]


or

Series[InverseErf[x]^2, {x, 1, 30}, Assumptions -> x < 1]


Mathematica always gives back the result:

$$\frac{1}{2} \left(-2 \log (1-x)-\log (-2 \log (1-x)+\log (2)-\log (\pi ))+\log \left(\frac{2}{\pi }\right)\right)+O\left((x-1)^2\right).~~~~~~~~~~~~(1)$$

However, the inverse error function has a well defined (albeit asymptotic) expansion around $$x=1$$, see for example the seminal paper by Blair et al on approximations of the error function, where in Eq.(2) they give the expression: $$(\textrm{InverseErf}[x])^2 \sim \eta-\frac12 \ln \eta +\eta^{-1} \left(\frac14 \ln \eta-\frac12\right)+\eta^{-2} \left(\frac1{16} \ln^2 \eta-\frac38 \ln \eta +\frac78\right)+...,~~~~~~~~~~~~(2)$$ where $$\eta=-\ln [\pi^{1/2}*(1-x)]$$ and in this notation $$\ln(x)=\log(x)=\log_e(x)$$ (ie the log base e). Also, $$\ln^2 \eta=[\ln(\eta)]^2$$ (thanks to @user293787).

Using the same variable $$\eta$$, Mathematica's series expansion, given originally by Eq. (1), can be written as

$$(\textrm{InverseErf}[x])^2 \sim \eta -\frac{1}{2} \ln \left(\eta +\frac{\ln (2)}{2}\right).~~~~~~~~~~~~(3)$$

So, my question is: why does Mathematica only stop at such low order in the expansion in terms of $$x$$ and how to coax it to give more terms?

Note 1: In fact the expression from the paper is much more accurate and by keeping up to the $$\eta^{-1}$$ term (inclusive) one can get an accuracy of $$0.00013\%$$ at $$x \rightarrow 1 - 10^{-30}$$ ($$\eta\sim68.5051878$$) with respect to the exact result, compared to the Series[] expansion which is "only" accurate to $$\sim 0.016 \%$$.

Note 2: I can in fact just code the expression from the paper and use that if I want or even use the exact arbitrary precision result from InverseErf[x], but the question is how to get Series to give more terms, mainly out of curiosity.

Note 3: The result from the paper also gives a complex $$i \pi$$ term, but I only consider the real part in the comparison.

Note 4: This question, only discusses the limit around $$x=0$$, which is already well-known via a MacLaurin series.

Note 5: Tested mainly on v11.2 on Windows but the same issue seems to persist in v13 as well.

• I assume this is because Mathematica isn't in fact calculating anything, but is simply taking this expansion from its tabulated data. Jul 22, 2022 at 12:38
• Ok, I wasn't aware of this, thanks! Jul 22, 2022 at 12:52
• While that data is in the functions.wolfram.com domain, it is not where the problematic series (lack of) expansion arises. Jul 22, 2022 at 18:32
• I mean that the Series code is not directly using that tabulated data. It is failing to expand for other reasons. This might be a bug in some way. Jul 22, 2022 at 20:10
• I filed a report about this by the way. It's a limitation I do not know how to address myself, but maybe someone else here will have some ideas for improving on what we now do. Jul 25, 2022 at 16:20

This is not strictly speaking an answer, but I thought I provide code that can be used to generate the kind of expansion that OP mentions. Perhaps it is useful for other people here.

I use $$y$$ as an abbreviation for InverseErf[x]^2. As a function of $$\eta$$ it satisfies the differential equation $$y'(\eta) = \sqrt{y} e^{y-\eta}$$ I will use this to determine the coefficients $$a_{km}$$ in the following ansatz for $$\eta \to \infty$$. This ansatz is motivated by the expansion that OP has given:

n = 6; (* order of the expansion *)
yansatz = Function[{η}, η - 1/2 * Log[η] +
Sum[1/η^k * Sum[a[k,m] * Log[η]^m, {m, 0, k}], {k, 1, n}]];


To plug this into the differential equation and solve for the coefficients, one can use the following code, where I am now using SolveAlways as per @J.M.'s suggestion (thanks!).

eq = {y'[η] - Sqrt[y[η]] * Exp[y[η] - η]};
eq2 = Normal[Series[eq /. {y -> yansatz}, {η, Infinity, n}]] /.
{Log[η] -> logeta} /. {η -> 1/inveta};
yansatz[η] /. First[SolveAlways[eq2 == 0, {logeta, inveta}]]


For example, with n = 6 this gives

The first few terms are consistent with the expansion given by OP.

• Very usuful answer, thanks! As the expression perfectly matches the one in the paper, it also resolves the question of @MichaelE2 about the logarithm, ie $\ln^2(\eta)=[\ln(\eta)]^2$. Jul 23, 2022 at 15:30
• This can be streamlined with SolveAlways[]: sol = yansatz[η] /. First[SolveAlways[eq2 == 0, {logeta, inveta}]] Jul 23, 2022 at 15:51