Need code for Inverse Error Function

Question

Does anyone know where I could find code for the "Inverse Error Function?" Freepascal/Delphi would be preferable but C/C++ would be fine too.

The TMath/DMath library did not have it :(

Answer 1

Here's an implementation of erfinv(). Note that for it to work well, you also need a good implementation of erf().

function erfinv(const y: Double): Double;

//rational approx coefficients
const
  a: array [0..3] of Double = ( 0.886226899, -1.645349621,  0.914624893, -0.140543331);
  b: array [0..3] of Double = (-2.118377725,  1.442710462, -0.329097515,  0.012229801);
  c: array [0..3] of Double = (-1.970840454, -1.624906493,  3.429567803,  1.641345311);
  d: array [0..1] of Double = ( 3.543889200,  1.637067800);

const
  y0 = 0.7;

var
  x, z: Double;

begin
  if not InRange(y, -1.0, 1.0) then begin
    raise EInvalidArgument.Create('erfinv(y) argument out of range');
  end;

  if abs(y)=1.0 then begin
    x := -y*Ln(0.0);
  end else if y<-y0 then begin
    z := sqrt(-Ln((1.0+y)/2.0));
    x := -(((c[3]*z+c[2])*z+c[1])*z+c[0])/((d[1]*z+d[0])*z+1.0);
  end else begin
    if y<y0 then begin
      z := y*y;
      x := y*(((a[3]*z+a[2])*z+a[1])*z+a[0])/((((b[3]*z+b[3])*z+b[1])*z+b[0])*z+1.0);
    end else begin
      z := sqrt(-Ln((1.0-y)/2.0));
      x := (((c[3]*z+c[2])*z+c[1])*z+c[0])/((d[1]*z+d[0])*z+1.0);
    end;
    //polish x to full accuracy
    x := x - (erf(x) - y) / (2.0/sqrt(pi) * exp(-x*x));
    x := x - (erf(x) - y) / (2.0/sqrt(pi) * exp(-x*x));
  end;

  Result := x;
end;

If you haven't got an implementation of erf() then you can try this one converted to Pascal from Numerical Recipes. It's not accurate to double precision though.

function erfc(const x: Double): Double;
var
  t,z,ans: Double;
begin
  z := abs(x);
  t := 1.0/(1.0+0.5*z);
  ans := t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
    t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
    t*(-0.82215223+t*0.17087277)))))))));
  if x>=0.0 then begin
    Result := ans;
  end else begin
    Result := 2.0-ans;
  end;
end;

function erf(const x: Double): Double;
begin
  Result := 1.0-erfc(x);
end;

Answer 2

Pascal Programs for Scientists and Engineers has the gaussian Error function (erf) and its complement erfc=(1-errf), but not the Inverse of the Error function. Obviously, you don't just take 1/ErrF. The inverse means x = erfinv(y) satisfies y = erf(x).

http://infohost.nmt.edu/~armiller/pascal.htm

Error function and its complement, are shown in this listing.

Again, the definition of Error Function Complement is 1-ErrF, not ErrF^-1, but this has got to be getting you close:

http://infohost.nmt.edu/~es421/pascal/list11-3.pas

I found this interesting implementation (language unknown, I'm guessing it's matlab). maybe it and its coefficients can help you:

http://w3eos.whoi.edu/12.747/mfiles/lect07/erfinv.m

Another PDF here: http://people.maths.ox.ac.uk/~gilesm/files/gems_erfinv.pdf

Relevant snippet:

Table 1: Pseudo-code to compute y = erfinv(x) , with p1(t)..p6(t) representing a 1st through 6th polynomial function of t :

a = |x|        
if a > 0.9375 then
t = sqrt( log(a) )
y = p1(t) / p2(t)
else if a > 0.75 then
y = p3(a) / p4(a)
else
y = p5(a) / p6(a)
end if
if x < 0 then
y = −y
end if

