ststat.pas

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

{$IFDEF Debug}
function SmallestSort(const Data: array of double; K : Integer) : Double;
var
  Size : Cardinal;
  NData : Integer;
  SD : PDoubleArray;
begin
  NData := High(Data)+1;
  if (NData <= 0) then
    RaiseStatError(stscStatBadCount);
  if (K <= 0) or (K > NData) then
    RaiseStatError(stscStatBadParam);

  {copy and sort array}
  Size := CopyAndSort(Data, NData, SD);
  try
    {K=1 returns smallest value, K=NData returns largest}
    Result := SD^[K-1];
  finally
    freemem(SD, Size);
  end;
end;
{$ENDIF}

function TrimMean(const Data: array of Double; Percent : Double) : Double;
begin
  Result := TrimMean16(Data, High(Data)+1, Percent);
end;

function TrimMean16(const Data; NData : Integer; Percent : Double) : Double;
var
  ntrim : Integer;
  Size : Cardinal;
  SD : PDoubleArray;
begin
  if (NData <= 0) then
    RaiseStatError(stscStatBadCount);
  if (Percent < 0.0) or (Percent > 1.0) then
    RaiseStatError(stscStatBadParam);

  {compute total number of trimmed points, rounding down to an even number}
  ntrim := trunc(Percent*NData);
  if odd(ntrim) then
    dec(ntrim);

  {take the easy way out when possible}
  if (ntrim = 0) then begin
    Result := Mean(Data, NData);
    exit;
  end;

  {copy and sort array}
  Size := CopyAndSort(Data, NData, SD);
  try
    {return Mean of remaining data points}
    Result := Mean(SD^[ntrim shr 1], NData-ntrim);
  finally
    freemem(SD, Size);
  end;
end;

{------------------------------------------------------------------------}

procedure LinEst(const KnownY: array of Double;
  const KnownX: array of Double; var LF : TStLinEst; ErrorStats : Boolean);
begin
  if (High(KnownY) <> High(KnownX)) then
    RaiseStatError(stscStatBadCount);
  LinEst16(KnownY, KnownX, High(KnownY)+1, LF, ErrorStats);
end;

procedure LinEst16(const KnownY; const KnownX; NData : Integer;
  var LF : TStLinEst; ErrorStats : Boolean);
var
  i : Integer;
  sx, sy, xmean, ymean, sxx, sxy, syy, x, y : Extended;
begin
  if (NData <= 2) then
    RaiseStatError(stscStatBadCount);

  {compute basic sums}
  sx := 0.0;
  sy := 0.0;
  sxx := 0.0;
  sxy := 0.0;
  syy := 0.0;
  for i := 0 to NData-1 do begin
    x := TDoubleArray(KnownX)[i];
    y := TDoubleArray(KnownY)[i];
    sx := sx+x;
    sy := sy+y;
    sxx := sxx+x*x;
    syy := syy+y*y;
    sxy := sxy+x*y;
  end;
  xmean := sx/NData;
  ymean := sy/NData;
  sxx := sxx-NData*xmean*xmean;
  syy := syy-NData*ymean*ymean;
  sxy := sxy-NData*xmean*ymean;

  {check for zero variance}
  if (sxx <= 0.0) or (syy <= 0.0) then
    RaiseStatError(stscStatBadData);

  {initialize returned parameters}
  fillchar(LF, sizeof(LF), 0);

  {regression coefficients}
  LF.B1 := sxy/sxx;
  LF.B0 := ymean-LF.B1*xmean;

  {error statistics}
  if (ErrorStats) then begin
    LF.ssr := LF.B1*sxy;
    LF.sse := syy-LF.ssr;
    LF.R2 := LF.ssr/syy;
    LF.df := NData-2;
    LF.sigma := sqrt(LF.sse/LF.df);
    if LF.sse = 0.0 then
      {pick an arbitrarily large number for perfect fit}
      LF.F0 := 1.7e+308
    else
      LF.F0 := (LF.ssr*LF.df)/LF.sse;
    LF.seB1 := LF.sigma/sqrt(sxx);
    LF.seB0 := LF.sigma*sqrt((1.0/NData)+(xmean*xmean/sxx));
  end;
end;

