Analytic error function and numeric inverse obtained by geometric means

TL;DR

This paper presents a geometric approach to deriving the error function's integral representation and its inverse, offering new formulas and insights useful for high-speed Monte Carlo simulations.

Contribution

It introduces a geometric method to derive the error function's integral form and inverse approximation formulas, enhancing computational efficiency.

Findings

Derived the Craig formula using geometric considerations

Proved convergence of the power series expansion

Developed systematic formulas for approximating the inverse error function

Abstract

Using geometric considerations, we provide a clear derivation of the integral representation for the error function, known as the Craig formula. We calculate the corresponding power series expansion and prove the convergence. The same geometric means finally help to systematically derive handy formulas that approximate the inverse error function. Our approach can be used for applications in e.g.\ high-speed Monte Carlo simulations where this function is used extensively.

Tables

Table 1. Table 1. Run time experiment for our algorithm under C for E = 0.8 𝐸 0.8 E=0.8 and different values of n 𝑛 n and p 𝑝 p (CPU time in seconds). As indicated, the errors are in the last displayed digit, i.e., ± 0.01 plus-or-minus 0.01 \pm 0.01 seconds.

	$n=0$	$n=1$	$n=2$	$n=3$	$n=4$	$n=5$
$p=0$	$0.07\left(1\right)$	$0.13\left(1\right)$	$0.17\left(1\right)$	$0.21\left(1\right)$	$0.31\left(1\right)$	$0.56\left(1\right)$
$p=1$	$0.14\left(1\right)$	$0.20\left(1\right)$	$0.24\left(1\right)$	$0.29\left(1\right)$	$0.39\left(1\right)$	$0.63\left(1\right)$
$p=2$	$0.25\left(1\right)$	$0.32\left(1\right)$	$0.35\left(1\right)$	$0.40\left(1\right)$	$0.50\left(1\right)$	$0.75\left(1\right)$

Table 2. Table 2. Results for p = 2 𝑝 2 p=2 and increasing values of n 𝑛 n for values of E 𝐸 E approaching E = 1 𝐸 1 E=1 . The last column shows the CPU time for 10 6 superscript 10 6 10^{6} runs according to the algorithm proposed in Ref.  [ 12 ] for the values of N 𝑁 N given in the last column of the table displayed in Ref.  [ 12 ] .

$E=$	$n=0$	$n=1$	$n=2$	$n=3$	$n=4$	$n=5$	[12]
$0.7$	$0.732995$	$0.732868$	$0.732869$	$0.732869$	$0.732869$	$0.732869$	$0.17$
$0.8$	$0.906326$	$0.906193$	$0.906194$	$0.906194$	$0.906194$	$0.906194$	$0.19$
$0.9$	$1.163247$	$1.163085$	$1.163087$	$1.163087$	$1.163087$	$1.163087$	$0.35$
$0.99$	$1.821691$	$1.821376$	$1.821387$	$1.821386$	$1.821386$	$1.821386$	$1.95$
$0.999$	$2.326608$	$2.326762$	$2.326752$	$2.326754$	$2.326754$	$2.326754$	$14.62$
$0.9999$	$2.749217$	$2.751197$	$2.751034$	$2.751076$	$2.751056$	$2.751971$	$128.30$

Equations

$$\mathrm{erf}(t)=\frac{1}{\sqrt{\pi }}{\int }_{-t}^t{e}^{-x^2}dx=\frac{1}{\sqrt{2\pi }}{\int }_{-\sqrt{2}t}^{\sqrt{2}t}{e}^{-x^2/2}dx=P(\sqrt{2}t).$$

$${\mathrm{erf}}_p(t{)}^2=1-\frac{1}{N}\sum _{n=1}^N{e}^{-k_{p,n}^2{t}^2},$$

$$\Big{|}\mathop{\rm erf}\nolimits(t)-\sqrt{1-e^{-k_{0,1}^{2}t^{2}}}\Big{|}<0.0033$$

$$\Big{|}\mathop{\rm erf}\nolimits(t)-\sqrt{1-\frac{1}{2}(e^{-k_{1,1}^{2}t^{2}}+e^{-k_{1,2}^{2}t^{2}})}\Big{|}<0.00024,$$