Apparently the library code functions by approximation, it's less work. Sometimes the approximations are to less than 6 decimal places accuracy, I read.

Fortran code that many people use for a reference, is here, it cites "Rational Chebyshev approximations for the error function" by W. J. Cody, Math. Comp., 1969, PP. 631-638.:

Answer 3

The math is pretty complex, but there's a decent approximation described here (warning: PDF) that includes Maple code. Unfortunately it involves a "solve for x" step that might make it useless to you.

Answer 4

Boost seems to have it as error_inv so look at the code.

Answer 5

I've used this, which I believe is reasonably accurate and quick (usually 2 iterations of the loop), but of course caveat emptor. NormalX assumes that 0<=Q<=1, and would likely give silly answers if that assumption doesn't hold.

/* return P(N>X) for normal N */
double  NormalQ( double x)
{   return 0.5*erfc( x/sqrt(2.0));
}

#define NORX_C0   2.8650422353e+00
#define NORX_C1   3.3271545598e+00
#define NORX_C2   2.7147548996e-01
#define NORX_D1   2.8716448975e+00
#define NORX_D2   1.1690926940e+00
#define NORX_D3   4.7994444496e-02
/* return X such that P(N>X) = Q for normal N */
double  NormalX( double Q)  
{
double  eps = 1e-12;
int signum = Q < 0.5;
double  QF = signum ? Q : (1.0-Q);
double  T = sqrt( -2.0*log(QF));
double  X = T - ((NORX_C2*T + NORX_C1)*T + NORX_C0)
                    /(((NORX_D3*T + NORX_D2)*T + NORX_D1)*T + 1.0);
double  SPI2 = sqrt( 2.0 * M_PI);
int i;
    /* newton's method */
    for( i=0; i<10; ++i)
    {
    double  dX  = (NormalQ(X) - QF)*exp(0.5*X*X)*SPI2;
            X += dX;
            if ( fabs( dX) < eps)   
            {   break;
            }
    }
    return signum ? X : -X;
}

Answer 6

function erf(const x: extended): extended;
var
  n: integer;
  z: extended;
begin
  Result := x;
  z := x;
  n := 0;

  repeat
    inc(n);
    z := -z * x * x * (2 * n - 1) / ((2 * n + 1) * n);
    Result := Result + z;
  until abs(z) < 1E-20;

  Result := Result * 2 / sqrt(pi);
end;

function erfinv(const x: extended): extended;
var
  n: integer;
  z: extended;
begin
  Result := 0;
  n := 0;

  repeat
    inc(n);
    z := (erf(Result) - x) * sqrt(pi) / (2 * exp(-Result * Result));
    Result := Result - z;
  until (n = 100) or (abs(z) < 1E-20);

  if abs(z) < 1E-20 then
    n := -20
  else
    n := Floor(Log10(abs(z))) + 1;

  Result := RoundTo(Result, n);
end;

Answer 7

Here is what I have from spe. To me it looks like they tried to speed it up by dragging all functions into one big polynomal. Keep in mind this is code from the 386 era

// Extract from package Numlib in the Free Pascal sources (http://www.freepascal.org)
// LGPL-with-static-linking-exception, see original source for exact license.
// Note this is usually compiled in TP mode, not in Delphi mode.

Const highestElement = 20000000;

Type ArbFloat = double;  // can be extended too.
     ArbInt   = Longint;
     arfloat0   = array[0..highestelement] of ArbFloat;

function spesgn(x: ArbFloat): ArbInt;

begin
  if x<0
  then
    spesgn:=-1
  else
    if x=0
    then
      spesgn:=0
    else
      spesgn:=1
end; {spesgn}

function spepol(x: ArbFloat; var a: ArbFloat; n: ArbInt): ArbFloat;
var   pa : ^arfloat0;
       i : ArbInt;
    polx : ArbFloat;
begin
  pa:=@a;
  polx:=0;
  for i:=n downto 0 do
    polx:=polx*x+pa^[i];
  spepol:=polx
end {spepol};

