unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs;
procedure CHSTWO(BINS1, BINS2:array of real;
                         NBINS, KNSTRN:integer;var DF, CHSQ, PROB:real);

implementation
Function GAMMLN(xx:real):real;
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(var GAMMCF:real; A:real; X:real; var GLN:real);
label 1;
const
  ITMAX=100;    EPS=0.0000003;
var
  N:integer; GOLD,G,FAC,B1,B0,ANF,ANA,AN,A1,A0:real;
begin
  GLN:=GAMMLN(A);
  GOLD:=0;
  A0:=1;
  A1:=X;
  B0:=0;
  B1:=1;
  FAC:=1;
  For N:=1 To ITMAX do
  begin
    AN:=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 Then
    begin
      FAC:=1 / A1;
      G:=B1 * FAC;
      If Abs((G - GOLD)/ G)< EPS Then GoTo 1;
      GOLD:=G;
    end
  end;
  ShowMessage('A too large, ITMAX too small');
1: GAMMCF:=Exp(-X + A * Ln(X)- GLN)* G;
end;

procedure GSER(var GAMSER:REAL;A:real;X:REAL;var GLN:REAL);
label 1;
const
  ITMAX=100;    EPS=0.3E-6;
var
  N:Integer;
  SUM,DEL,AP:Real;
begin
  GLN:=GAMMLN(A);
  If X <= 0 Then
    If X < 0 Then
      begin
        ShowMessage('警告：X<0，退出程序。');
        goto 1;
      end
  else
    begin
      GAMSER:=0;
      goto 1;
    end;
  AP:=A;
  Sum:=1 / A;
  DEL:=Sum;
  For N:=1 To ITMAX do
    begin
      AP:=AP + 1;
      DEL:=DEL * X / AP;
      Sum:=Sum + DEL;
      If Abs(DEL) < Abs(Sum) * EPS Then GoTo 1;
    end;
  ShowMessage('A too large, ITMAX too small');
1: GAMSER:=Sum * Exp(-X + A * Ln(X) - GLN);
end;

Function GAMMQ(A,X:REAL):REAL;
VAR
  GAMSER,GLN:REAL;
BEGIN
    If (X < 0) Or (A <= 0) Then  ShowMessage('PAUSE');
    If X < A + 1 Then
      begin
        GSER(GAMSER, A, X, GLN);
        GAMMQ:=1-GAMSER;
      end
    Else
      begin
        GCF(GAMSER, A, X, GLN);
        GAMMQ:=GAMSER;
      end;
END;

procedure CHSTWO(BINS1, BINS2:array of real;
                         NBINS, KNSTRN:integer;var DF, CHSQ, PROB:real);
var
    J:integer;
begin
    DF:=NBINS - 1 - KNSTRN;
    CHSQ:=0; 
    For J:=1 To NBINS do
    begin
        If (BINS1[J] = 0) And (BINS2[J] = 0)  Then
          DF:=DF - 1 
        Else
          CHSQ:=CHSQ + Sqr(BINS1[J] - BINS2[J]) / (BINS1[J] + BINS2[J]);
    end; 
    PROB:=GAMMQ(0.5 * DF, 0.5 * CHSQ);
end;
end.