• Synopsis
• Copyright
• Procedure mnlfit
  • Original Procedure
  • Adaptations of Original Procedure
  • Definitions
  • Procedure mnlfit
• Logistic Cost Curve
  • Problem
  • Data
  • Fitting the Logistic Equation to the Data
  • Checks
• Standardised Logistic Cost Curve
  • Standardised Data and Fitting the Curve
  • Checking the Standardised Logistic Curve
  • Spreadsheet
• Observations

Logistic Cost Curve by the Levenberg-Marquardt Method

By J.M. R edwood

Synopsis

Procedure mnlfit , based on the Levenberg-Marquardt method, is used to fit a Logistic function to a project's cost data and provide a means for estimating cash flows in similar projects.

Copyright

The procedure used in this worksheet is an adaptation of Dr. D.E. Holmgren's procedure mnlfit 2000. His permission to reproduce and adapt the procedure is gratefully acknowledged.

Procedure mnlfit

Original Procedure

This worksheet implements the Levenberg-Marquardt method of fitting a non-linear function to a set of data. The general form of the fitting function is:

y = f(a[0] ...a[n] ;x[0] ...x[m] )

i.e., it is possible to have more than one independent variable. The user must supply a Maple routine to compute the function f (see example below). The procedure mnlfit will then compute the required derivatives. The user must also supply the data in the form of vectors of x and y values (or a matrix of x values in the case of more than one independent variable), as well as weights for the data points. The weights are usually the reciprocals of the standard deviations for each of the observed data points, but if not known, they may be set to unity. The routine below will iterate the solution either until the tolerance tol is reached or until the maximum number of iterations (50) is reached.

If the procedure requires a full 50 iterations, the data and formulation of the problem should be re-examined. In any case, the best way of learning how to use mnlfit is to run through the example below.

If there is sufficient interest and/or demand, a future version of this code will include automatic differentiation for more complicated model functions.

All of the usual qualifications concerning free software apply here (i.e., no claims of suitability for any particular application). Responsibility lies with the user to ensure that correct results are obtained!

Any questions about mnlfit should be addressed to D.E. Holmgren, together with any suggestions for modifications.

This procedure was developed jointly by:

D.E. Holmgren (holmgren@brandonu.ca) Brandon University, Brandon, MB, Canada
J.F. Ogilvie (ogilvie@munkebo.chem.ou.dk)
M. Monagan (monagan@cecm.sfu.ca) Center for Experimental and Constructive Mathematics, Simon Fraser University, Vancouver BC, Canada.

Adaptations of Original Procedure

Dr. D.E. Holmgren's original worksheet was downloaded from www.maplesoft.com/apps/categories/data_analysis_stats/data/html/genfit_6.html, and was adapted with Maple V Release 5.1. It can be used with both Maple V R5.1 and Maple 6.

The adaptations were:

(a) Removal of the redundant procedure fitstats .

(b) Removal of statements that attempted to calculate the covariance ratio, F (referred to as fstat in the original).

(c) Removal of all commands preceded by the # sign.

(d) Removal of the square root from weighting operations.

(e) Provision of means to control the output.

(f) Provision of an analysis of covariance and tests of the regression.

(g) Improved presentation of output.

Tests of the regression include G.W. Snedecor's test of the covariance ratio, F . This test is only valid for data drawn from a (multi-variable) Normal distribution. Users must satisfy themselves that this is the case, before placing any credence in this test. Naturally, the same applies to any other tests based, like Snedecor's, on Sir Ronald Fisher's z distribution and his transforms of z .

Definitions

mnlfit := proc(f, X, x::array, y::vector, wt::vector, P::{name,list(name)},p::vector,tol,Prin,digi::integer)

f is the name of the function to be fitted to the data,

X is the name of the independent variable appearing in this function,

x is the array of observed values of X ,

y is the array of observed values of the dependent variable or function,

wt is an array of weights (usually the inverses of the standard deviation of each value of dependent variable, but can be set to unity if weights are not known),

P is the name of the parameters of X as a vectorial quantity,

p is a vector of initial estimates of values of the parameters of X (rough guesses are usually good enough - the Levenberg-Marquardt method is not very sensitive to the initial estimates),

tol is a tolerance that serves as a criterion of convergence of the solution,
Prin controls the output of supplementary information - insert:

S in place of Prin if the fit statistics are required,

R if a table of residuals is required,

B if both statistics and residuals are required,

It if the residuals after each iteration are required,

All causes all the above items to be printed,

