逆错误功能需要代码 [英] Need code for Inverse Error Function

问题描述

有人知道我在哪里可以找到反错误函数"的代码吗? Freepascal/Delphi是更好的选择，但C/C ++也可以.

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.

TMath/DMath库没有它:(

The TMath/DMath library did not have it :(

推荐答案

这是erfinv()的实现.请注意，要使其正常运行，您还需要一个很好的erf()实现.

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;

如果您尚未实现erf()，则可以尝试将其从《数字食谱》转换为Pascal.不过，将精度提高一倍是不准确的.

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;

这篇关于逆错误功能需要代码的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！