mrqmin.dem

Delphi Pascal 数据挖掘领域算法包 数值算法大全

PROGRAM D14R8(input,output);
(* driver for routine MRQMIN *)
LABEL 1;
CONST
   npt=100;
   mma=6;
   spread=0.001;
TYPE
   glmma = ARRAY [1..mma] OF real;
   glnparam = glmma;
   gllista = ARRAY [1..mma] OF integer;
   glcovar = ARRAY [1..mma,1..mma] OF real;
   glnalbynal = glcovar;
   glncabynca = glcovar;
   glnpbynp = glcovar;
   glnpbymp = glcovar;
   glndata = ARRAY [1..npt] OF real;
VAR
   glochisq : real;
   glatry,glbeta : glmma;
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   gliset : integer;
   glgset : real;
   alamda,chisq,ochisq : real;
   i,idum,itst,j,jj,k,mfit : integer;
   x,y,sig : glndata;
   lista : gllista;
   a,gues : glmma;
   covar,alpha : glcovar;

(*$I MODFILE.PAS *)
(*$I COVSRT.PAS *)

(*$I RAN3.PAS *)

(*$I GASDEV.PAS *)

PROCEDURE funcs(x: real; a: glnparam; VAR y: real;
       VAR dyda: glnparam; na: integer);
(* Programs using routine FGAUSS must define the type
TYPE
   glnparam = ARRAY [1..na] OF real;
where na is the number of parameters. *)
VAR
   i,ii : integer;
   fac,ex,arg : real;
BEGIN
   y := 0.0;
   FOR ii := 1 to (na DIV 3) DO BEGIN
      i := 3*ii-2;
      arg := (x-a[i+1])/a[i+2];
      ex := exp(-sqr(arg));
      fac := a[i]*ex*2.0*arg;
      y := y+a[i]*ex;
      dyda[i] := ex;
      dyda[i+1] := fac/a[i+2];
      dyda[i+2] := fac*arg/a[i+2]
   END
END;

(*$I GAUSSJ.PAS *)

(*$I MRQCOF.PAS *)

(*$I MRQMIN.PAS *)

BEGIN
   gliset := 0;
   a[1] := 5.0; a[2] := 2.0; a[3] := 3.0;
   a[4] := 2.0; a[5] := 5.0; a[6] := 3.0;
   gues[1] := 4.5; gues[2] := 2.2; gues[3] := 2.8;
   gues[4] := 2.5; gues[5] := 4.9; gues[6] := 2.8;
   idum := -911;
   FOR i := 1 to 100 DO BEGIN
      x[i] := 0.1*i;
      y[i] := 0.0;
      FOR jj := 1 to 2 DO BEGIN
         j := 3*jj-2;
         y[i] := y[i]+a[j]*exp(-sqr((x[i]-a[j+1])/a[j+2]))
      END;
      y[i] := y[i]*(1.0+spread*gasdev(idum));
      sig[i] := spread*y[i]
   END;
   mfit := 6;
   FOR i := 1 to mfit DO lista[i] := i;
   alamda := -1;
   FOR i := 1 to mma DO a[i] := gues[i];
   mrqmin(x,y,sig,npt,a,mma,lista,mfit,covar,alpha,
      mma,chisq,alamda);
   k := 1;
   itst := 0;
1:   writeln;
   writeln('Iteration #',k:2,'chi-squared:':17,chisq:10:4,
      'alamda:':10,alamda:9);
   writeln('a[1]':7,'a[2]':8,'a[3]':8,'a[4]':8,'a[5]':8,'a[6]':8);
   FOR i := 1 to 6 DO write(a[i]:8:4);
   writeln;
   k := k+1;
   ochisq := chisq;
   mrqmin(x,y,sig,npt,a,mma,lista,mfit,covar,alpha,
      mma,chisq,alamda);
   IF (chisq > ochisq) THEN BEGIN
      itst := 0
   END ELSE IF (abs(ochisq-chisq) < 0.1) THEN BEGIN
      itst := itst+1
   END;
   IF (itst < 2) THEN GOTO 1;
   alamda := 0.0;
   mrqmin(x,y,sig,npt,a,mma,lista,mfit,covar,alpha,
      mma,chisq,alamda);
   writeln('Uncertainties:');
   FOR i := 1 to 6 DO write(sqrt(covar[i,i]):8:4);
   writeln
END.