Exponential Curve Fitting

Introduction

This new article describes the exponential curve fitting method implemented in Graphics-Explorer,
my equations grapher program.
It replaces the old article, which can be found [here].
New is an exerciser program allowing step by step observation of the curve fitting process.

The curve fitter calculates the best fitting exponential function given a set of points.
This function is
y = a.bx + c
where a,b,c are called the parameters.

Contents

Exerciser description

Below is a reduced picture of the exerciser at work
 a,b,c calculated parameters. (black:valid; red:invalid) preset params click radio button to preset parameters preset x1..x5 click radio button to preset points x, y list of points, maximum is 10 [OK] validate click to check/organize points prior to calculation [GO] start calculation of parameters a,b,c [ ] testmode mark to watch step by step calculation

Operation

1. selecting 1 of 4 preset point sets (click radio button)
2. mouseclick on the coordinate system (2nd click removes point)
3. typing into the [x,y] points list
After changing points in any way, click the [OK] button to check and to sort the points.
To watch the process step by step, mark the [testmode] checkbox.
Then click the [GO] button to calculate the a,b,c parameters and paint the graph.
To clear parameters or points, click the appropriate [X] button or [RESTART].

Curve fitting theory

Input to the curve fitter is a set of points [x1,y1]....[xn,yn]
The minimal required number of points is 3.
The maximum number of points is 10.
The fitter calculates parameters a,b,c such that the curve y = a.bx + c has the smallest distance to these points.

For internal computations y = a.ebx + c is used (eb is reported for b)
The method has a numerical and an analytical part.
These are the steps:
1. obtain an approximation of b
2. using b, obtain an approximation of a
3. using a,b, obtain an approximation of c
4. discard a,b and replace y = aebx + c ....by.... y - c = aebx
5. ln(y - c) = ln(a) + bx ...........a linear function
6. call the least square method to find a,b of the best fitting line

approximation of b

array dx holds the x differences between adjacent points
dx1 = x2 - x1....dx3 = x4 - x3

array dy holds the y differences
dy1 = y2 - y1....dy3 = y4 - y3

array cx holds the centers between points
cx1 = (x1+x2)/2....cx3 = (x3+x4)/2

array dq holds the differential quotients
dq1 = dy1/dx1....dq3 = dy3/dx3

Note: in the Delphi language, cx1 is written as cx[1]...etc.

Starting with
y = aebx + c..........differentiation
y' = abebx
which eliminates c.
So, we may write
y1' = abebx1
y2' = abebx2
y3' = abebx3
Parameter a can be eliminated by division of subsequent y':

with above arrays in mind:

The center values (cx) are used, because the tangent (differential quotient dq) is most accurate
between two x values.
For subsequent points, the average value of b is calculated.
This is the approximation of b.

approximation of a

Starting with
y1 = aebx1 + c
y2 = aebx2 + c
......
we may eliminate c by subtraction of subsequent rows:
dy1 = y2 - y1 = a(ebx2 - ebx1)
so....

An average value of a is calculated.
This is the approximation of a.

approximation of c

y = aebx + c .....so
c = y - aebx
We may write
c = y1 - aebx1
c = y2 - aebx2
....
An avarage value of c is calculated.
This is the approximation of c.

Intermezzo

In equation
y = gx
the value of y will increase by a factor g if x increases by 1.
g is called the growth factor.

Exponential functions (having c=0, the x-axis is the horizontal asymptote), are the result
of constant relative growth.

In equation
y = ebx
we notice that g = eb.
b is called the relative growth speed because

So, if a function has a constant relative growth speed of b, it's growth factor g = eb
where e is the base of the natural logarithm. (e = 2.71828....)

Before, we have obtained an approximation of c so we may write ...y - c = aebx
If the points (x1,y1)...(xn,yn) are on a true exponential function,
we will find the same growth factor, comparing subsequent points.
In the next part, fine tuning of c, we are going to vary c to obtain the smallest
deviation in growth factors.

