如何提高标准正态分布的计算的精度

Question

通常都是查标准正态分布表的，这次需要写一个算法，使得数据更准确。但是通过验算，发现似乎还不如查表来得准确，和MiniTab计算的值也有些差距：（
    不知道Delphi有没有相关的函数，以及如何提高计算精度。下面是我用的算法。请各位不吝赐教！
＝＝＝＝＝
unit UnitMath;interfacetype
    TMath=class
        function factorial(n:Integer):Int64;
        function power(x:Double;n:Integer):Extended;
        function getNormal(x:Double;precision:Double):Double;
    end;implementation//计算阶乘
function TMath.factorial(n:Integer):Int64;
var
        i:Integer;begin
        Result:=1;        if(n=0) then
                exit;        for i:=1 to n do
                Result:=Result*i;
end;//计算乘方
function TMath.power(x:Double;n:Integer):Extended;
var
        i:Integer;
begin
        Result:=1;
        if n=0 then exit;        for i:=1 to n do
                Result:=Result*x;
end;//计算正态分布
function TMath.getNormal(x:Double;precision:Double):Double;
var
        i:Integer;
        iD:Integer;
        f:Double;
        fCommon:Double;
begin
        f:=x/sqrt(2);
        if(f<0) then f:=-f;
        Result:=f;
        iD:=1;
        i:=0;
        while(true) do begin
                i:=i+1;
                iD:=iD*(-1);
                fCommon:=self.power(f,2*i+1) / (2*i+1) / self.factorial(i);
                Result:=Result + iD * fCommon;
                if(fCommon/sqrt(pi)<precision) then break;
        end;
        Result:=0.5+Result/sqrt(pi);
        if(Result>1) then Result:=1;
        if(x<0) then Result:=1-Result;
end;end.

Answer 1

http://fundementals.sourceforge.net/cStatistics.htmlfunction  CummNormal01(const X: Extended): Extended算出来的结果却不对！

Answer 2

unit UnitMath;interfaceFUNCTION gammp(a,x: Extended): Extended;
FUNCTION gammq(a,x: Extended): Extended;
PROCEDURE gser(a,x: Extended; VAR gamser,gln: Extended);
PROCEDURE gcf(a,x: Extended; VAR gammcf,gln: Extended);
FUNCTION gammln(xx: Extended): Extended;
FUNCTION erf(x: Extended): Extended;
FUNCTION erfc(x: Extended): Extended;
function CummNormal01(const X: Extended): Extended;implementation
const Sqrt2     = 1.414213562373095048801688724209698;
function CummNormal01(const X: Extended): Extended;
begin
  CummNormal01 := erfc(-X / Sqrt2) / 2.0;
end;FUNCTION erf(x: Extended): Extended;
BEGIN
   IF (x < 0.0) THEN BEGIN
      erf := -gammp(0.5,sqr(x))
   END ELSE BEGIN
      erf := gammp(0.5,sqr(x))
   END
END;
FUNCTION erfc(x: Extended): Extended;
BEGIN
   IF (x < 0.0) THEN BEGIN
      erfc := 1.0+gammp(0.5,sqr(x))
   END ELSE BEGIN
      erfc := gammq(0.5,sqr(x))
   END
END;FUNCTION gammq(a,x: Extended): Extended;
VAR
   gamser,gln: Extended;
BEGIN
   IF ((x < 0.0) OR (a <= 0.0)) THEN BEGIN
      writeln('pause in GAMMQ - invalid arguments'); readln
   END;
   IF (x < a+1.0) THEN BEGIN
      gser(a,x,gamser,gln);
      gammq := 1.0-gamser
   END ELSE BEGIN
      gcf(a,x,gamser,gln);
      gammq := gamser
   END
END;
PROCEDURE gser(a,x: Extended; VAR gamser,gln: Extended);
LABEL 1;
CONST
   itmax=100;
   eps=3.0e-7;
VAR
   n: integer;
   sum,del,ap: Extended;
BEGIN
   gln := gammln(a);
   IF (x <= 0.0) THEN BEGIN
      IF (x < 0.0) THEN BEGIN
         writeln('pause in GSER - x less than 0'); readln
      END;
      gamser := 0.0
   END ELSE BEGIN
      ap := a;
      sum := 1.0/a;
      del := sum;
      FOR n := 1 TO itmax DO BEGIN
         ap := ap+1.0;
         del := del*x/ap;
         sum := sum+del;
         IF (abs(del) < abs(sum)*eps) THEN GOTO 1
      END;
      writeln('pause in GSER - a too large, itmax too small'); readln;
1:      gamser := sum*exp(-x+a*ln(x)-gln)
   END
END;FUNCTION gammln(xx: Extended): Extended;
CONST
   stp = 2.50662827465;
   half = 0.5;
   one = 1.0;
   fpf = 5.5;
VAR
   x,tmp,ser: double;
   j: integer;
   cof: ARRAY [1..6] OF double;
BEGIN
   cof[1] := 76.18009173;
   cof[2] := -86.50532033;
   cof[3] := 24.01409822;
   cof[4] := -1.231739516;
   cof[5] := 0.120858003e-2;
   cof[6] := -0.536382e-5;
   x := xx-one;
   tmp := x+fpf;
   tmp := (x+half)*ln(tmp)-tmp;
   ser := one;
   FOR j := 1 TO 6 DO BEGIN
      x := x+one;
      ser := ser+cof[j]/x
   END;
   gammln := tmp+ln(stp*ser)
END;PROCEDURE gcf(a,x: Extended; VAR gammcf,gln: Extended);
LABEL 1;
CONST
   itmax=100;
   eps=3.0e-7;
VAR
   n: integer;
   gold,g,fac,b1,b0,anf,ana,an,a1,a0: Extended;
BEGIN
   gln := gammln(a);
   gold := 0.0;
   a0 := 1.0;
   a1 := x;
   b0 := 0.0;
   b1 := 1.0;
   fac := 1.0;
   FOR n := 1 TO itmax DO BEGIN
      an := 1.0*n;
      ana := an-a;
      a0 := (a1+a0*ana)*fac;
      b0 := (b1+b0*ana)*fac;
      anf := an*fac;
      a1 := x*a0+anf*a1;
      b1 := x*b0+anf*b1;
      IF (a1 <> 0.0) THEN BEGIN
         fac := 1.0/a1;
         g := b1*fac;
         IF (abs((g-gold)/g) < eps) THEN GOTO 1;
         gold := g
      END
   END;
   //writeln('pause in GCF - a too large, itmax too small'); readln;
1:   gammcf := exp(-x+a*ln(x)-gln)*g
END;FUNCTION gammp(a,x: Extended): Extended;
VAR
   gammcf,gln: Extended;
BEGIN
   IF ((x < 0.0) OR (a <= 0.0)) THEN BEGIN
      writeln('pause in GAMMP - invalid arguments'); readln
   END;
   IF (x < (a+1.0)) THEN BEGIN
      gser(a,x,gammcf,gln);
      gammp := gammcf
   END ELSE BEGIN
      gcf(a,x,gammcf,gln);
      gammp := 1.0-gammcf
   END
END;end.