有谁知道我在哪里可以找到“逆误差函数”的代码？ 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()，那么你可以试试这个从数字食谱转换为帕斯卡的。但是双精度并不准确。
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;
【讨论】：
erf() 或者您正在使用我在上面从 NR 提供的那个？你有NR的副本吗？如果没有，你应该买一个。
Pascal Programs for Scientists and Engineers 具有高斯误差函数 (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
相关sn-p：
表 1：计算 y = erfinv(x) 的伪代码，其中 p1(t)..p6(t) 表示 t 的第 1 到第 6 多项式函数：
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 代码，许多人用作参考，是 here，它引用了 W. J. Cody, Math. Comp., 1969, PP. 631-638 的 "Rational Chebyshev approximations for the error function"。:
【讨论】：
Boost 似乎将其命名为 error_inv，因此请查看代码。
【讨论】：
我已经使用过这个，我认为它相当准确和快速（通常是循环的 2 次迭代），但当然要注意购买者。 NormalX 假设 0
/* 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};
【讨论】：