procedure (xy,ycol,fcol,mit,maxtry,b,init,npar,a);

Last revision 5 October 1990 CGE

Marquardt-Levenberg Non-Linear Least Squares Regression Procedure

Four items are needed or created to operate this program:

1. A function which must have the name func and returns a single value

2. A patameter table

3. A correlations table

4. A data table

The procedure passes several arguments and is called in the following manner;

MRQFIT('DATA',YCOL,FCOL,MIT,MAXTRY{,'PARAMETERS',INITIALIZE,NUMPAR,'

CORRELATIONS'})

The first five (5) arguments ('DATA',YCOL,FCOL,MIT AND MAXTRY) are

required. The remaining arguments are optional if not supplied are provided

with defalt values.

'DATA' is the name of the data table containing the dependent and

independent variable as columns. This argument may be passes as an RS/1

variable containing the table name as text or the table name enclosed in quotes

(") or apostrophes (').

YCOL is the name (enclosed in quotes (of the column in 'data'

corresponding to the dependent variable) or colunm number.

FCOL is the name (enclosed in quotes (of the column in 'data') or

column number into which the computed (FUNC) values will be stored. The 'data'

table must have space for the calculated values.

MIT is the maximum number of iterations. 5 is suggested.

MAXTRY is the maximum trial solutions per iteration. A range of 10 to

50 is suggested. Smaller valures will reduce exicution time but may not

converge.

The last four arguments are optional.

'PARAMETERS' is the name of a table which contains (or will contain)

the parameters as discussed previously. If the name is not provided, the name

PARAMETERS will be used. If the table does not exist it will be created and

initialized.

The argument initialize is set to true when the table 'PARAMETERS' is

CREATED by MRQFIT. When the table already exists, it will be initialized if

initialized is true, or used as is if initialized is false.

NUPAR is the number of parameters. This value is required only when

the table 'PARAMETERS' is to be initialized. If NUMPAR is not provided, it

will be requested.

'CORRELATIONS' is the name of table which will be stored the results of

the least squares data reduction. This table is an information table which is

restructured and overwritten by MRQFIT. If the name is not provided the name

CORRELATIONS will be used.

The Function: The function is a set of statements which contain the

variables to fit the data. The function must be named FUNC.

FUNC has the form:

procedure(x,p,fval);

ind = 1;

fval=b6(p[1],p[2],p[3],p[4],p[5],p[6],p[7],x[ind]);

end;

p refers to the parameter table the numbers in the [] refer to

the row number. x is the table containing the independent

variables the number in the [] is the column containing the

independent variable. The order that the parameters or

independent variables occur in the function ie. b6 is the same

as in the procedure for that function

DATA EXAMPLE:

Sample No Void Vol C/Co C/Co comp

1 .1 .001 .001

2 .2 .02 .019

etc

PARAMETERS EXAMPLE:

PARAMETER VALUE MINIMUM MAXIMUM STD DEV STUDENT T LOWER UPPER

KD 0 0 0

VD 4.5 4.5 4.5

DW .7 .01 1 .03 20 .6 .8

THETA .4 .4 .4

RHO 1.6 1.6 1.6

L 10 10 10

Xo 4 4 4

The parameter table needs to be in the same order as the function

arguments.

stopcr=0.0005;

nobs=lastrow(xy);

/* ---- Set-up Parameter Table ---- */

begin;

if empty(b) then b='parameters';

if not tableexists(b) then begin;

if empty(npar) then

npar=getnumber('Number of parameters:',false);

call maketable(b,npar,7,100,'Function Parameters');

b[0,1]='Parameter';

b[0,2]='Minimum';

b[0,3]='Maximum';

b[0,4]='Std Dev.';

b[0,5]="Student's t";

b[0,6]='lower';

b[0,7]='upper';

init=true;

end;

else if empty(init) then init=yesanswer('Initialize parameters:',

true,'Yes or No');

if init then begin;

if empty(npar) then

npar=getnumber('Number of parameters:',false);

/* INPUT PARAMETERS FOR MODEL */

erase;

type 'Initial parameter estimates and bounds.';

type 'Respond to the request for minimum and maximum parameter';

type 'values with:';

type ' 1) the same value if the parameter is to be constant,';

type ' 2) the bounds of the parameter if it is to bounded,';

type ' 3) EMPTY if the parameter is unconstrained.';

do i=1 to npar;

type nocr "For parameter #",i;

b[i,0]=gettext('Name:');

/* ----- READ CONSTRAINTS ON PARAMETER VALUES:

Columns 2 and 3 of B are minimum and maximum constraints.

B[I,2] or B[I,3] empty WILL OMIT CONSTRAINT ON PARAMETER I.

b[i,3]=getnumber('maximum value:',true,empty,

'Return (enter) if parameter is unconstrained above');

if notempty(b[i,3]) then begin;

b[i,2]=getnumber('minimum value (=empty if no constraint):',

true,b[i,3],

'Type =EMPTY if parameter is unconstrained below.');

if notempty(b[i,2]) then begin;

b[i,1]=(b[i,2]+b[i,3])/2.0;

if(b[i,2] eq b[i,3]) then donext;

end;

else b[i,1] = b[i,3];

end;

else begin;

b[i,2]=getnumber('minimum value:',true,empty,

'Return (enter) if parameter is unconstrained');

if notempty(b[i,2]) then

b[i,1] = b[i,2];

else b[i,1] = 0.0;

end;

/* ----- READ INITIAL ESTIMATES ----- */

b[i,1]=getnumber('Initial value:',false,b[i,1],

'First guess of parameter value.');

end;

else if empty(npar) then npar=lastrow(b);

/* ----- Allocate Parameter ARRAYS ----- */

notes of table(b) = empty;

set col 4 of table(b) = empty;

set col 5 of table(b) = empty;

set col 6 of table(b) = empty;

set col 7 of table(b) = empty;

r=vec$alloc(nobs);

th=vec$alloc(npar);

tb=vec$alloc(npar);

NP=0;

DO I=1 to Npar;

if notempty(b[i,2]) and notempty(b[i,3]) and

(b[i,2] EQ b[i,3]) then donext;

NP=NP+1;

iii[np]=i;

TB[NP]=B[I,1];

TH[NP]=B[I,1];

end;

/* Allocate work arrays */

q=vec$alloc(np);

e=vec$alloc(np);

p=vec$alloc(np);

phi=vec$alloc(np);

d=mat$alloc(np,np);

delz=mat$alloc(nobs,np);

if empty(a) then a='correlations';

if tableexists(a) then delete table(a);

call maketable(a,np,np,100);

do i=1 to np;

j=iii[i];

a[i,0]=b[j,0];

a[0,i]=a[i,0];

end;

/* ---- Initialize Marquadt Sum of Squares ---- */

GA=0.02;

NIT=0;

NP2=2*NP;

SSQ=0.0;

DO I=1 to Nobs;

call func(vec$t2v(row i of table(xy)),

vec$t2v(col 1 of table(b)),fvalue);

xy[i,fcol]=fvalue;

R[I]=xY[I,ycol]-Fvalue;

SSQ=SSQ+R[I]*R[I];

end;

IF (MIT EQ 0) GOTO s140;

/* ----- BEGIN ITERATION ----- */

s34:

NIT=NIT+1;

NTRIAL=0;

GA=0.1*GA;

/* ---- Estimate Parameter Derivatives ---- */

DO J=1 to NP;

k=iii[j];

b[k,1]=1.01*TH[J];

Q[J]=0.0;

DO I=1 to Nobs;

call func(vec$t2v(row i of table(xy)),

vec$t2v(col 1 of table(b)),fvalue);

DELZ[I,J]=fvalue-xy[I,fcol];

Q[J]=Q[J]+DELZ[I,J]*R[I];

end;

Q[J]=100.*Q[J]/th[j];

b[k,1]=th[j];

end;

/* ----- Q=XT*R (STEEPEST DESCENT) ----- */

DO I=1 to NP;

DO J=1 to I;

total=0.0;

DO K=1 to Nobs;

total=total+DELZ[K,I]*DELZ[K,J];

end;

D[I,J]=10000.*total/(TH[I]*TH[J]);

D[J,I]=D[I,J];

end;

E[I]=max(1.0e-30,SQRT(D[I,I]));

/* IF(E[I] EQ 0.) E[I]=1.E-30; */

end;

/* ----- A IS THE SCALED MOMENT MATRIX ----- */

s50:

DO I=1 to NP;

/* if e[i] < 1.1e-30 then donext; */

DO J=1 to NP;

/* if e[j] < 1.1e-30 then donext; */

A[I,J]=D[I,J]/max(1.0e-30,E[I]*E[J]);

end;

P[I]=Q[I]/E[I];

PHI[I]=P[I];

A[I,I]=A[I,I]+GA;

end;

CALL invmat (A,NP,P);

/* ----- P/E IS THE CORRECTION VECTOR ----- */

stp=1.0;

s56:

DO i=1 to np;

j=iii[i];

TB[I]=P[I]*stp/E[I]+TH[I];

if empty(b[j,2]) then

if empty(b[j,3]) then donext;

else tb[i] = min(tb[i],b[j,3]);

else if empty(b[j,3]) then

tb[i] = max(tb[i],b[j,2]);

else tb[i] = min(b[j,3],max(tb[i],b[j,2]));

P[I]=(TB[I]-TH[I])*E[I]/stp;

end;

s60:

DO I=1 to NP;

IF(TH[I]*TB[I] LE 0.0) goto s66;

end;

do i=1 to np;

j=iii[i];

b[j,1]=tb[i];

end;

totalB=0.0;

DO I=1 to Nobs;

call func(vec$t2v(row i of table(xy)),

vec$t2v(col 1 of table(b)),fvalue);

xy[i,fcol]=fvalue;

R[I]=xY[I,ycol]-Fvalue;

totalB=totalB+R[I]*R[I];

end;

s66: total1=0.0;

total2=0.0;

total3=0.0;

DO I=1 to NP;

total1=total1+P[I]*PHI[I];

total2=total2+P[I]*P[I];

total3=total3+PHI[I]*PHI[I];

end;

IF(total2*total3 NE 0.0) then ARG0=total1/SQRT(total2*total3);

else arg0=0.0;

IF (ARG0 GT 1.0) THEN ARG0 = 1;

if (arg0 ne 0.0) then

begin;

ARG0=SQRT(1.-ARG0*ARG0)/arg0;

ANGLE=57.29578*ATAN(ARG0);

end;

else angle=90.0;

/* ---------- */

DO I=1 to NP;

IF(TH[I]*TB[I] le 0.0) goto s74;

end;

NTRIAL=NTRIAL+1;

IF(NTRIAL GT MAXTRY) GOTO s95;

IF (totalB/SSQ-1.0) le 0.0 then goto s80;

s74:

IF (ANGLE-30.0) GT 0.0 then goto s78;

stp=0.5*stp;

GOTO s56;

s78:

GA=10.*GA;

GOTO s50;

s80:

DO I=1 to NP;

TH[I]=TB[I];

end;

do i=1 to np;

IF(ABS(P[I]*stp/E[I])/(1.0e-20+ABS(TH[I])) > STOPCR) goto s94;

end;

GOTO s96;

s94: SSQ=totalB;

IF(NIT LT MIT) GOTO s34;

s95:

notes of table(b) = "No convergence in ".nit." iterations, "

.ntrial." trials.";

/* ----- END OF ITERATION LOOP ----- */

s96:

/* ----- WRITE CORRELATION MATRIX ----- */

do i = 1 to np;

do j = 1 to np;

a[i,j]=d[i,j];

end;

CALL invmat (a,NP,P);

DO I=1 to NP;

E[I]=SQRT(a[I,I]);

IF(E[I] EQ 0.0) E[I]=1.0E-30;

end;

IF(NP EQ 1) GOTO s104;

DO I=1 to NP;

DO J=1 to I-1;

A[J,I]=a[J,I]/max(1.0e-30,E[I]*E[J]);

a[i,j]=a[j,i];

end;

/* ----- CALCULATE REGRESSION COEFFICIENT ----- */

s104:

totalC=0.0;

totalC2=0.0;

DO I=1 to Nobs;

totalC=totalC+xY[I,ycol];

totalC2=totalC2+xY[I,ycol]*xY[I,ycol];

end;

totalc2=(totalC2-totalC*totalC/Nobs);

RSQ=1.-totalB/totalc2;

title of table(a) = 'Index of determination = '.rsq;

/* ----- CALCULATE 95% CONFIDENCE INTERVAL ----- */

Ndf=Nobs-Np;

Z=1./Ndf;

SDEV=SQRT(Z*totalB);

TVAR=1.96+Z*(2.3779+Z*(2.7135+Z*(3.187936+2.466666*Z**2)));

notes of table(a) = 'Standard Deviation = '.sdev.' with '.

ndf.' degrees of freedom.';

DO I=1 to NP;

j=iii[i];

b[j,4]=E[I]*SDEV;

b[j,5]=TH[I]/b[j,4];

TSEC=TVAR*b[j,4];

b[j,6]=TH[I]-TSEC;

b[j,7]=TH[I]+TSEC;

end;

s140:

/* ----- END OF PROBLEM RING BELL ----- */

/* DO I = 1 TO 10;

CALL $BELL;

END;

display table(b);*/

return;

end;