ststat.pas

条码控件： 一维条码控件 二维条码控件 PDF417Barcode MaxiCodeBarcode

  Result := Combinations(N, K)*IntPower(p, K)*IntPower(1.0-p, N-K);
end;

function BinomDist(NumberS, Trials : Integer; ProbabilityS : Single;
  Cumulative : Boolean) : Single;
begin
  if (Trials < 0) or (NumberS < 0) or (NumberS > Trials) or
  (ProbabilityS < 0.0) or (ProbabilityS > 1.0) then
    RaiseStatError(stscStatBadParam);

  if (Cumulative) then
    Result := 1.0+BinomDistPmf(Trials, NumberS, ProbabilityS)-
      BetaDist(ProbabilityS, NumberS, Trials-NumberS+1, 0.0, 1.0)
  else
    Result := BinomDistPmf(Trials, NumberS, ProbabilityS);
end;

function CritBinom(Trials : Integer; ProbabilityS, Alpha : Single) : Integer;
var
  s : Integer;
  B : Double;
begin
  if (Trials < 0) or (ProbabilityS < 0.0) or (ProbabilityS > 1.0) or
  (Alpha < 0.0) or (Alpha > 1.0) then
    RaiseStatError(stscStatBadParam);

  B := 0.0;
  for s := 0 to Trials do begin
    B := B+BinomDistPmf(Trials, s, ProbabilityS);
    if (B >= Alpha) then begin
      Result := s;
      exit;
    end;
  end;
  {in case of roundoff problems}
  Result := Trials;
end;

function GammSer(a, x : Single) : Single;
  {-Returns the series computation for GammP}
const
  MAXIT = 100;
  EPS = 3.0E-7;
var
  N : Integer;
  sum, del, ap : Double;