Inserting any other name, or leaving it as Prin , prevents their output,

digi is an integer that controls the number of digits printed by mnlfit .

mnlfit prints the number of iterations, the Std Deviation of Residuals (Goodness of Fit) statistic, the supplementary information required, and returns f .

Procedure mnlfit

> restart: Digits := 40:

> mnlfit := proc(f, X, x::array, y::vector, wt::vector, P::{name,list(name)},p::vector,
tol,Prin,digi::integer)
local gradf,nobs,npar,k,yck,xk,r,a,drv,niter,j,n,corr,corrprint,alpha,beta,cvm,b,p1,errors,dp,
res,nu,res1,sigma1,reso,nxvar,yc,finals,hdg,ssr,ssyc,ssy,cf,msr,msyc,msy,F,pF,Rsq,Tests, anova,resi,rout;
# Initialize.
if type(P,name) then
RETURN(mnlfit(f,X,x,y,wt,[seq(P[k],k=1..linalg[vectdim](p))],p,tol,Prin,digi));
fi;
gradf := linalg[grad](f,P);
reso := 1.0*10^9;
nobs := linalg[vectdim](y);
npar := linalg[vectdim](p);
nxvar := linalg[coldim](x);
xk := vector(nxvar);
r := vector(nobs);
a := matrix(nobs,npar,0.0);
corr := matrix(npar,npar,0.0);
drv := vector(npar);
p1 := vector(npar);
errors := vector(npar);
niter := 50;
nu := 1000;
# Begin iterations.
for n from 1 to niter do
# Compute initial fit and residuals ( = observed - calculated).
for k from 1 to nobs do
for j from 1 to nxvar do xk[j] := x[k,j]
od; yck,drv := gry(f,gradf,X,P,npar,xk,p);
r[k] := evalf(wt[k]*(y[k] - yck));
for j from 1 to npar do
a[k,j] := wt[k]*drv[j];
od:
od:
# Sum of squares of residuals
res1 := evalf(linalg[dotprod](r,r));
if Prin = It or Prin = All then print(` Iteration No `, n ,` initial SSR = `, evalf(res1,digi));
fi; # Compute matrix alpha of normal equations.
alpha := evalf(evalm(linalg[transpose](a) &* a));
# Compute vector beta for right sides of normal equations.
beta := evalf(evalm(linalg[transpose](a) &* r));
# Note use of factor nu in this next loop; this increases as one approaches solution.
for j from 1 to npar do
alpha[j,j] := evalf(alpha[j,j]*(1. + 1./nu));
od:
# Solve a system of normal equations to correct parameters.
dp := evalf(linalg[linsolve](alpha, beta));
# Compute sum of squares of residuals at new point p1 = p + dp.
p1 := evalf(evalm(p + dp));
for k from 1 to nobs do
for j from 1 to nxvar do xk[j] := x[k,j] od;
yck,drv := gry(f,gradf,X,P,npar,xk,p1);
r[k] := evalf(wt[k]*(y[k]-yck));
# Store yc[k] ie calculated y.
yc[k] := yck;
od:
res := evalf(linalg[dotprod](r,r));
if Prin = It or Prin = All then print(` new SSR = `, evalf(res,digi));
fi; # Modify nu according to whether the sum of squares increases or decreases after this correction. # For the wrong direction, elements are not corrected.
if res1 <= res then
nu := nu/10.;
else
# For the right direction, correct the elements.
nu := nu*10.;
p := evalf(evalm(p + dp));
fi:
# Test for convergence.
if abs(res - reso) <= tol then
break
fi;
reso := res;
od:
# Covariance matrix
cvm := evalf(linalg[inverse](alpha));
# Standard deviation of residuals
sigma1 := sqrt(res/(nobs - npar));
# Errors of parameters
for j from 1 to npar do
errors[j] := evalf(sigma1*sqrt(cvm[j,j]));
od:
for j from 1 to npar do
for k from 1 to npar do
corr[j,k] := cvm[j,k]/sqrt(cvm[j,j]*cvm[k,k])
od:
od: # Print results.
print (array([[`No of Iterations`,`Std Deviation of Residuals`],[n,evalf(sigma1, digi)]]));
print();
if Prin = S or Prin = B or Prin = All then
print();
# Calculate Sums of Squares, Covariances & Tests.
cf := (add(wt[j]*y[j],j=1..nobs))^2/nobs;
ssyc := add((wt[j]*yc[j])^2,j=1..nobs) - cf;
ssr := add(wt[j]^2*(yc[j] - y[j])^2,j=1..nobs);
ssy := add((wt[j]*y[j])^2,j=1..nobs) - cf;
msyc := ssyc/(npar-1);
msy := ssy/(nobs - 1);
msr := ssr/(nobs - npar);
Rsq := ssyc/ssy;
pF := 1 - stats[statevalf,cdf,fratio[(npar - 1),(nobs - npar)]](msyc/msr);
anova := array([[`Source`,Sigma(`Squares`),`DF`,`Mean Square`], [`Regression`,trunc(ssyc),(npar-1),evalf(msyc,digi)], [`Residuals`,trunc(ssr),(nobs-npar),evalf(msr,digi)],[`Total`,trunc(ssy),(nobs - 1),` `]]);
Tests := array([[R^2,`F`,`Prob of F by chance`],[` `,` `,`for Normal data`],
[evalf(Rsq,3),evalf(msyc/msr,3),evalf(pF,3)]]);
# Print Anova and Tests
print(`Analysis of Covariance`);
print();
print(anova);
print();
print(`Tests of Covariance`);
print();
print(Tests);
print();
print (`Final values of parameters`); hdg := array([`Parameter`,`Value`,`Standard Error`]); finals := linalg[augment]([seq(j,j=1..npar)], convert(p,list), [seq(errors[j],j=1..npar)]); finals := evalf(evalm(finals),digi):
finals := linalg[stackmatrix](hdg,finals);
print(finals);
print();
corr := evalf(evalm(corr),digi):
cvm := evalf(evalm(cvm),digi):
print(`Matrix of Covariances`);
print();
print(cvm);
print();
print(` Matrix of Correlation Coefficients `);
print();
print(corr);
print();
fi;
if Prin =R or Prin = B or Prin = All then print(` Table of residuals `);
print();
resi := convert([seq([(seq(evalf(x[j,k],digi),k=1..nxvar),evalf(y[j],digi),evalf(yc[j],digi),
evalf(y[j]-yc[j],digi), evalf(wt[j]*(y[j]-yc[j])/sigma1, digi))],j=1..nobs)],matrix);
hdg := array([(seq(`x obs`[k],k=1..nxvar),` y obs `,` y calc `,` difference `,`wt*difference / sd of residuals`)]);
resi := linalg[stackmatrix](hdg,resi);
print(resi);
print();
fi;
RETURN(subs(seq(P[j]=evalf(p[j],digi),j=1..npar),f));
end:

