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]

Sqrt[2] InverseErf[0.9545]

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]


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.

  • 2
    $\begingroup$ I assume this is because Mathematica isn't in fact calculating anything, but is simply taking this expansion from its tabulated data. $\endgroup$
    – Domen
    Jul 22, 2022 at 12:38
  • $\begingroup$ Ok, I wasn't aware of this, thanks! $\endgroup$
    – Hans Olo
    Jul 22, 2022 at 12:52
  • 1
    $\begingroup$ While that data is in the functions.wolfram.com domain, it is not where the problematic series (lack of) expansion arises. $\endgroup$ Jul 22, 2022 at 18:32
  • 2
    $\begingroup$ 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. $\endgroup$ Jul 22, 2022 at 20:10
  • 3
    $\begingroup$ 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. $\endgroup$ Jul 25, 2022 at 16:20

1 Answer 1


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.

  • 1
    $\begingroup$ 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$. $\endgroup$
    – Hans Olo
    Jul 23, 2022 at 15:30
  • 3
    $\begingroup$ This can be streamlined with SolveAlways[]: sol = yansatz[η] /. First[SolveAlways[eq2 == 0, {logeta, inveta}]] $\endgroup$ Jul 23, 2022 at 15:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.