begin
  Result := 0.0;
  if (x > 0.0) then begin
    {x shouldn't be < 0.0, tested by caller}
    ap := a;
    sum := 1.0/a;
    del := sum;
    for N := 1 to MAXIT do begin
      ap := ap+1;
      del := del*x/ap;
      sum := sum+del;
      if (abs(del) < abs(sum)*EPS) then begin
        Result := sum*exp(-X+a*ln(X)-GammaLn(a));
        exit;
      end;
    end;
    RaiseStatError(stscStatNoConverge);
  end;
end;

function GammCf(a, x : Single) : Single;
  {-Returns the continued fraction computation for GammP}
const
  MAXIT = 100;
  EPS = 3.0e-7;
  FPMIN = 1.0e-30;
var
  i : Integer;
  an, b, c, d, del, h : Double;
begin
  b := x+1.0-a;
  c := 1.0/FPMIN;
  d := 1.0/b;
  h := d;
  for i := 1 to MAXIT do begin
    an := -i*(i-a);
    b := b+2.0;
    d := an*d+b;
    if (abs(d) < FPMIN) then
      d := FPMIN;
    c := b+an/c;
    if (abs(c) < FPMIN) then
      c := FPMIN;
    d := 1.0/d;
    del := d*c;
    h := h*del;
    if (abs(del-1.0) < EPS) then
      break;
    if i = MAXIT then
      RaiseStatError(stscStatNoConverge);
  end;
  Result := h*exp(-x+a*ln(x)-GammaLn(a));
end;

function GammP(a, x : Single) : Single;
  {-Returns the incomplete gamma function P(a, x)}
begin
  if (x < 0.0) or (a <= 0.0) then
    RaiseStatError(stscStatBadParam);
  if (x < a+1.0) then
    {use the series representation}
    Result := GammSer(a, x)
  else
    {use the continued fraction representation}
    Result := 1.0-GammCf(a, x);
end;

function ChiDist(X : Single; DegreesFreedom : Integer) : Single;
begin
  if (DegreesFreedom < 1) or (X < 0.0) then
    RaiseStatError(stscStatBadParam);
  Result := 1.0-GammP(DegreesFreedom/2.0, X/2.0);
end;

function ChiInv(Probability : Single; DegreesFreedom : Integer) : Single;
const
  MAXIT = 100;
  UNUSED = -1.11e30;
  XACC = 3e-7;
  FACTOR = 1.6;
var
  j : Integer;
  ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew, dsign : Double;
begin
  if (Probability < 0.0) or (Probability > 1.0) or (DegreesFreedom < 1) then
    RaiseStatError(stscStatBadParam);

  if (Probability = 0.0) then
    Result := 0.0
  else begin
    {bracket the interval}
    xl := 0.0;
    xh := 100.0;
    fl := ChiDist(xl, DegreesFreedom)-Probability;
    fh := ChiDist(xh, DegreesFreedom)-Probability;
    for j := 1 to MAXIT do begin
      if (fl*fh < 0.0) then
        {bracketed the root}
        break;
      if (abs(fl) < abs(fh)) then begin
        xl := xl+FACTOR*(xl-xh);
        fl := ChiDist(xl, DegreesFreedom)-Probability;
      end else begin
        xh := xh+FACTOR*(xh-xl);
        fh := ChiDist(xh, DegreesFreedom)-Probability;
      end;
    end;
    if (fl*fh >= 0.0) then
      {couldn't bracket it, means Probability too close to 1.0}
      RaiseStatError(stscStatNoConverge);

    {Ridder's method of finding the root of ChiDist = Probability}
    ans := UNUSED;

    for j := 1 to MAXIT do begin
      xm := 0.5*(xl+xh);
      fm := ChiDist(xm, DegreesFreedom)-Probability;
      s := sqrt(fm*fm-fl*fh);
      if (s = 0.0) then begin
        Result := ans;
        exit;
      end;
      if (fl >= fh) then
        dsign := 1.0
      else
        dsign := -1.0;
      xnew := xm+(xm-xl)*(dsign*fm/s);
      if (abs(xnew-ans) <= XACC) then begin
        Result := ans;
        exit;
      end;
      ans := xnew;

      fnew := ChiDist(ans, DegreesFreedom)-Probability;
      if (fnew = 0.0) then begin
        Result := ans;
        exit;
      end;

      {keep root bracketed on next iteration}
      if (Sign(fm, fnew) <> fm) then begin
        xl := xm;
        fl := fm;
        xh := ans;
        fh := fnew;
      end else if (Sign(fl, fnew) <> fl) then begin
        xh := ans;
        fh := fnew;
      end else if (Sign(fh, fnew) <> fh) then begin
        xl := ans;
        fl := fnew;
      end else
        {shouldn't get here}
        RaiseStatError(stscStatNoConverge);

      if (abs(xh-xl) <= XACC) then begin
        Result := ans;
        exit;
      end;
    end;
    Result := 0.0; {avoid compiler warning}
    RaiseStatError(stscStatNoConverge);
  end;
end;

function ExponDist(X, Lambda : Single; Cumulative : Boolean) : Single;
begin
  if (X < 0.0) or (Lambda <= 0.0) then
    RaiseStatError(stscStatBadParam);

  if (Cumulative) then
    Result := 1.0-exp(-Lambda*X)
  else
    Result := Lambda*exp(-Lambda*X);
end;

function FDist(X : Single;
  DegreesFreedom1, DegreesFreedom2 : Integer) : Single;
begin
  if (X < 0) or (DegreesFreedom1 < 1) or (DegreesFreedom2 < 1) then
    RaiseStatError(stscStatBadParam);

  Result := BetaDist(DegreesFreedom2/(DegreesFreedom2+DegreesFreedom1*X),
    DegreesFreedom2/2.0, DegreesFreedom1/2.0, 0.0, 1.0);
end;

function FInv(Probability : Single;
  DegreesFreedom1, DegreesFreedom2 : Integer) : Single;
const
  MAXIT = 100;
  UNUSED = -1.11e30;
  XACC = 3e-7;
  FACTOR = 1.6;
var
  j : Integer;
  ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew, dsign : Double;
begin
  if (Probability < 0.0) or (Probability > 1.0) or
  (DegreesFreedom1 < 1) or (DegreesFreedom2 < 1) then
    RaiseStatError(stscStatBadParam);

  if (Probability = 1.0) then
    Result := 0.0
  else begin
    {bracket the interval}
    xl := 0.0;
    xh := 100.0;
    fl := FDist(xl, DegreesFreedom1, DegreesFreedom2)-Probability;
    fh := FDist(xh, DegreesFreedom1, DegreesFreedom2)-Probability;
    for j := 1 to MAXIT do begin
      if (fl*fh < 0.0) then
        {bracketed the root}
        break;
      if (abs(fl) < abs(fh)) then begin
        xl := xl+FACTOR*(xl-xh);
        fl := FDist(xl, DegreesFreedom1, DegreesFreedom2)-Probability;
      end else begin
        xh := xh+FACTOR*(xh-xl);
        fh := FDist(xh, DegreesFreedom1, DegreesFreedom2)-Probability;
      end;
    end;
    if (fl*fh >= 0.0) then
      {couldn't bracket it, means Probability too close to 0.0}
      RaiseStatError(stscStatNoConverge);

    {Ridder's method of finding the root of FDist = Probability}
    ans := UNUSED;

    for j := 1 to MAXIT do begin
      xm := 0.5*(xl+xh);
      fm := FDist(xm, DegreesFreedom1, DegreesFreedom2)-Probability;
      s := sqrt(fm*fm-fl*fh);
      if (s = 0.0) then begin
        Result := ans;
        exit;
      end;
      if (fl >= fh) then
        dsign := 1.0
      else
        dsign := -1.0;
      xnew := xm+(xm-xl)*(dsign*fm/s);
      if (abs(xnew-ans) <= XACC) then begin
        Result := ans;
        exit;
      end;
      ans := xnew;

      fnew := FDist(ans, DegreesFreedom1, DegreesFreedom2)-Probability;
      if (fnew = 0.0) then begin
        Result := ans;
        exit;
      end;

      {keep root bracketed on next iteration}
      if (Sign(fm, fnew) <> fm) then begin
        xl := xm;
        fl := fm;
        xh := ans;
        fh := fnew;
      end else if (Sign(fl, fnew) <> fl) then begin
        xh := ans;
        fh := fnew;
      end else if (Sign(fh, fnew) <> fh) then begin
        xl := ans;
        fl := fnew;
      end else
        {shouldn't get here}
        RaiseStatError(stscStatNoConverge);

      if (abs(xh-xl) <= XACC) then begin
        Result := ans;
        exit;
      end;
    end;
    Result := 0.0; {avoid compiler warning}
    RaiseStatError(stscStatNoConverge);
  end;
end;

function LogNormDist(X, Mean, StandardDev : Single) : Single;
begin
  if (X <= 0.0) or (StandardDev <= 0) then
    RaiseStatError(stscStatBadParam);
  Result := NormSDist((ln(X)-Mean)/StandardDev);
end;

function LogInv(Probability,  Mean, StandardDev : Single) : Single;
begin
  if (Probability < 0.0) or (Probability > 1.0) or (StandardDev <= 0) then
    RaiseStatError(stscStatBadParam);
  Result := exp(Mean+StandardDev*NormSInv(Probability));
end;

function NormDist(X, Mean, StandardDev : Single;
  Cumulative : Boolean) : Single;
var
  Z : Extended;
begin
  if (StandardDev <= 0) then
    RaiseStatError(stscStatBadParam);
  Z := (X-Mean)/StandardDev;
  if (Cumulative) then
    Result := NormSDist(Z)
  else
    Result := exp(-Z*Z/2.0)/(StandardDev*sqrt2pi);
end;

function NormInv(Probability, Mean, StandardDev : Single) : Single;
begin
  if (Probability < 0.0) or (Probability > 1.0) or (StandardDev <= 0) then
    RaiseStatError(stscStatBadParam);
  Result := Mean+StandardDev*NormSInv(Probability);
end;

function Erfc(X : Single) : Single;
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
    Result := ans
  else
    Result := 2.0-ans;
end;

function NormSDist(Z : Single) : Single;
const
  sqrt2 = 1.41421356237310;
begin
  Result := 1.0-0.5*Erfc(Z/sqrt2);
end;

function NormSInv(Probability : Single) : Single;
const
  MAXIT = 100;
  UNUSED = -1.11e30;
  XACC = 3e-7;
  FACTOR = 1.6;
var
  j : Integer;
  ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew, dsign : Double;
begin
  if (Probability < 0.0) or (Probability > 1.0) then
    RaiseStatError(stscStatBadParam);
  Result := 0.0;

  {bracket the interval}
  xl := -2.0;
  xh := +2.0;
  fl := NormSDist(xl)-Probability;
  fh := NormSDist(xh)-Probability;
  for j := 1 to MAXIT do begin
    if (fl*fh < 0.0) then
      {bracketed the root}
      break;
    if (abs(fl) < abs(fh)) then begin
      xl := xl+FACTOR*(xl-xh);
      fl := NormSDist(xl)-Probability;
    end else begin
      xh := xh+FACTOR*(xh-xl);
      fh := NormSDist(xh)-Probability;
    end;
  end;
  if (fl*fh >= 0.0) then
    {couldn't bracket it, means Probability too close to 0.0 or 1.0}
    RaiseStatError(stscStatNoConverge);

  {Ridder's method of finding the root of NormSDist = Probability}
  ans := UNUSED;

  for j := 1 to MAXIT do begin
    xm := 0.5*(xl+xh);
    fm := NormSDist(xm)-Probability;
    s := sqrt(fm*fm-fl*fh);
    if (s = 0.0) then begin
      Result := ans;
      exit;
    end;
    if (fl >= fh) then
      dsign := 1.0
    else
      dsign := -1.0;
    xnew := xm+(x