$$|\mathrm{erf}(t)-{\mathrm{erf}}_p(t)|<\frac{\exp (-t^2)}{2N}\sqrt{1-\exp (-t^2)}.$$

$$\mathrm{erf}(t{)}^2=1-\frac{4}{\pi }{\int }_0^{\pi /4}\exp \left(\frac{-t^2}{{\cos }^2\varphi }\right)d\varphi .$$

$$\mathrm{erf}(t{)}^2=1-\frac{4}{\pi }\sum _{n=0}^{\infty }{c}_n\frac{(-1{)}^n}{n!}(t^2{)}^n$$

$$c_n$$

$$\left(\begin{array}{c}n\\ k\end{array}\right)+\left(\begin{array}{c}n\\ k-1\end{array}\right)=\frac{n!}{k!(n-k)!}+\frac{n!}{(k-1)!(n-k+1)!}$$

$$\sum _{k=0}^n\left(\begin{array}{c}n\\ k\end{array}\right)=2^n$$

$$\sum _{k=0}^{n+1}\left(\begin{array}{c}n+1\\ k\end{array}\right)=1+\sum _{k=1}^n\left(\begin{array}{c}n+1\\ k\end{array}\right)+1=$$

$$\frac{1}{2n-1}\le \frac{1}{2k+1}\le 1.$$

$$\frac{1}{2n-1}\sum _{k=0}^{n-1}\left(\begin{array}{c}n-1\\ k\end{array}\right)\le c_n\le \sum _{k=0}^{n-1}\left(\begin{array}{c}n-1\\ k\end{array}\right)=2^{n-1},$$

$${\mathrm{erf}}_p(t{)}^2=1-\frac{4}{\pi }\sum _{n=0}^N{c}_n\frac{(-1{)}^n}{n!}(t^2{)}^n=1-\frac{2}{\pi }\sum _{n=0}^N{c}_n^{\ast }\frac{(-2t^2{)}^n}{n!},$$

$$\frac{2}{\pi }\sum _{n=0}^N{c}_n^{\ast }\frac{(-2t^2{)}^n}{n!}=c_N^{\ast \ast }\frac{2}{\pi }\sum _{n=0}^N\frac{(-2t^2{)}^n}{n!}<\frac{2}{\pi }\sum _{n=0}^N\frac{(-2t^2{)}^n}{n!}.$$

$$\mathrm{erf}(t{)}^2=\sum _{n=1}^{\infty }{c}_n\frac{(-1{)}^{n-1}}{n!}(t^2{)}^n,c_n=\sum _{k=0}^{n-1}\frac{1}{2k+1}\left(\begin{array}{c}n-1\\ k\end{array}\right).$$

$$T_0=\sqrt{-\ln (1-E^2)}/k_{0,1}^2.$$

$${\mathrm{erf}}_1(T_0)=\sqrt{1-\frac{1}{2}(e^{-k_{1,1}^2{T}_0^2}+e^{-k_{1,2}^2{T}_0^2})},$$

$$\frac{\mathop{\rm erf}\nolimits(t)-\mathop{\rm erf}\nolimits(T_{0})}{t-T_{0}}=\frac{\Delta\mathop{\rm erf}\nolimits(t)}{\Delta t}\Big{|}_{t=T_{0}}\approx\frac{d\mathop{\rm erf}\nolimits(t)}{dt}\Big{|}_{t=T_{0}}=\frac{2}{\sqrt{\pi}}e^{-T_{0}^{2}},$$

$$t=t_0+A(E-\mathrm{erf}(t_0))$$

$$t=t_0+A(E-\mathrm{erf}(t_0))+B(E-\mathrm{erf}(t_0){)}^2$$

$$\Delta t=A\Delta E_1+B\Delta E_1^2,2\Delta t=A\Delta E_2+B\Delta E_2^2$$

$$A=-\frac{(2\Delta E_1^2-\Delta E_2^2)\Delta t}{\Delta E_1\Delta E_2(\Delta E_1-\Delta E_2)},B=\frac{(-2\Delta E_1+\Delta E_2)\Delta t}{\Delta E_1\Delta E_2(\Delta E_1-\Delta E_2)}.$$

