unit2.pas

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

unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls,unit1, Forms, Dialogs;
Function GASDEV:real;
Procedure LFIT(X, Y, SIG:array of real; NDATA:integer;var A:array of real;
      MA:integer; LISTA:array of integer; MFIT:integer;var COVAR:matrx2;
      NCVM:integer; var CHISQ: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 GAUSSJ(VAR A:matrx2; N:integer; VAR B:array of real);
var
    IPIV,INDXR,INDXC:array[1..50] of integer;
    J,I,K,L,LL:integer;
    BIG,PIVINV,DUM:real; IROW,ICOL:integer;
begin
    For J:=1 To N do
      IPIV[J]:=0;
    For I:=1 To N do
    begin
        BIG:=0;
        For J:=1 To N do
        begin
            If IPIV[J] <> 1 Then
            begin
                For K:=1 To N do
                begin
                    If IPIV[K] = 0 Then
                    begin
                        If Abs(A[J, K]) >= BIG Then
                        begin
                            BIG:=Abs(A[J, K]);
                            IROW:=J;
                            ICOL:=K;
                        end;
                    end
                    Else if IPIV[K] > 1 Then
                        ShowMessage('Singular matrix.');
                end;
            end;
        end;
        IPIV[ICOL]:=IPIV[ICOL] + 1;
        If IROW <> ICOL Then
        begin
            For L:=1 To N do
            begin
                DUM:=A[IROW, L];
                A[IROW, L]:=A[ICOL, L];
                A[ICOL, L]:=DUM;
            end;
            DUM:=B[IROW];
            B[IROW]:=B[ICOL];
            B[ICOL]:=DUM;
        end;
        INDXR[I]:=IROW;
        INDXC[I]:=ICOL;
        If A[ICOL, ICOL] = 0 Then ShowMessage('Singular matrix.');
        PIVINV:=1 / A[ICOL, ICOL];
        A[ICOL, ICOL]:=1;
        For L:=1 To N do
          A[ICOL, L]:=A[ICOL, L] * PIVINV;
        B[ICOL]:=B[ICOL] * PIVINV;
        For LL:=1 To N do
        begin
            If LL <> ICOL Then
            begin
                DUM:=A[LL, ICOL];
                A[LL, ICOL]:=0;
                For L:=1 To N do
                    A[LL, L]:=A[LL, L] - A[ICOL, L] * DUM;
                B[LL]:=B[LL] - B[ICOL] * DUM;
            end;
        end;
    end;
    For L:=N DownTo 1 do
    begin
        If INDXR[L] <> INDXC[L] Then
        begin
            For K:=1 To N do
            begin
                DUM:=A[K, INDXR[L]];
                A[K, INDXR[L]]:=A[K, INDXC[L]];
                A[K, INDXC[L]]:=DUM;
            end;
        end;
    end;
end;

Procedure COVSRT(var COVAR:matrx2; NCVM, MA:integer;
                        LISTA:array of integer; MFIT:integer);
var
    I,J:integer;  SWAP:real;
begin
    For J:=1 To MA - 1 do
        For I:=J + 1 To MA do
            COVAR[I, J]:=0;
    For I:=1 To MFIT - 1 do
    begin
        For J:=I + 1 To MFIT do
        begin
            If LISTA[J] > LISTA[I] Then
                COVAR[LISTA[J], LISTA[I]]:=COVAR[I, J]
            Else
                COVAR[LISTA[I], LISTA[J]]:=COVAR[I, J];
        end;
    end;
    SWAP:=COVAR[1, 1];
    For J:=1 To MA do
    begin
        COVAR[1, J]:=COVAR[J, J];
        COVAR[J, J]:=0;
    end; 
    COVAR[LISTA[1], LISTA[1]]:=SWAP;
    For J:=2 To MFIT do
        COVAR[LISTA[J], LISTA[J]]:=COVAR[1, J];
    For J:=2 To MA do
        For I:=1 To J - 1 do
            COVAR[I, J]:=COVAR[J, I];
End;

Procedure LFIT(X, Y, SIG:array of real; NDATA:integer;var A:array of real;
      MA:integer; LISTA:array of integer; MFIT:integer;var COVAR:matrx2;
      NCVM:integer; var CHISQ:real);
var
    BETA,AFUNC:array[0..50] of real;
    I,J,K,KK,IHIT:integer;  YM,SIG2I,WT,SUM:real;
begin
    KK:=MFIT + 1;
    For J:=1 To MA do
    begin
        IHIT:=0;
        For K:=1 To MFIT do
            If LISTA[K] = J Then IHIT:=IHIT + 1;
        If IHIT = 0 Then
        begin
            LISTA[KK]:=J;
            KK:=KK + 1;
        end
        Else If IHIT > 1 Then
            ShowMessage(' Improper set in LISTA');
    end; 
    If KK <> (MA + 1) Then ShowMessage(' Improper set in LISTA');
    For J:=1 To MFIT do
    begin
        For K:=1 To MFIT do
            COVAR[J, K]:=0; 
        BETA[J]:=0; 
    end; 
    For I:=1 To NDATA do
    begin
        FUNCS(X[I], AFUNC, MA);
        YM:=Y[I];
        If MFIT < MA Then
        begin
            For J:=MFIT + 1 To MA do
                YM:=YM - A[LISTA[J]] * AFUNC[LISTA[J]];
        end;
        SIG2I:=1 / Sqr(SIG[I]);
        For J:=1 To MFIT do
        begin
            WT:=AFUNC[LISTA[J]] * SIG2I;
            For K:=1 To J do
                COVAR[J, K]:=COVAR[J, K] + WT * AFUNC[LISTA[K]];
            BETA[J]:=BETA[J] + YM * WT;
        end; 
    end; 
    If MFIT > 1 Then
    begin
        For J:=2 To MFIT do
            For K:=1 To J - 1 do
                COVAR[K, J]:=COVAR[J, K];
    end;
    GAUSSJ(COVAR, MFIT, BETA);
    For J:=1 To MFIT do
        A[LISTA[J]]:=BETA[J];
    CHISQ:=0; 
    For I:=1 To NDATA do
    begin
        FUNCS(X[I], AFUNC, MA);
        Sum:=0; 
        For J:=1 To MA do
            Sum:=Sum + A[J] * AFUNC[J];
        CHISQ:=CHISQ + Sqr((Y[I] - Sum) / SIG[I]);
    end; 
    COVSRT(COVAR, NCVM, MA, LISTA, MFIT);
end;
end.