function speerf(x : ArbFloat) : ArbFloat;
const

        xup = 6.25;
     sqrtpi = 1.7724538509055160;
     c : array[1..18] of ArbFloat =
     ( 1.9449071068178803e0,  4.20186582324414e-2, -1.86866103976769e-2,
       5.1281061839107e-3,   -1.0683107461726e-3,   1.744737872522e-4,
      -2.15642065714e-5,      1.7282657974e-6,     -2.00479241e-8,
      -1.64782105e-8,         2.0008475e-9,         2.57716e-11,
      -3.06343e-11,           1.9158e-12,           3.703e-13,
      -5.43e-14,             -4.0e-15,              1.2e-15);

     d : array[1..17] of ArbFloat =
     ( 1.4831105640848036e0, -3.010710733865950e-1, 6.89948306898316e-2,
      -1.39162712647222e-2,   2.4207995224335e-3,  -3.658639685849e-4,
       4.86209844323e-5,     -5.7492565580e-6,      6.113243578e-7,
      -5.89910153e-8,         5.2070091e-9,        -4.232976e-10,
       3.18811e-11,          -2.2361e-12,           1.467e-13,
      -9.0e-15,               5.0e-16);

  var t, s, s1, s2, x2: ArbFloat;
         bovc, bovd, j: ArbInt;
begin
  bovc:=sizeof(c) div sizeof(ArbFloat);
  bovd:=sizeof(d) div sizeof(ArbFloat);
  t:=abs(x);
  if t <= 2
  then
    begin
      x2:=sqr(x)-2;
      s1:=d[bovd]; s2:=0; j:=bovd-1;
      s:=x2*s1-s2+d[j];
      while j > 1 do
        begin
          s2:=s1; s1:=s; j:=j-1;
          s:=x2*s1-s2+d[j]
        end;
      speerf:=(s-s2)*x/2
    end
  else
    if t < xup
    then
      begin
        x2:=2-20/(t+3);
        s1:=c[bovc]; s2:=0; j:=bovc-1;
        s:=x2*s1-s2+c[j];
        while j > 1 do
          begin
            s2:=s1; s1:=s; j:=j-1;
            s:=x2*s1-s2+c[j]
          end;
        x2:=((s-s2)/(2*t))*exp(-sqr(x))/sqrtpi;
        speerf:=(1-x2)*spesgn(x)
      end
    else
      speerf:=spesgn(x)
end;  {speerf}

function spemax(a, b: ArbFloat): ArbFloat;
begin
  if a>b
  then
    spemax:=a
  else
    spemax:=b
end; {spemax}

function speefc(x : ArbFloat) : ArbFloat;
const

   xlow = -6.25;
  xhigh = 27.28;
      c : array[0..22] of ArbFloat =
      ( 1.455897212750385e-1, -2.734219314954260e-1,
        2.260080669166197e-1, -1.635718955239687e-1,
        1.026043120322792e-1, -5.480232669380236e-2,
        2.414322397093253e-2, -8.220621168415435e-3,
        1.802962431316418e-3, -2.553523453642242e-5,
       -1.524627476123466e-4,  4.789838226695987e-5,
        3.584014089915968e-6, -6.182369348098529e-6,
        7.478317101785790e-7,  6.575825478226343e-7,
       -1.822565715362025e-7, -7.043998994397452e-8,
        3.026547320064576e-8,  7.532536116142436e-9,
       -4.066088879757269e-9, -5.718639670776992e-10,
        3.328130055126039e-10);

  var t, s : ArbFloat;
begin
  if x <= xlow
  then
    speefc:=2
  else
  if x >= xhigh
  then
    speefc:=0
  else
    begin
      t:=1-7.5/(abs(x)+3.75);
      s:=exp(-x*x)*spepol(t, c[0], sizeof(c) div sizeof(ArbFloat) - 1);
      if x < 0
      then
        speefc:=2-s
      else
        speefc:=s
    end
end {speefc};