procedure LogEst(const KnownY: array of Double;
  const KnownX: array of Double; var LF : TStLinEst; ErrorStats : Boolean);
begin
  if (High(KnownY) <> High(KnownX)) then
    RaiseStatError(stscStatBadCount);
  LogEst16(KnownY, KnownX, High(KnownY)+1, LF, ErrorStats);
end;

procedure LogEst16(const KnownY; const KnownX; NData : Integer;
  var LF : TStLinEst; ErrorStats : Boolean);
var
  i : Integer;
  Size : LongInt;
  lny : PDoubleArray;
begin
  if (NData <= 2) then
    RaiseStatError(stscStatBadCount);

  {allocate array for the log-transformed data}
  Size := LongInt(NData)*sizeof(Double);
  {f (Size > MaxBlockSize) then}
  { RaiseStatError(stscStatBadCount);}
  getmem(lny, Size);
  try
    {initialize transformed data}
    for i := 0 to NData-1 do
      lny^[i] := ln(TDoubleArray(KnownY)[i]);

    {fit transformed data}
    LinEst16(lny^, KnownX, NData, LF, ErrorStats);

    {return values for B0 and B1 in exponential model y=B0*B1^x}
    LF.B0 := exp(LF.B0);
    LF.B1 := exp(LF.B1);
    {leave other values in LF in log form}
  finally
    freemem(lny, Size);
  end;
end;

function Forecast(X : Double; const KnownY: array of Double;
  const KnownX: array of Double) : Double;
begin
  if (High(KnownY) <> High(KnownX)) then
    RaiseStatError(stscStatBadCount);
  Result := Forecast16(X, KnownY, KnownX, High(KnownY)+1);
end;

function Forecast16(X : Double; const KnownY; const KnownX;
  NData : Integer) : Double;
var
  LF : TStLinEst;
begin
  LinEst16(KnownY, KnownX, NData, LF, false);
  Result := LF.B0+LF.B1*X;
end;

function ForecastExponential(X : Double; const KnownY: array of Double;
  const KnownX: array of Double) : Double;
begin
  if (High(KnownY) <> High(KnownX)) then
    RaiseStatError(stscStatBadCount);
  Result := ForecastExponential16(X, KnownY, KnownX, High(KnownY)+1);
end;

function ForecastExponential16(X : Double; const KnownY; const KnownX;
  NData : Integer) : Double;
var
  LF : TStLinEst;
begin
  LogEst16(KnownY, KnownX, NData, LF, false);
  Result := LF.B0*Power(LF.B1, X);
end;

function Intercept(const KnownY: array of Double;
  const KnownX: array of Double) : Double;
begin
  if (High(KnownY) <> High(KnownX)) then
    RaiseStatError(stscStatBadCount);
  Result := Intercept16(KnownY, KnownX, High(KnownY)+1);
end;

function Intercept16(const KnownY; const KnownX; NData : Integer) : Double;
var
  LF : TStLinEst;
begin
  LinEst16(KnownY, KnownX, NData, LF, false);
  Result := LF.B0;
end;

function RSquared(const KnownY: array of Double;
  const KnownX: array of Double) : Double;
begin
  if (High(KnownY) <> High(KnownX)) then
    RaiseStatError(stscStatBadCount);
  Result := RSquared16(KnownY, KnownX, High(KnownY)+1);
end;

function RSquared16(const KnownY; const KnownX; NData : Integer) : Double;
var
  LF : TStLinEst;
begin
  LinEst16(KnownY, KnownX, NData, LF, true);
  Result := LF.R2;
end;

function Slope(const KnownY: array of Double;
  const KnownX: array of Double) : Double;
begin
  if (High(KnownY) <> High(KnownX)) then
    RaiseStatError(stscStatBadCount);
  Result := Slope16(KnownY, KnownX, High(KnownY)+1);
end;

function Slope16(const KnownY; const KnownX; NData : Integer) : Double;
var
  LF : TStLinEst;
begin
  LinEst16(KnownY, KnownX, NData, LF, false);
  Result := LF.B1;
end;

function StandardErrorY(const KnownY: array of Double;
  const KnownX: array of Double) : Double;
begin
  if (High(KnownY) <> High(KnownX)) then
    RaiseStatError(stscStatBadCount);
  Result := StandardErrorY16(KnownY, KnownX, High(KnownY)+1);