fine tuning c

At this point we forget the a and b values and continue with c and the points [x1,y1]...[xn,yn].

For a certain value of c we write
y1-c = aebx1
y2-c = aebx2
.....
division eliminates a

Two successive points produce a value of ee.
The values of ee are summed in variable AVG, for now AVG = ee1 + ee2 + ....
The growthfactors are stored in array growthfactor:
growthfactor1 = eee/dx1..........

By summing the ee values, we actually sum values b.dx1 + b.dx2 + b.dx3... = b(last X - first X)
So, b = AVG / (last X - first X)
Average growth factor is eb
AVG = eAVG/(lastX-firstX)
Next, we have to compare the values in the growthfactor array with AVG.
Deviation = abs(AVG-growthFactor1) + abs(AVG-growthfactor2) + ......
Deviation is the result of procedure deviation which is called for a certain value of c.

least square method

The least-square-method is an algorithm to calculate the best fitting polynomial y = c0 + c1x + c2x2+ ....
given a set of points.
With a final value of c, we may write

which represents a linear function like y = a + bx.
We present the points (ln(yi-c),xi) to the least square procedures to obtain the best a,b values.
Resulting a is replaced by ea to cancel the ln( ) function.

The least square method is a nice application of linear algebra.
Look [here] for a description.

Aborts

The deviation procedure produces an abort if c is between two y values.
In that case y = c intersects the curve.
c must be outside the range y1, y2,y3.....
ee = (y2 - c)/(y1 - c)
abort = ee <= 0
A negative or zero ee value does not allow for an ln function.
Below is the flowchart for the fine tuning of c
DeltaC is the increment value of c.
The fine tuning ends when deltaC has become very small.

Watchdog timer

Procedure morepointsexp controls the calculation of the a,b,c parameters.
Floating point operation are between try...except directives.
There are 2 cases for an abort:
1.Illegal Floating point operation
This situation occurs if the points are too far off to make an exponential function
an example being (1,1) , (5,5) , (8,2)

2.C runaway
If C is too far out, an even further position may yield a smaller growth factor deviation.
So, C "runs away".
The watchdog counter w ends this process.
Fine tuning of c must be completed in maximal 100 steps, otherwise the process is aborted.
This situation occurs in case of sharp bends.

The Delphi-7 project

There are 4 units and 1 form:
form1,unit1: paintbox, events, controls
unit2: data and procedures for the curve fitting algorithm
unit3: procedures for testmode reporting
unit4: least square procedures

Unit1

Graphs are painted in bitmap map, points are painted on paintbox1.
The dimension is 800*800 pixels, coordinate (0,0) is at map position (400,400).
x domain is -10< x < 10,.........y range is -10< y <10.
The scale is 40:1.
There are 2 important flags:
pointsOK : all points in the list are OK and sorted in ascending order.
paramsOK : parameters are valid
If paramsOK = true, pressing [GO] just plots the function.
If paramsOK = false, pressing [GO] calculates a,b,c and then plots the function.
if parameters a,b,c do not represent the actual situation, they are painted in red.

Please refer to the source code for details.

Unit2

Main call is made to procedure makeExpFunction.
Points are collected from unit1 and stored in array Epoints[ ] of type Tpnt.
```type TPnt = record
x,y : double;
end;
```

procedure deviation calculates the variation of growthfactors while fine tuning c.
These values are stored in variables Bvar1,Bvar2.

Unit3

Unit2 calls procecures in unit3 to display the steps in the curve fitting process.
Information is written in form1 paintbox1.
The waitflag : boolean (in unit1) controls the step-by-step action.
While the waitflag is true, the program executes the processmessages procedure.

Unit4

Has the necessary procedures for the least square method.
Main call is to procedure buildexp12(var pp : Tpoints; maxp : byte),
where Tpoints is the array of points and maxp the number of points supplied.

This concludes the description.