$$t=t_0+A(E-\mathrm{erf}(t_0))+B(E-\mathrm{erf}(t_0){)}^2+C(E-\mathrm{erf}(t_0){)}^3$$

$$A$$

$$B$$

$$C$$

Full text

MSC Class: 62E15, 62E17, 60E15, 26D15.

eestidima@gmail.com, stefan.groote@ut.ee

1. Introduction

High-speed Monte Carlo simulations are used for a large spectrum of applications from mathematics to economy. As input for such simulations, the probability distribution are usually generated by pseudo-random number sambling, a method going back to a work of John von Neumann from 1951 [1]. In the era of “big data” such methods have to be fast and reliable, a sign of such neccessity being the release of the very first randomness Quside processing unit in 2023 [2]. Still, these samblings need to be cross-checked by exact methods, and for these the knowledge of analytical functions to describe the stochastic processes, among those the error function, are of tremendous importance.

By definition, a function is called analytic if it locally given by a converging Taylor series expansion. Even if a function itself turns out not to be analytic, its inverse can be analytic. The error function can be given analytically, one of these analytic expressions is the integral representation given by Craig in 1991 [3]. Craig mentioned this representation only in passing and did not give a derivation of it. In the following, there have been a couple of derivations of this formula [4, 5, 6]. In Sec. 2 we add a further one which is based on the same geometric considerations as employed in Ref. [7]. In Sec. 3 we give the series expansion for Craig’s integral representation and show the fast convergence of this series.

