我看见一个例子,不知如何?与我想要是不是一样.要怎样改
unit Unit1;
interface
uses
Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
StdCtrls;
procedure FUNCS(X:real;var AFUNC:array of real
MA:integer);
type
TForm1 = class(TForm)
Button1: TButton;
Memo1: TMemo;
procedure Button1Click(Sender: TObject);
private
{ Private declarations }
public
{ Public declarations }
end;
type
matrx2=array of array of real;
var
Form1: TForm1;
ISET:^integer
gset:^real;
implementation
//PROGRAM D7R2
//Driver for routine LFIT
uses
unit2;
{$R *.DFM}
procedure FUNCS(X:real;var AFUNC:array of real
MA:integer);
var
I:integer;
begin
AFUNC[1]:= 1;
For I:= 2 To MA do
AFUNC:= X * AFUNC[I - 1];
end;
procedure TForm1.Button1Click(Sender: TObject);
const
s1='%14.6f'
s2='.#####0e+00'
s3='#'
s4='.#0e+00'
s5='%9.6f';
NPT = 100
SPREAD = 0.1
NTERM = 3;
var
COVAR:matrx2;
A,X,Y,SIG:array[0..100] of real
LISTA:ARRAY[0..100] of integer;
I,J,MFIT,II:integer
AAA,CHI2,Q,CHISQ:real;
F:TextFile;
begin
SetLength(COVAR,4,4);
New(ISET)
New(GSET);
Iset^:=0
GSET^:=0;
Randomize;
For I:=1 To NPT do
begin
X:=0.1 * I;
Y:=NTERM;
For J:=NTERM-1 DownTo 1 do
Y:= J + Y*X;
Y:= Y+SPREAD*GASDEV;
SIG:=SPREAD;
end;
MFIT:=NTERM;
For I:= 1 To MFIT do
LISTA:= I;
LFIT(X, Y, SIG, NPT, A, NTERM, LISTA, MFIT, COVAR, NTERM, CHISQ);
//输出计算结果到文件
AssignFile(F, 'd:/delphi_shu/p7/d7r2.dat');
Rewrite(F);
Writeln(F,'Parameter Uncertainty');
For I:= 1 To NTERM do
Writeln(F,'A(',FormatFloat(s3,I), ')= ', Format(s5,[A]),
Format(s1,[Sqrt(COVAR[I, I])]));
Writeln(F);
Writeln(F,'Chi-squared = ', FormatFloat(s2,CHISQ));
Writeln(F);
Writeln(F,'Full covariance matrix');
For I:= 1 To NTERM do
Writeln(F,FormatFloat(s4,COVAR[I, 1]),' ',FormatFloat(s4,COVAR[I, 2]),
' ',FormatFloat(s4,COVAR[I, 3]));
//Now test the LISTA feature
For I:= 1 To NTERM do
LISTA:= NTERM + 1 - I;
LFIT(X, Y, SIG, NPT, A, NTERM, LISTA, MFIT, COVAR, NTERM, CHISQ);
Writeln(F);
Writeln(F,'Parameter Uncertainty');
For I:= 1 To NTERM do
Writeln(F,'A(',FormatFloat(s3,I), ')= ', Format(s5,[A]),
Format(s1,[Sqrt(COVAR[I, I])]));
Writeln(F);
Writeln(F,'Chi-squared = ', FormatFloat(s2,CHISQ));
Writeln(F);
Writeln(F,'Full covariance matrix');
For I:= 1 To NTERM do
Writeln(F,FormatFloat(s4,COVAR[I, 1]),' ',FormatFloat(s4,COVAR[I, 2]),
' ',FormatFloat(s4,COVAR[I, 3]));
//Now check results of restricting fit s
II:= 1;
For I:= 1 To NTERM do
begin
AAA:= I - (I Div 2) * 2;
If AAA = 1 Then
begin
LISTA[II]:= I;
II:= II + 1;
end;
end;
MFIT:= II - 1;
LFIT(X, Y, SIG, NPT, A, NTERM, LISTA, MFIT, COVAR, NTERM, CHISQ);
Writeln(F);
Writeln(F,'Parameter Uncertainty');
For I:= 1 To NTERM do
Writeln(F,'A(',FormatFloat(s3,I), ')= ', Format(s5,[A]),
' ',Format(s1,[Sqrt(COVAR[I, I])]));
Writeln(F);
Writeln(F,'Chi-squared = ', FormatFloat(s2,CHISQ));
Writeln(F);
Writeln(F,'Full covariance matrix');
For I:= 1 To NTERM do
Writeln(F,FormatFloat(s4,COVAR[I, 1]),' ',FormatFloat(s4,COVAR[I, 2]),
' ',FormatFloat(s4,COVAR[I, 3]));
CloseFile(F);
//屏幕显示计算结果
memo1.Lines.LoadFromFile('d:/delphi_shu/p7/d7r2.dat');
end;
end.
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:=IROW;
INDXC:=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 Then
COVAR[LISTA[J], LISTA]:=COVAR[I, J]
Else
COVAR[LISTA, 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, AFUNC, MA);
YM:=Y;
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);
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, AFUNC, MA);
Sum:=0
For J:=1 To MA do
Sum:=Sum + A[J] * AFUNC[J];
CHISQ:=CHISQ + Sqr((Y - Sum) / SIG);
end
COVSRT(COVAR, NCVM, MA, LISTA, MFIT);
end;
end.
