Блог

Кто-нибудь знает, где я мог бы найти код для «Функции обратной ошибки»? Freepascal / Delphi был бы предпочтительнее, но C / C тоже подошел бы.

В библиотеке TMath / DMath его не было : (

Вот реализация erfinv() . Обратите внимание, что для ее эффективной работы вам также нужна хорошая реализация 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;
  

Если у вас нет реализации erf() , то вы можете попробовать эту, преобразованную в Pascal из числовых рецептов. Однако удвоенная точность не является точной.

 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;
  

В программах Pascal для ученых и инженеров есть гауссова функция ошибки (erf) и ее дополнение erfc =(1-errf), но не обратная функция ошибки. Очевидно, что вы не просто берете 1 / ErrF. Обратное означает, что x = erfinv(y) удовлетворяет y = erf (x).

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

Функция ошибки и ее дополнение показаны в этом списке.

Опять же, определение функции дополнения к ошибке является 1-ErrF , а не ErrF^-1 , но это должно подвести вас ближе:

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

Я нашел эту интересную реализацию (язык неизвестен, я предполагаю, что это matlab). возможно, это и его коэффициенты могут вам помочь:

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

Другой PDF здесь:http://people.maths.ox.ac.uk /~gilesm/files/gems_erfinv.pdf

Соответствующий фрагмент:

Таблица 1: Псевдокод для вычисления y = erfinv(x), где p1 (t) ..p6 (t) представляют полиномиальную функцию с 1-го по 6-й из 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
  

По-видимому, библиотечный код функционирует путем аппроксимации, это требует меньше работы. Я читал, что иногда приближения выполняются с точностью менее 6 знаков после запятой.

Код Fortran, который многие люди используют в качестве справочного, находится здесь, в нем цитируется «Рациональные приближения Чебышева для функции ошибки» У. Дж. Коди, математик. Сост., 1969, СТР. 631-638.:

Математика довольно сложная, но есть приличное приближение, описанное здесь (предупреждение: PDF), которое включает код Maple. К сожалению, это включает в себя шаг «решить для x», который может сделать его бесполезным для вас.

Похоже, что Boost имеет это как error_inv, поэтому посмотрите на код.

Я использовал это, которое, на мой взгляд, является достаточно точным и быстрым (обычно 2 итерации цикла), но, конечно, будьте осторожны. NormalX предполагает, что 0<=Q <=1, и, вероятно, даст глупые ответы, если это предположение не выполняется.

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

 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;
  

Вот что я получил от spe. Мне кажется, что они пытались ускорить это, перетащив все функции в один большой полином. Имейте в виду, что это код из эпохи 386

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