gry := proc(f,G,x,a,n,xval,p::vector)
# Procedure to evaluate derivatives and formula
local k,y,i,drv,point;
point := {x=xval, seq(a[i]=p[i], i=1..n)};
drv := vector(n);
for k from 1 to n do
drv[k] := evalf(eval(G[k],point))
od:
y := evalf(eval(f,point));
RETURN(y,drv);
end:

Logistic Cost Curve

Problem

A method was needed for estimating the cash flows of engineering development projects undertaken by a certain company. One such project took 27 months to complete and the cumulative cost returns were collected throughout the life of the project. The accounts were closed 3 months after completion of the project, when the last bills were brought to account.

The cumulative costs were collected at the end of each month up to the final fixed price of 1,000,000. The data therefore comprises 30 pairs of (end) of month numbers and the cumulative costs in k.

The cumulative logistic distribution function often fits data from growth situations that are limited by a finite resource. In this case, the costs are limited by the fixed price for the job. They grow slowly at first as just a few, then more and more people on the development team become involved. They then increase more rapidly as parts are bought in and manufacturing, assembly and test proceed, and then taper off as the manufacturing and development teams reduce with final evaluation and delivery to the customer, followed by settlement of the last bills from suppliers.

Data

The data recorded was

> cost := [[1,2],[2,7],[3,11],[4,20],[5,36],[6,58],[7,87],[8,106],[9,138],[10,186],[11,234], [12,268],[13,352],[14,396],[15,489],[16,545],[17,643],[18,672],[19,750],[20,819], [21,854],[22,875],[23,924],[24,938],[25,948],[26,968],[27,975],[28,982],[29,995],[30,1000]]:

n := nops(cost):

As a first step, the data is plotted.

> with(plots):
pointplot(cost,view=[0..30,0..1000],labels=[`Months`,`Cost k`],title=`An Engineering Development Project`);

Warning, the name changecoords has been redefined

The shape of the plot is typical of the cumulative expenditure for a fixed price project, and the sigmoidal form suggests that the logistic equation should fit the data. (Note the familiar slow start because engineers and draughtsmen were still involved with a different project!)

