unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics,unit1, Controls, Forms, Dialogs;

Function GASDEV:real;
procedure AVEVAR(DATA:array of real; N:integer; var AVE, VAR1:real);
procedure TPTEST(DATA1, DATA2:array of real; N:integer;var T, PROB:real);

implementation
Function GASDEV:real;
var
   V1,V2,FAC,R:real;
begin
    If ISET^= 0 Then
    begin
      repeat
        V1:=2 * Random - 1;
        V2:=2 * Random - 1;
        R:=Sqr(V1) + Sqr(V2);
      until (R < 1);
      FAC:=Sqrt(-2 * Ln(R) / R);
      GSET^:=V1 * FAC;
      GASDEV:=V2 * FAC;
      ISET^:=1;
    end
    Else
    begin
      GASDEV:=GSET^;
      ISET^:=0;
    end;
end;

procedure AVEVAR(DATA:array of real; N:integer; var AVE, VAR1:real);
var
    J:integer;  S:real;
begin
    AVE:=0; 
    VAR1:=0;
    For J:=1 To N do
        AVE:=AVE + DATA[J];
    AVE:=AVE / N;
    For J:=1 To N do
    begin
        S:=DATA[J] - AVE;
        VAR1:=VAR1 + S * S;
    end;
    VAR1:=VAR1 / (N - 1);
end;

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;

Function BETACF( A, B, X:real):real;
label 1;
const
  ITMAX=100;    EPS=0.0000003;
var
  TEM,QAP,QAM,QAB,EM,D,BZ,BP,BPP,BM,AZ,AAP,AM,AOLD,AP:real;
  M:INTEGER;
begin
  AM:=1;
  BM:=1;
  AZ:=1;
  QAB:=A + B;
  QAP:=A + 1;
  QAM:=A - 1;
  BZ:=1 - QAB * X / QAP;
  For M:=1 To ITMAX do
  begin
    EM:=M;
    TEM:=EM + EM;
    D:=EM * (B - M) * X / ((QAM + TEM) * (A + TEM));
    AP:=AZ + D * AM;
    BP:=BZ + D * BM;
    D:=-(A + EM) * (QAB + EM) * X / ((A + TEM) * (QAP + TEM));
    AAP:=AP + D * AZ;
    BPP:=BP + D * BZ;
    AOLD:=AZ;
    AM:=AP / BPP;
    BM:=BP / BPP;
    AZ:=AAP / BPP;
    BZ:=1;
    If Abs(AZ - AOLD) < EPS * Abs(AZ) Then GoTo 1;
  end;
  ShowMessage('A or B too big, or ITMAX too small');
1: BETACF:=AZ;
end;

Function BETAI(A,B,X:real):real;
var
  AAA,BT:real;
begin
  If (X < 0) Or (X > 1) Then ShowMessage('bad argument X in BETAI');
  If (X=0) Or (X=1) Then
    BT:=0
  Else
    begin
      AAA:=GAMMLN(A + B) - GAMMLN(A) - GAMMLN(B);
      BT:=Exp(AAA + A * Ln(X) + B * Ln(1 - X));
    end;
  If X < (A + 1) / (A + B + 2) Then
    BETAI:=BT * BETACF(A, B, X) / A
  Else
    BETAI:=1 - BT * BETACF(B, A, 1 - X) / B;
end;

procedure TPTEST(DATA1, DATA2:array of real; N:integer;var T,PROB:real);
var
    COV,DF,AVE1,VAR1,AVE2,VAR2,SD:real;  J:integer;
begin
    AVEVAR(DATA1, N, AVE1, VAR1);
    AVEVAR(DATA2, N, AVE2, VAR2);
    COV:=0;
    For J:=1 To N do
        COV:=COV + (DATA1[J] - AVE1) * (DATA2[J] - AVE2);
    DF:=N - 1;
    COV:=COV / DF;
    SD:=Sqrt((VAR1 + VAR2 - 2  * COV) / N);
    T:=(AVE1 - AVE2) / SD;
    PROB:=BETAI(0.5 * DF, 0.5, DF / (DF + T*T));
end;

end.