Required code for the inverse error function

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

-, , . 6 , .

Fortran, , , " WJ Cody, Math. Comp., 1969, PP. 631-638:.

+2

Warren P 12 '11 1:41

math , (: PDF), Maple. , "solve for x", .

+1

Tim Sylvester 12 '11 0:01

Boost, , error_inv, .

+1

lhf May 12, '11 at 0:05

source share

I used this, which, in my opinion, is reasonably accurate and fast (usually 2 iterations of a loop), but of course emptor warns. NormalX assumes that 0 <= Q <= 1 and is likely to give silly answers if this assumption is not fulfilled.

/* 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 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;
}

+1

dmuir May 12, '11 at 8:21

source share

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;

0

Kyle xy May 13, '11 at 12:08

source share

Here is what I have from spe. It seems to me that they tried to speed it up by dragging all the functions into one large polynomial. Keep in mind that 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};

0

Marco van de voort May 13, '11 at 20:04

source share

David heffernan · Accepted Answer · 2011-05-12T08:43:07+0000

This is where the implementation is done erfinv(). Note that in order for it to work well, you also need a good implementation 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 do not have an implementation erf(), you can try converting this to Pascal from Numerical Recipes. However, it is not so accurate to double the accuracy.

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;

Required code for the inverse error function

More articles: