Non-linear fitting is defined here to be the minimization of
In equation (1) fi is the i'th data point associated with the position xi which is frequently one dimensional, but need not be. The relative "error" associated with this point is . The derivations are in Extremal.htm. The key routine is annotated in the document ROBMIN.doc. Plot routines are described in gplot\WELCOME.htm. For windows, these are windows\winplot.zip works with 256 colors and win32\gplot.zip for those with more colors. In general NLFIT reads a file of the form x,y and creates a file fit.dat of the form x,y. It also produces err.out, abserr.out and res.out files which are respectively (data – fit), abs(data-fit), and (data-fit)/e.
The user needs to write and link to the main code a routine poly.for that returns fA and (p(j), j=1,Nhat), where p is the partial of fA with respect to each of the Nhat constants. A sample listing of possible approximating function is in fA(c,x) . Poly subroutines, annotated below, using the analytic derivatives of the generalized Fermi function F(x) and of F’(x) are in the routines FPOLY.FOR, DFPOLY.FOR while polynomials using numerical derivatives of these same functions are in FNDPOLY.FOR, and DFNDPOLY.FOR. The use of numerical derivatives slows the code by approximately a factor of 2 times the number of derivatives involved.
– 1/4/2008 Numerical multiple precision derivatives will be much slower, but will have minimal subtraction problems and hence much more accuracy.
The main Fortran routine is formp\NLFIT.FOR which is annotated below in main. This routine ends with a number of includes. It can be run from the watfor editor by the simple command run. An executable nlfit.exe can be made by the command run/exe. This will then be able to use a command line name for the dir file as described in Wsystem.doc htm.
The routine is not even the fastest Watfor routine owing to its array bounds checking. The routine W87 will need to be configured as described in the reference above to look for the library. Then W87 /exe nlfit will produce a reasonably fast executable. The file README.htm mentions the bat file mkdffit.bat.txt for producing MS fort 5.2 files àThe text version now has fl32 to produce Fortran Power Station files. The last file in the compilation list is formp\PSSYS.FOR. This file and similar ones for other compilers are discussed in ..\..\MultiplePrecision\Bob\introduction.doc htm. Bat files are somewhat machine dependent and should be considered as illustrative only.
Once the appropriate poly routine has been entered, the fastest route into nlfit is described in the file read.me which can also be reached by entering the commands nlfit ? or nlfit help.
The following is a few portions of the code formp\NLFIT.FOR with annotations.
C *** THIS ROUTINE DOES A WEIGHTED NON-LINEAR LEAST SQUARES FIT TO
C *** USER SUPPLIED DATA IN L DIMENSIONS. THE ROUTINE POLY MUST
C *** RETURN THE PARTIALS OF THE FITTING FUNCTIONS WITH RESPECT TO
C *** THE CONSTANTS BEING VARIED.
C *** DATA IS SUPPLIED IN FREE FORMAT AS X1,FXY,X2,X3,...
This form makes multi-dimensional data plot as a series of one dimensional plots of fxy, versus x1. Option 12 in ..\..\gplot\WELCOME.DOC .htm stops the connection of points sufficiently far apart to achieve this.
PARAMETER(MAXCON=20,MAXDIM=2,IT=63)
DIMENSION P(MAXCON),CONS(MAXCON),IVORF(MAXCON)
DIMENSION MPP(IT+1,MAXCON),MPTEMP(IT+1)
DIMENSION X(MAXDIM),XS(MAXDIM),XB(MAXDIM),IFL(MAXCON)
CHARACTER*64 NA,NADIR,NADAT,STATUS,NAOUT,CTEMP
CHARACTER*1 CSTOP
COMMON/PASSRM/VCONS(MAXCON),PC(MAXCON),MPPC(IT+1,MAXCON),
2 PPCC(MAXCON*(MAXCON+1)/2),DUM(2*MAXCON+1),
3 MPPPCC(IT+1,MAXCON*(MAXCON+1)/2),MP
C *** VCONS --> CONS IN ROBMIN, PC, PPCC, ARE USED IN FOFLAM AND THUS NEED TO
C *** BE PASSED IN COMMON - DUM --> SM IS DEFINED IN ROBMIN AND USED BY FOFLAM
C *** DUM - DOES NOT APPEAR IN NLFIT
DATA XS,XB/MAXDIM*1.D32,MAXDIM*-1.D32/
DATA CONS/MAXCON*0.D0/
c testing mpsminv
c OPEN(3,FILE='TESTA.OUT')
OPEN(4,FILE='FAIL.OUT')
C OPEN(1,FILE='SMINV.OUT')
C SCRATCH FILES
CALL MAOPUF(9)
CALL MAOPUF(17)
CALL MAOPUF(19)
These files are supposed to be deleted by the system on completion of the code. Frequently they are left over as zz.. files or for09 files. The user may need to delete these manually now and then.
CTEST OPEN(21,FILE='TEST.OUT')
CHLAST=1D67
C *** READING THE NAME OF THE INPUT FILE FROM THE COMMAND LINE
NUMARGS=NARGS()
IF(NUMARGS.GT.1)THEN
NTEMP=1
CALL BGETARG(NTEMP,NADIR)
PRINT*,' NADIR=',NADIR
ELSE
PRINT*,' ENTER THE NAME OF THE DIRECTION FILE or ? for help'
READ(*,'(A)')NADIR
ENDIF
C *** GIVING SCREEN HELP
IF(NADIR(1:4).EQ.'HELP'.OR.NADIR(1:4).EQ.'help'.OR.
2 NADIR(1:1).EQ.'?')THEN
NADIR='ReadMe.htm'
CALL FSYSTEM(NADIR)
STOP 'requested help'
Fsystem acts as though it were a command line. It is different for every compiler. It is contained in the files formp\intelsys.for, formp\openwsys.for, formp\PSSYS.FOR, formp\unixsys.for and formp\WSYSTEM.FOR. These also contain the various opens and other commands needed to operate with different compilers.
ENDIF
C *** READING THE DIRECTION FILE
The code at this point reads the direction file a sample of which is in formp\aigau.htm
Data is read as formatted from unit 8 and written unformatted to unit 9 for later use. Error coefficients are discussed in Welcome.doc#DirectionFile and in The error form.doc .htm.
8 CONTINUE
IF(IEP.EQ.1)THEN
READ(8,*,END=10)X(1),FXY,(X(I),I=2,ND)
EP=SQRT(AERR+BERR*ABS(FXY))+CERR*ABS(FXY)
ELSE
READ(8,*,END=10)X(1),FXY,(X(I),I=2,ND),EP
ENDIF
WRITE(9)(X(I),I=1,ND),FXY,EP
GOTO 8
10 CLOSE(8)
The minimization loop follows. The routine formp\poly.for returns FA and the partials of FA with respect to each of the constants. The parameter MP, this is the option on line 6 of the direction file 400,1 in formp\ReadMe.htm. If MP is ³ 3, POLY should return a multiple precision FA and multiple precision partials. The coefficients always remain double precision.
After the minimization has been completed the matrix ppcc is inverted without exp(l) to form an “honest” error matrix.
The quote refers to the fact that this inversion continues even if the matrix is technically singular with a zero or less diagonal element. When this singularity occurs parts of the array are linearly dependent which means that the coefficients associated with these parts of the array can cancel each other leading to infinite errors. Changes in these coefficients are technically made zero and the code proceeds to estimate the rest of the errors. I have found that these estimates are reasonable reliable in giving what we like to think of as error estimates. The coefficients leading to the singularities are written out with errors of 10-37. The presence of these values implies that all of the error estimates are suspicious and that the fit contains linear dependencies at this level. LinDep.doc.
IF(MP.GE.1)THEN
CALL MPHSMINV(PPCC,MPPPCC,NVMIN,IFL,IFLT,1)
ELSE
CALL SMINV(PPCC,NVMIN,IFL,IFLT,1)
ENDIF
The call to SMINV makes PPCC the error matrix. This is used to create the next direction file with the name given on line 1 of the direction file. In principle, this can be the same name as the input direction file, but I recommend using a different name so that it is easy to take a single step backwards.
The section beginng with
NA='FIT.DAT'
Forms fit.dat, err.out and res.out from the data and the fit. It is designed to be such that plotting uses the command
gplot Laigau.out fit.dat
The data is in Laigau.out and the fit is in fit.dat. A successful fit shows one on top of the other. In fact many fits need more information than this.
gplot res.out
shows ((fi -fa)/erri)2Npow while
gplot err.out
shows fi -fa.