unit2.pas

用于开发税务票据管理的软件

unit Unit2;

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

Function GASDEV:real;
Procedure FIT(X, Y:array of real; NDATA:integer; SIG:array of real;
               var MWT:integer; var A, B, SIGA, SIGB, CHI2, Q:real);
procedure MEDFIT(X, Y:array of real; NDATA:integer;var A, B, ABDEV: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;

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 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;

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;

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 FIT(X, Y:array of real; NDATA:integer; SIG:array of real;
     var MWT:integer; var A, B, SIGA, SIGB, CHI2, Q:real);
var
    I,J:integer;   SX,SY,ST2,SS,WT,SXOSS,T,SIGDAT:real;
begin
    SX:=0;
    SY:=0;
    ST2:=0;
    B:=0;
    If MWT <> 0 Then
    begin
        SS:=0;
        For I:=1 To NDATA do
        begin
            WT:=1  / (Sqr(SIG[I]));
            SS:=SS + WT;
            SX:=SX + X[I] * WT;
            SY:=SY + Y[I] * WT;
        end;
    end
    Else
    begin
        For I:=1 To NDATA do
        begin
            SX:=SX + X[I];
            SY:=SY + Y[I];
        end;
        SS:=NDATA;
    end;
    SXOSS:=SX / SS;
    If MWT <> 0 Then
    begin
        For I:=1 To NDATA do
        begin
            T:=(X[I] - SXOSS) / SIG[I];
            ST2:=ST2 + T * T;
            B:=B + T * Y[I] / SIG[I];
        end;
    end
    Else
    begin
        For I:=1 To NDATA do
        begin
            T:=X[I] - SXOSS;
            ST2:=ST2 + T * T;
            B:=B + T * Y[I];
        end;
    end;
    B:=B / ST2;
    A:=(SY - SX * B) / SS;
    SIGA:=Sqrt((1  + SX * SX / (SS * ST2)) / SS);
    SIGB:=Sqrt(1 / ST2);
    CHI2:=0;
    If MWT = 0 Then
    begin
        For I:=1 To NDATA do
            CHI2:=CHI2 + Sqr(Y[I] - A - B * X[I]);
        Q:=1;
        SIGDAT:=Sqrt(CHI2 / (NDATA - 2));
        SIGA:=SIGA * SIGDAT;
        SIGB:=SIGB * SIGDAT
    end
    Else
    begin
        For I:=1 To NDATA do
            CHI2:=CHI2 + Sqr((Y[I] - A - B * X[I]) / SIG[I]);
        Q:=GAMMQ(0.5 * (NDATA - 2), 0.5 * CHI2);
    end;
end;

Procedure SORT(N:integer; var RA:array of real);
Label 99;
var
    I,J,L,IR:integer;
    RR,RRA:real;
begin
    L:= N div 2 + 1;
    IR:=N;
    While true do
    begin
        If L > 1 Then
        begin
            L:=L - 1;
            RRA:=RA[L];
        end
        Else
        begin
            RRA:=RA[IR];
            RA[IR]:=RA[1];
            IR:=IR - 1;
            If IR = 1 Then
            begin
              RA[1]:=RRA;
              goto 99;
            end;
        end;
        I:=L;
        J:=L + L;
        While J <= IR do
        begin
            If J < IR Then
                If RA[J] < RA[J + 1] Then J:=J + 1;
            If RRA < RA[J] Then
            begin
                RA[I]:=RA[J];
                I:=J;
                J:=J + J;
            end
            Else
                J:=IR + 1;
        end;
        RA[I]:=RRA;
    end;
99: end;

Function ROFUNC(B:real):real;
var
    N1,NM2,NMH,NML,J:integer; SUM,D,ZZ:real;
begin
    N1:=NDATAT + 1;
    NML:=N1 DIV 2;
    NMH:=N1 - NML;
    For J:=1 To NDATAT do
        ARR[J]:=YT[J] - B * XT[J];
    SORT(NDATAT, ARR);
    AA:=0.5 * (ARR[NML] + ARR[NMH]);
    Sum:=0; 
    ABDEVT:=0;
    For J:=1 To NDATAT do
    begin
        D:=YT[J] - (B * XT[J] + AA);
        ABDEVT:=ABDEVT + Abs(D);
        If D >= 0 then
            ZZ:=1
        else
            ZZ:=-1; 
        Sum:=Sum + XT[J] * ZZ;
    end; 
    ROFUNC:=Sum;
end;

procedure MEDFIT(X, Y:array of real; NDATA:integer; var A, B, ABDEV:real);
Label 1,2,3;
var
    J:integer;
    SX,SY,SXY,SXX,DEL,AA,BB,CHISQ,SIGB,B1,F1,B2,F2,ZZ,F:real;
begin
    SX:=0; 
    SY:=0;
    SXY:=0; 
    SXX:=0; 
    For J:=1 To NDATA do
    begin
        XT[J]:=X[J];
        YT[J]:=Y[J];
        SX:=SX + X[J];
        SY:=SY + Y[J];
        SXY:=SXY + X[J] * Y[J];
        SXX:=SXX + Sqr(X[J]);
    end; 
    NDATAT:=NDATA;
    DEL:=NDATA * SXX - Sqr(SX);
    AA:=(SXX * SY - SX * SXY) / DEL;
    BB:=(NDATA * SXY - SX * SY) / DEL;
    CHISQ:=0;
    For J:=1 To NDATA do
        CHISQ:=CHISQ + Sqr(Y[J] - (AA + BB * X[J]));
    SIGB:=Sqrt(CHISQ / DEL);
    B1:=BB;
    F1:=ROFUNC(B1);
    If f1 >= 0 then
        ZZ:=1
    Else
        ZZ:=-1;
    B2:=BB + Abs(3 * SIGB) * ZZ;
    F2:=ROFUNC(B2);
1:  If F1 * F2 > 0  Then
    begin
        BB:=2 * B2 - B1;
        B1:=B2;
        F1:=F2;
        B2:=BB;
        F2:=ROFUNC(B2);
        GoTo 1;
    end;
    SIGB:=0.01 * SIGB;
2:  If Abs(B2 - B1) > SIGB Then
    begin
        BB:=0.5 * (B1 + B2);
        If (BB = B1) Or (BB = B2) Then GoTo 3;
        F:=ROFUNC(BB);
        If F * F1 >= 0  Then
        begin
            F1:=F;
            B1:=BB;
        end
        Else
        begin
            F2:=F;
            B2:=BB;
        end;
        GoTo 2;
    end;
3:  A:=AA;
    B:=BB;
    ABDEV:=ABDEVT / NDATA;
end;
end.