For the inverse error function, handbooks for special functions (cf. e.g. Ref. [8]) do not unveil such an analytic property. Instead, this function have to be approximated. Known approximations are dating back to the late 1960s and early 1970s [9, 10]) and reach up to semi-analytical approximations by asymptotic expansion (cf., e.g., Refs. [11, 12, 13, 14, 15, 16]. Using the same geometric considerations, in Sec. 4 we develop a couple of handy approximations which can easily be implemented in different computer languages, indicating the deviations from an exact treatment. In Sec. 5 we discuss our result and test the CPU time. Sec. 6 contains our conclusions.

2. Derivation of Craig’s integral representation

Ref. [7] provides an approximation for the integral over the Gaussian standard normal distribution obtained by geometric considerations that is related to the cumulative distribution function via $P(t)=\Phi(t)-\Phi(-t)$, where $\Phi(t)$ is the Laplace function. The same considerations apply to the error function $\mathop{\rm erf}\nolimits(t)$ which is related to $P(t)$ via

[TABLE]

Translating the results of Ref. [7] to the error function, one obtains the approximation of order $p$ to be

[TABLE]

where the $N=2^{p}$ values $k_{p,n}$ ($n=1,2,\ldots,N$) are found in the intervals between $1/\cos(\pi(n-1)/(4N))$ and $1/\cos(\pi n/(4N))$. The way of selecting those values is extensively described in Ref. [7] where it is shown that

[TABLE]

for $k_{0,1}=1.116$, and with $14\approx 0.0033/0.00024$ times larger precision

[TABLE]

for $k_{1,1}=1.01$, $k_{1,2}=1.23345$. For the parameters taking the values $k_{p,n}=1/\cos(\pi n/(4N))$ of the upper limits of those intervals, it can be shown that the deviation is given by

[TABLE]

Given the values $k_{p,n}=1/\cos\phi(n)$ with $\phi(n)=\pi n/(4N)$, in the limit $N\to\infty$ the sum over $n$ in Eq. (2) can be replaced by an integral with measure $dn=(4N/\pi)d\phi(n)$ to obtain

[TABLE]

3. Power series expansion

The integral in Eq. (6) can be expanded into a power series in $t^{2}$,

[TABLE]

with

[TABLE]

where $y=\tan\phi$. The coefficients $c_{n}$ can be expressed by the hypergeometric function, $c_{n}={}_{2}F_{1}(1/2,1-n;3/2;-1)$, also known as Barnes’ extended hypergeometric function. On the other hand, we can derive a constraint for the explicit finite series expression for $c_{n}$ that renders the series in Eq. (7) to be convergent for all values of $t$. In order to be self-contained, intermediate steps to derive this constraint and to show the convergence are shown in the following. Necessary is Pascal’s rule

[TABLE]

and the sum over the rows of Pascal’s triangle,

[TABLE]

which can be shown by mathematical induction. The base case $n=0$ is obvious, as $(\begin{smallmatrix}0\\ 0\end{smallmatrix})=1=2^{0}$. For the induction step from $n$ to $n+1$ we write the first and last elements $(\begin{smallmatrix}n+1\\ 0\end{smallmatrix})=1$ and $(\begin{smallmatrix}n+1\\ n+1\end{smallmatrix})=1$ separately and use Pascal’s rule to obtain

[TABLE]

This proves Eq. (10). Returning to Eq. (8), one has $0\leq k\leq n-1$ and, therefore,

[TABLE]

For the result in Eq. (8) this means that

[TABLE]

i.e., the existence of a real number $c_{n}^{*}$ between $1/(2n-1)$ and $1$ such that $c_{n}=c_{n}^{*}2^{n-1}$. One has

[TABLE]

and because of $0\leq c_{n}^{*}\leq 1$ there is again a real number $c_{N}^{**}$ in the corresponding open interval so that

[TABLE]

As the latter is the power series expansion of $(2/\pi)e^{-2t^{2}}$ which is convergent for all values of $t$, also the original series is convergent and, therefore, $\mathop{\rm erf}\nolimits_{p}(t)^{2}$ with the limiting value shown in Eq. (7). A more compact form of the power series expansion is given by

[TABLE]

4. Approximations for the inverse error function

Based on the geometric approach from Ref. [7], in the following we describe how to find simple, handy formulas that, guided by higher and higher orders of the approximation (2) for the error function lead to more and more advanced approximation of the inverse error function. Starting point is the degree $p=0$, i.e., the approximation in Eq. (3). Inverting $E=\mathop{\rm erf}\nolimits_{0}(t)=(1-e^{-k_{0,1}^{2}t^{2}})^{1/2}$ leads to $t^{2}=-\ln(1-E^{2})/k_{0,1}^{2}$, and using the parameter $k_{0,1}=1.116$ from Eq. (3) gives

[TABLE]

For $0\leq E\leq 0.92$ the relative deviation $(T_{(0)}-t)/t$ from the exact value $t$ is less than $1.11\%$, for $0\leq E<1$ the deviation is less than $10\%$. Therefore, for $E>0.92$ a more precise formula has to be used. As such higher values for $E$ appear only in $8\%$ of the cases, this will not essentially influence the CPU time.

Continuing with $p=1$, we insert $T_{0}=\sqrt{-\ln(1-E^{2})}/k_{0,1}^{2}$ into Eq. (2) to obtain

[TABLE]

where $k_{1,1}=1.01$ and $k_{1,2}=1.23345$ are the same as for Eq. (4). Taking the derivative of Eq. (1) and approximating this by the difference quotient, one obtains

[TABLE]

leading to $t\approx T_{1}=T_{0}+\frac{1}{2}\sqrt{\pi}e^{T_{0}^{2}}(E-\mathop{\rm erf}\nolimits_{1}(T_{0}))$. In this case, for in the larger interval $0\leq E\leq 0.995$ the relative deviation $(T_{1}-t)/t$ is less than $0.1\%$. Using $\mathop{\rm erf}\nolimits_{2}(t)$ instead of $\mathop{\rm erf}\nolimits_{1}(t)$ and inserting $T_{1}$ instead of $T_{0}$ one obtains $T_{2}$ with a relative deviation of maximally $0.01\%$ for the same interval. The results are shown in Fig. 1.

The method to can be optimised by a method similar to the shooting method in boundary problems, giving dynamics to the calculation. Suppose that following one of the previous methods, for a particular argument $E$ we have found an approximation $t_{0}$ for the value of the inverse error function at this argument. Using $t_{1}=1.01t_{0}$, one can adjust the improved result

[TABLE]

by inserting $E=\mathop{\rm erf}\nolimits(t)$ and and calculating $A$ for $t=t_{1}$. In general, this procedure gives a vanishing deviation close to $E=0$. In this case and for $t_{0}=T_{1}$, in the interval $0\leq E\leq 0.7$ the maximal deviation is slightly larger than $10^{-6}=0.0001\%$ while up to $E=0.92$ the deviation is restricted to $10^{-5}=0.001\%$. A more general ansatz

[TABLE]

can be adjusted by inserting $E=\mathop{\rm erf}\nolimits(t)$ for $t=1.01t_{0}$ and $t=1.02t_{0}$, and the system of equations

[TABLE]

with $\Delta t=0.01t_{0}$, $\Delta E_{i}=\mathop{\rm erf}\nolimits(t_{i})-\mathop{\rm erf}\nolimits(t_{0})$ can be solved for $A$ and $B$ to obtain

[TABLE]

For $0\leq E\leq 0.70$ one obtains a relative deviation of $1.5\cdot 10^{-8}$, for $0\leq E\leq 0.92$ the maximal deviation is $5\cdot 10^{-7}$. Finally, the adjustment of

[TABLE]

leads to

[TABLE]

where $D=\Delta E_{1}\Delta E_{2}\Delta E_{3}(\Delta E_{1}-\Delta E_{2})(\Delta E_{1}-\Delta E_{3})(\Delta E_{2}-\Delta E_{3})$. For $0\leq E\leq 0.70$ the relative deviation is restricted to $5\cdot 10^{-10}$ while up to $E=0.92$ the maximal relative deviation is $4\cdot 10^{-8}$. The results for the deviations of $T_{(n)}$ ($n=1,2,3$) for linear, quadratic and cubic dynamical approximation are shown in Fig. 2.

5. Discussion

In order to test the feasibility and speed, we have coded our algorithm in the computer language C under Slackware 15.0 (linux 5.15.19) on an ordinary hp laptop with Intel® Core™2 Duo CPU P8600 @ 2.4GHz with 3MiB memory used. The dependence of the CPU time for the calculation is estimated by calculating the value $10^{6}$ times in sequence. The speed of the calculation of course does not depend on the value for $E$, as the precision is not optimised. This would have to be neccessary for a practical application. For an arbitrary starting value $E=0.8$ we perform this test, and the results are given in Table 1. An analysis of this table shows that a further step in the degree $p$ doubles the run time while the dynamics for increasing $n$ adds a constant value of approximately $0.06$ seconds to the result. Despite the fact that the increase of the dynamics needs the solution of a linear system of equations and the coding of the result, this endeavour is justified, as by using the dynamics one can increase the precision of the result without loosing calculational speed.

The results for the deviations in Figs. 1 and 2 are multiplied by increasing decimal powers in order to make the result comparable. This fact indicates that the convergence is improved in each of the steps for $p$ or $n$ at least by the corresponding inverse power. While the approximations in the statical approximations $n=0$ in Fig. 1 show both deviations close to $E=0$ and for higher values of $E$, the dynamical approximations in Fig. 2 show no deviation at $E=0$ and moderate deviations for higher values. On the other hand, the costs for an improvement step in either $p$ or $n$ is at most a factor of two higher CPU time. This means that the calculation and coding of expressions like Eqs. (4) is justified by the increase of precision gained. Given the goals for the precision, the user can decide to which degrees $p$ and $n$ the algorithm should be developed. In order to prove the precision, in Table 2 we show the convergence of our procedure for $p=2$ fixed and increasing values of $n$. The last column shows the CPU times for $10^{6}$ runs of the algorithm proposed in Ref. [12] with $N$ given in the last column of the table in Ref. [12], as coded in C.

6. Conclusions

In this paper we developed and described an approximative algorithm for the determination of the error function which is based on geometric considerations. Along the lines explained in this paper, the algorithm can be easily implemented and extended. We have shown that each improvement step gains an improvement of the precision of at least a factor of ten, at the cost of at most a factor of two more CPU time. As a bonus, we have given a geometric derivation of Craig’s integral representation of the error function and a converging power series expansion for this formula.

Acknowledgments

This research was funded by the European Regional Development Fund under Grant No. TK133.