end;

function StandardErrorY16(const KnownY; const KnownX;
  NData : Integer) : Double;
var
  LF : TStLinEst;
begin
  LinEst16(KnownY, KnownX, NData, LF, true);
  Result := LF.sigma;
end;

{------------------------------------------------------------------------}

function BetaCf(a, b, x : Single) : Single;
  {-Evaluates continued fraction for incomplete beta function}
const
  MAXIT = 100;
  EPS = 3.0E-7;
  FPMIN = 1.0E-30;
var
  m, m2 : Integer;
  aa, c, d, del, h, qab, qam, qap : Double;
begin
  qab := a+b;
  qap := a+1.0;
  qam := a-1.0;
  c := 1.0;
  d := 1.0-qab*x/qap;
  if (abs(d) < FPMIN) then
    d := FPMIN;
  d := 1.0/d;
  h := d;
  for m := 1 to MAXIT do begin
    m2 := 2*m;
    aa := m*(b-m)*x/((qam+m2)*(a+m2));
    d := 1.0+aa*d;
    if (abs(d) < FPMIN) then
      d := FPMIN;
    c := 1.0+aa/c;
    if (abs(c) < FPMIN) then
      c := FPMIN;
    d := 1.0/d;
    h := h*d*c;
    aa := -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));

    d := 1.0+aa*d;
    if (abs(d) < FPMIN) then
      d := FPMIN;
    c := 1.0+aa/c;
    if (abs(c) < FPMIN) then
      c := FPMIN;
    d := 1.0/d;
    del := d*c;
    h := h*del;

    {check for convergence}
    if (abs(del-1.0) < EPS) then
      break;
    if m = MAXIT then
      RaiseStatError(stscStatNoConverge);
  end;
  Result := h;
end;

function BetaDist(X, Alpha, Beta, A, B : Single) : Single;
var
  bt : Double;
begin
  if (X < A) or (X > B) or (A = B) or (Alpha <= 0.0) or (Beta <= 0.0) then
    RaiseStatError(stscStatBadParam);

  {normalize X}
  X := (X-A)/(B-A);

  {compute factors in front of continued fraction expansion}
  if (X = 0.0) or (X = 1.0) then
    bt := 0.0
  else
    bt := exp(GammaLn(Alpha+Beta)-GammaLn(Alpha)-GammaLn(Beta)+
      Alpha*ln(X)+Beta*ln(1.0-X));

  {use symmetry relationship to make continued fraction converge quickly}
  if (X < (Alpha+1.0)/(Alpha+Beta+2.0)) then
    Result := bt*BetaCf(Alpha, Beta, X)/Alpha
  else
    Result := 1.0-bt*BetaCf(Beta, Alpha, 1.0-X)/Beta;
end;

function Sign(a, b : Double) : Double;
  {-Returns abs(a) if b >= 0.0, else -abs(a)}
begin
  if (b >= 0.0) then
    Result := abs(a)
  else
    Result := -abs(a);
end;

function BetaInv(Probability, Alpha, Beta, A, B : Single) : Single;
const
  MAXIT = 100;
  UNUSED = -1.11e30;
  XACC = 3e-7;
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
  (A = B) or (Alpha <= 0.0) or (Beta <= 0.0) then
    RaiseStatError(stscStatBadParam);

  if (Probability = 0.0) then
    Result := A
  else if (Probability = 1.0) then
    Result := B
  else begin
    {Ridder's method of finding the root of BetaDist = Probability}
    fl := -Probability; {BetaDist(A, Alpha, Beta, A, B)-Probability}
    fh := 1.0-Probability; {BetaDist(B, Alpha, Beta, A, B)-Probability}
    xl := A;
    xh := B;
    ans := UNUSED;

    for j := 1 to MAXIT do begin
      xm := 0.5*(xl+xh);
      fm := BetaDist(xm, Alpha, Beta, A, B)-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 := BetaDist(ans, Alpha, Beta, A, B)-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;
    BetaInv := A; {avoid compiler warning}
    RaiseStatError(stscStatNoConverge);
  end;
end;

function BinomDistPmf(N, K : Integer; p : Extended) : Double;
  {-Returns the Probability mass function of the binomial distribution}
begin