lfit.txt

《Delphi常用数值算法集》的配书源码

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;