The (cumulative) logistic distribution function is given by where are constants, and t is the independent variable.

Fitting the Logistic Equation to the Data

Starting values for are chosen from previous experience, and all weights are set to unity.

> p := vector(3,[1000,5,1]);
w := vector(n,1):

> f := a[1]/(1+exp(a[2]-a[3]*t[1]));

Choosing Prin = B and Digi = 6 , mnlfit is run.

> costv := linalg[col](cost,2):
month := matrix(n,1,linalg[col](cost,1)):
f1 := mnlfit(f,t,month,costv,w,a,p,10^(-20),B,6);

> F1:=unapply(subs(t[1]=t,f1),t);

F1 is the desired logistic function that fits the data. It is now plotted against the data from which it was obtained.

> aa := pointplot(cost,view=[0..30,0..1000],labels=[`Months`,`Cost k`],title=`An Engineering Development Project`): bb := plot(F1(x),x=0..n):
display(aa,bb);

Checks

By eye, the fit appears very good and this is supported by the relatively small Std Deviation of Residuals statistic.

However, one should plot the residuals to check that there is no clear trend,

> pts := zip((x,y)->[x,y],[seq(j,j=1..n)],[seq(cost[i,2]-F1(i),i=1..n)]):
pointplot(pts,labels=["Month","Actual - Est Cost k"]);

The above plot shows the residuals are scattered at random, and the Table of Residuals shows they are all less than about twice the Std Deviation of Residuals . There is no reason to believe that the cost data, or the residuals, should be Normally distributed. Snedecor's test of F produced by mnlfit is therefore meaningless. However, the test shows that the regression accounts for 99.9% of the total covariance and F is very large. These reinforce the impression obtained by eye, and that given by the small Std Deviation of Residuals , that the fit is very good.

Standardised Logistic Cost Curve

By standardising the data as percentages of the total time and cost of the project, the resulting function would be more useful to accountants estimating the cash flow of other similar, planned projects.

The usual financial health warnings about past performance being no guide to future performance apply with a vengeance!!

Standardised Data and Fitting the Curve

> costvpc := evalm(100*costv/costv[n]):

> timepc := evalm(100*month/month[n,1]):

> costpc := zip((x,y)->[x,y],convert(timepc,vector),evalm(costvpc)):

> f100 := mnlfit(f,t,timepc,costvpc,w,a,p,10^(-20),S,6);

> F100:=unapply(subs(t[1]=t,f100),t);

> aa := pointplot(costpc,labels=[`% Total Time`,`% Total Cost`],title=`Fixed Price Engineering Development Project Estimates`):
bb := plot(F100(x),x=0..100):
with(plots): display(aa,bb);

Checking the Standardised Logistic Curve

Checking, as before, that there is no clear trend in the residuals,

> pts := zip((x,y)->[x,y],[seq(j,j=1..30)],[seq((costvpc[i]-F100(i*100/30)),i=1..30)]): pointplot(pts,labels=["%Total Time x 30/100","%Actual - %Est Cost "]);

As expected, the plot of errors is identical to the previous one.

Finally, a sample spreadsheet is produced for the benefit of the accountants who have to estimate the cash flows of similar projects in future.

Spreadsheet

The missing 3402 from the Spent to Date after 24 months is reserved for the Project Manager's celebrations...hic, haec, hoc!

>

>

Observations

a. Fitting the very non-linear logistic (cumulative) distribution curve to a set of data is widely acknowledged to be extremely demanding. Nevertheless, Maple and Dr. D.E. Holmgren's procedure mnlfit produced an excellent fit to the cost data in this example. The author has obtained similarly good results with other sets of data, e.g. population of the USA.

b. It is quite remarkable that a procedure based on user-programmable symbolic logic should be able to perform the kind of computations normally done with specially coded numerical programs.

c. The very ease with which mnlfit produces results should warn us that care is essential when interpreting, or using, the results. In this example, applying the time-cost function to projects that are not very similar to the one from which the data was obtained, is a recipe for disaster.

d. The test of covariance shows that the regression accounts for 99.9% of the total covariation, indicating an excellent fit, and the value of F indicates the same. These reinforce the goodness of fit indicated by the Std Deviation of Residuals .

Disclaimer : While every effort has been made to validate the solutions in this worksheet, Waterloo Maple Inc. and the contributors are not responsible for any errors contained and are not liable for any damages resulting from the use of this material.