Curve Fitting 2
 Volume Number: 1 Issue Number: 12 Column Tag: Forth Forum

# Curve Fitting, Part II

By Jörg Langowski, Chemical Engineer, Grenoble, France, MacTutor Editorial Board

After last month's refresher on accessing floating point routines from Forth, this column is now going to show you the main part of the curve fitter program.

Let's state our objective again: We have a series of data points (measurements) yi at certain time points ti. The measured data is supposed to follow some quantitative law, so that we can state a theoretical relationship between time t and data y:

The ai are parameters that determine the exact form of the function f.

Lets assume we have estimated initial values for the parameters ai° somehow so that they are not too far away from the true values. Then we need a method that creates some correction terms ai, which when added to the initial values give new, better estimates of the parameters ai':

As you could read in last month's column, the ai are eventually obtained from the solution of a system of linear equations, and we had defined the word gauss to implement the Gauss algorithm that solves such a system.

This month I am going to show you how one sets up the equations (example given for three parameters),

starting from the data and the initial estimate of the function that one wants to fit to it. The full details of the method are given in an appendix to this column.

At this point it is only important to know that the coefficients cij, as defined in last month's column, look like the following:

where N is the total number of data points and fk/ai is the first derivative of the theoretical function at the time tk with respect to the parameter ai.

The terms on the right hand side of the equations, bi, are:

Ri are the residuals (as defined last month),

the differences between the theoretical and measured values at the time points ti.

Therefore, the curve fitting algorithm will consist of the following major parts:

- one routine that calculates the theoretical function value f(ti), given a set of parameters (this will be the 'model' that you use to fit your data),

- one routine that calculates the derivative of this function with respect to one of the parameters ai,

- calculation of all the derivatives fk/ai with k ranging from 0 to (# of data points - 1) and i from 0 to (# of parameters - 1) and saving those derivatives in a matrix,

- computation of the coefficients cij and bi (and setting them up in a matrix),

- solution of the linear equation system thus obtained, giving correction values for the parameters ai,

- changing the parameters by the correction terms and repetition of the whole algorithm if the change is still larger than some predefined (small) number.

The algorithm is implemented in this month's program example (listing 1). The first part contains some additions to the floating point and Gauss algorithm routines that you've already seen last month. One bug had to be corrected in gauss, which left a number on the stack that was supposed to be dropped (added drop in italics in the listing). The floating point output routine has been slightly modified so that real numbers are always printed with one digit in front of the decimal.

In order to be able to check whether the iteration has been completed or not, we will have to compare real numbers, so first we define the

## Floating point comparison operator

The SANE package provides a routine for the comparison of two floating point numbers. Four cases are distinguished and the processor status register flags set accordingly:

X N Z V C

x < y 1 1 0 0 1

x = y 0 0 1 0 0

x > y 0 0 0 0 0

unordered0 0 0 1 0

where the address of x is on top of the address of y on the stack.

The operation code for floating point comparison in SANE is 8, so we define

```hex : f> 8 fp68k @sr 10 and ;
```

for comparing two real numbers on the stack. Bit 5 of the result (the X flag) is and-ed out, the reason being that this is the only one that is left unchanged after the SANE trap has been called and the Forth interpreter has taken over again. All the other flags cannot be used, so we won't be able to compare two real numbers for being equal. With 80-bit precision, however, such a test does not make much sense, you will almost never check for equality but rather if a number is less than one other or whether their difference is smaller than some predefined value.

The comparison for 'unorderedness', which sets the V flag, is also trashed by the Forth interpreter. This would only be important if we wanted to test for a NaN (Not a Number), which for instance results from a 0/0 division.

This leaves us with the X flag as the only one usable after the fp68k call. It is set when the real number whose address is on top of the stack is less than the number whose address is below it.

## Random number generator - floating point version

For initializing the 'data points' as an input to the curve fitting program, we are going to use simulated data from the same function that we are fitting. We add some 'random noise' to it and need a random number generator for that. SANE contains a floating point random number generator, which computes from an 80-bit floating point input a new random 80-bit floating point number in the range 0 < x < 231-1. We want a number between 0 and 1 (for convenience) so the word ranf scales this number after calling the random number routine. ranf leaves the address of the result on the stack. For initializing the random number generator with an arbitrary seed, ranset is defined.

You will notice one important concept in the definition of ranf and ranset: both words operate on variables (rr, rs and sf) that are axed after the definition is completed. The variables are erased from the vocabulary that way and stay local to the random number routines. This is the way the 'nameless' tokens are created that you read about some issues ago.

The function to be fitted is defined as the word func. 50 bytes are reserved for the array par, so that at maximum 5 parameters (extended precision) can be used in the function. The data points are stored in the single precision arrays xdat and ydat. init_pars initializes the parameters to some arbitrary values, in this case par[1] = 1.0, par[2] = 2.0, and par[3] = -0.1. Calling func with these values, init then creates a set of simulated data (including random noise).

deriv calculates the first derivatives of the function with respect to the parameters from differences,

This definition, too, uses the local variables da1, da2, da3 and da4, which are axed afterwards.

make_derivmat computes all necessary derivative values and sets them up in a single precision matrix. This matrix is then used by the word make_resmat to compute the sums that make up the coefficient matrix. For the right-hand-side coefficients one also needs the residuals, which are calculated and stored into a matrix by residuals.

One iteration of the curve fitting process is done by the word one_iter, which sets up the derivative matrix and computes the residuals by calling the appropriate words, prints the sum of error squares (so you can check the quality of the fit) then sets up the coefficient matrix and calls gauss to solve the linear equation system. The solution is stored in the array delta; these are the correction terms that have to be added to the parameters to get improved estimates. new_pars does this correction and prints out the values of the new parameters, leaving false (zero) on the stack if any parameter has changed by more than one part in 10-5.

The actual curve fitter, nlsqfit, then loops through the iteration until the best fit is obtained.

You can check the fitting process by calling init first and then, for instance, setting par[1] = 2.0, par[2] = 1.0, par[3] = -0.05:

```two par x2x  one par 10+ x2x
two par 20 + f/
```

and start nlsqfit. This will bring you back close to the simulated values, par[1] = 1.0, par[2] = 2.0, and par[3] = -0.1 in about 4 to 5 iterations. The fitted values will not be exactly equal to the simulated ones because of the random noise added to the data.

So far, this is only a skeleton of a curve fitting program because we cannot input floating point numbers manually or from the clipboard (e.g. it would be nice to transfer data from a spreadsheet); also a graphical output will be needed to display the data points and the fitted curve. Next column will deal with those problems.

```Listing 1: Non-linear least squares curve fitting routine
( © 1985 J.Langowski by MacTutor )
( Note that this is not stand-alone but needs some definitions from last
month's example. Only the changed parts are printed here )

hex
: f> 8 fp68k @sr 10 and ;  : fabs f fp68k ;
: lnx 0 elems68k ; : log2x 2 elems68k ;
: ln1x 4 elems68k ; : log21x 6 elems68k ;
: expx 8 elems68k ; : exp2x a elems68k ;
: exp1x c elems68k ; : exp21x e elems68k ;
: x^i 8010 elems68k ; : x^y 8012 elems68k ;
: compoundx c014 elems68k ;
: annuityx c016 elems68k ;
: sinx 18 elems68k ; : cosx 1a elems68k ;
: tanx 1c elems68k ; : atanx 1e elems68k ;
: randomx 20 elems68k ;   decimal

: dec. ( float\format# -- )
zzformat ! zzformat swap zzs1 b2d
zzs1 dup w@ 255 > if ." -" else ."  " then
dup 4+ count over 1 type ." ."
swap 1+ swap 1- type ( mantissa )
2+ w@ ( get exponent )
1 w* zzformat @ + 1-
." E" . ;

( define constants )
float one  float -one  float zero  float two  float four
1 sp@ one in2x drop  -1 sp@ -one in2x drop
0 sp@ zero in2x drop
2 sp@ two in2x drop  4 sp@ four in2x drop
( define some floating accumulators)
float fa1   float fa2   float fa3   float fa4

( Gauss algorithm for linear equations)
float dg    float fk    float ee
variable nv   variable coeff variable solution
( addresses for storing actual parameters)
: gauss ( z\x\n | --)  nv !  8- coeff !  solution !
nv @ 1- 0 do  ( i-loop)
i dup coeff @ calc.offset dg s2x ( diag elem)
nv @ i 1+ do  ( j-loop)
i j coeff @ calc.offset fk s2x   dg fk f/
nv @ 1+ j do  ( k-loop)
k i coeff @ calc.offset fa1 s2x
fk fa1 f*  fa1 fneg  ( -fk*x[i,k])
j i coeff @ calc.offset dup fa1 s+
fa1 swap x2s
loop
loop
loop
nv @ dup 0 do i over coeff @ calc.offset  fa1 s2x
fa1 solution @ i 4* + x2s loop drop
1 nv @ 1- do
i dup coeff @ calc.offset dg s2x
solution @ i 4* + ee s2x  dg ee f/
0 i 1- do i j coeff @ calc.offset fa1 s2x
ee fa1 f* fa1 fneg
solution @ i 4* + dup fa1 s+ fa1 swap x2s
-1 +loop
-1 +loop
nv @ 0 do  solution @ i 4* +  fa1 s2x
i dup coeff @ calc.offset  fa1 s/
fa1 solution @ i 4* + x2s
loop ;

( declarations for curve fitter )
create ydat 400 allot   create xdat 400 allot
create residues 400 allot
100 10 matrix derivmat    10 11 matrix resmat
3 constant npars        10 constant npts
create par 50 allot   create delta 20 allot
float eps  float errsum
1 sp@ eps in2x drop  10000 sp@ eps in/ drop
float onehundred  100 sp@ onehundred in2x drop
float ten  10 sp@ ten in2x drop
( define function )
: func ( x -- f[x] = par[1] + par[2] * exp[par[3]*x] )
par 20 + over f* dup expx
par 10 + over f*  par over f+ ;
: test 10 0 do i sp@ fa1 in2x . 2 spaces
fa1 func 10 dec. cr loop ;
: >fa1  fa1 s2x ;
: init_pars
one par x2x  two par 10+ x2x
-one par 20 + x2x  ten par 20 + f/ ;
init_pars

( derivative, matrix of derivs )
float da1 float da2 float da3 float da4  ( local vars )
: deriv ( par \ x -- d-func/d-par at x )
dup da1 x2x da2 x2x  dup da4 x2x  eps da4 f*
da4 over f+          da2 func  da3 x2x
da4 over 2dup f- f-  da1 func  da3 f-
da4 da3 f/   two da3 f/  da4 swap f+  da3  ;
axe da1 axe da2 axe da3 axe da4

: make_derivmat
npts 0 do  npars 0 do
xdat j 4* +  >fa1
par i 10 * +   fa1 deriv  j i derivmat x2s
loop  loop ;

( calculate residuals )
: residuals
zero errsum x2x
npts 0 do
xdat i 4* + >fa1   fa1 func   ydat i 4* +  swap  s-
fa1 residues i 4* + x2s  fa1 dup f*  fa1 errsum f+
loop  ;
: .resid
npts 0 do  residues i 4* + >fa1  fa1 7 dec. cr loop ;

( make matrix of residuals )
make_resmat
npars 0 do   npars 0 do    zero fa1 x2x
npts 0 do
i k derivmat fa2 s2x    i j derivmat fa2 s*
fa2 fa1 f+  loop
fa1  i j resmat x2s    fa1  j i resmat x2s
loop  loop
npars 0 do  zero fa1 x2x
npts 0 do
i j derivmat fa2 s2x  residues i 4* + fa2 s*
fa2 fa1 f-  loop
fa1  i npars  resmat x2s loop ;

( calculate correction terms)
: one_iter
make_derivmat    residuals
." sum of error squares: " errsum 7 dec. cr
make_resmat      delta 0 0 resmat npars gauss ;

: new_pars 16 ( true if no significant changes )
npars 0 do par i 10 * +
delta i 4* +  over  s+
." par[" i . ." ] = " dup  7 dec. cr
delta i 4* + fa1 s2x  fa1 f/
fa1 fabs  eps fa1 f> and loop ;

( ranf, initialize data matrices )
float rr float rs float sf ( local to ranf )
1 31 scale 1 - sp@ sf in2x drop
: ranset rr x2x ;
: ranf rr randomx  rr rs x2x  sf rs f/  rs ;
axe rr axe rs axe sf
12345678 sp@ fa1 in2x drop fa1 ranset
80 ' npts !
: init npts 0 do i sp@ fa1 in2x 4*
xdat over + fa1 swap x2s
ydat over + fa1 func  ranf fa2 x2x
ten fa2 f/  fa2 over f+  swap x2s
i .  xdat over + >fa1 fa1 7 dec. 2 spaces
ydat  +   >fa1   fa1 7 dec. cr    loop ;

( print matrices for debugging )
: .dmat
npts 0 do
npars 0 do  j i derivmat >fa1 fa1 5 dec. loop
cr loop ;
: .rmat
npars 0 do
npars 1+ 0 do  j i resmat >fa1 fa1 5 dec. loop
cr loop ;

( nonlinear fit, core routine)
: nlsqfit cr   begin  one_iter cr  new_pars cr  until ;
```

## Appendix: Theoretical background of the curve fitting routine

We want to determine the values of the ai in such a way that the differences between the theoretical function and the measured yk values at times tk become a minimum. These differences are called the residuals rk:

rk = f (tk, a1, a2, a3, .... , an) - yk

and one usually tries to minimize the sum of the squared residuals of all data points.

Lets assume ri are the 'true' residuals that one obtains with the exact ai values. If we estimate the parameters by some initial values ai°, then 'computed' residuals

Rk = f (tk, a1°, a2°, a3°, .... , an°) - yk

can be calculated, which are usually larger than the true ones. To get a correction term that brings the ai° closer to the 'true' ai, one now linearly expands the function f around the estimated value:

f(tk,a1,a2,....,an)

f(tk,a1°,a2°,....,an°) + fk/a1(a1-a1°)

+ fk/a2 (a2-a2°)

. . . . .

+ fk/an(an-an°)

The differences, (ai-ai°), are denoted by ai, now we can write

```f(tk,a1,a2,....,an) - yk
f(tk,a1°,a2°,....,an°) - yk
+  fk/a1 a1 +  fk/a2 a2
. . . . .
+  fk/an an
```

which gives us a relationship between the true and the computed residuals

```rk   Rk +  fk/a1 a1 +  fk/a2 a2
. . . . .
+  fk/an an .
```

It is the sum of the squares of the true residuals (N being the number of data points)

that has to be minimized with respect to changes in ai, this means all the derivatives Q/(ai) have to be zero simultaneously. When you evaluate the expressions for the Q/(ai) and set them to zero, you arrive at the equation system that was desribed in the main article.

Community Search:
MacTech Search:

Microsoft Office 2016 16.11 - Popular pr...
Microsoft Office 2016 - Unmistakably Office, designed for Mac. The new versions of Word, Excel, PowerPoint, Outlook, and OneNote provide the best of both worlds for Mac users - the familiar Office... Read more
Adobe Photoshop CC 2018 19.1.2 - Profess...
Photoshop CC 2018 is available as part of Adobe Creative Cloud for as little as \$19.99/month (or \$9.99/month if you're a previous Photoshop customer). Adobe Photoshop CC 2018, the industry standard... Read more
Adobe Dreamweaver CC 2018 18.1.0.10155 -...
Dreamweaver CC 2018 is available as part of Adobe Creative Cloud for as little as \$19.99/month (or \$9.99/month if you're a previous Dreamweaver customer). Adobe Dreamweaver CC 2018 allows you to... Read more
Adobe Flash Player 29.0.0.113 - Plug-in...
Adobe Flash Player is a cross-platform, browser-based application runtime that provides uncompromised viewing of expressive applications, content, and videos across browsers and operating systems.... Read more
Drive Genius 5.2.0 - \$79.00
Drive Genius features a comprehensive Malware Scan. Automate your malware protection. Protect your investment from any threat. The Malware Scan is part of the automated DrivePulse utility. DrivePulse... Read more
MegaSeg 6.0.6 - Professional DJ and radi...
MegaSeg is a complete solution for pro audio/video DJ mixing, radio automation, and music scheduling with rock-solid performance and an easy-to-use design. Mix with visual waveforms and Magic... Read more
ffWorks 1.0.7 - Convert multimedia files...
ffWorks (was iFFmpeg), focused on simplicity, brings a fresh approach to the use of FFmpeg, allowing you to create ultra-high-quality movies without the need to write a single line of code on the... Read more
Dash 4.1.5 - Instant search and offline...
Dash is an API documentation browser and code snippet manager. Dash helps you store snippets of code, as well as instantly search and browse documentation for almost any API you might use (for a full... Read more
Evernote 7.0.3 - Create searchable notes...
Evernote allows you to easily capture information in any environment using whatever device or platform you find most convenient, and makes this information accessible and searchable at anytime, from... Read more
jAlbum Pro 15.3 - Organize your digital...
jAlbum Pro has all the features you love in jAlbum, but comes with a commercial license. You can create gorgeous custom photo galleries for the Web without writing a line of code! Beginner-friendly... Read more

## Latest Forum Discussions

Around the Empire: What have you missed...
Oh hi nice reader, and thanks for popping in to check out our weekly round-up of all the stuff that you might have missed across the Steel Media network. Yeah, that's right, it's a big ol' network. Obviously 148Apps is the best, but there are some... | Read more »
All the best games on sale for iPhone an...
It might not have been the greatest week for new releases on the App Store, but don't let that get you down, because there are some truly incredible games on sale for iPhone and iPad right now. Seriously, you could buy anything on this list and I... | Read more »
Everything You Need to Know About The Fo...
In just over a week, Epic Games has made a flurry of announcements. First, they revealed that Fortnite—their ultra-popular PUBG competitor—is coming to mobile. This was followed by brief sign-up period for interested beta testers before sending out... | Read more »
The best games that came out for iPhone...
It's not been the best week for games on the App Store. There are a few decent ones here and there, but nothing that's really going to make you throw down what you're doing and run to the nearest WiFi hotspot in order to download it. That's not to... | Read more »
Death Coming (Games)
Death Coming 1.1.1.536 Device: iOS Universal Category: Games Price: \$1.99, Version: 1.1.1.536 (iTunes) Description: --- Background Story ---You Died. Pure and simple, but death was not the end. You have become an agent of Death: a... | Read more »
Hints, tips, and tricks for Empires and...
Empires and Puzzles is a slick match-stuff RPG that mixes in a bunch of city-building aspects to keep things fresh. And it's currently the Game of the Day over on the App Store. So, if you're picking it up for the first time today, we thought it'd... | Read more »
What You Need to Know About Sam Barlow’s...
Sam Barlow’s follow up to Her Story is #WarGames, an interactive video series that reimagines the 1983 film WarGames in a more present day context. It’s not exactly a game, but it’s definitely still interesting. Here are the top things you should... | Read more »
Pixel Plex Guide - How to Build Better T...
Pixel Plex is the latest city builder that has come to the App Store, and it takes a pretty different tact than the ones that came before it. Instead of being in charge of your own city by yourself, you have to work together with other players to... | Read more »
Fortnite Will Be Better Than PUBG on Mob...
Before last week, if you asked me which game I prefer between Fortnite Battle Royale and PlayerUnknown’s Battlegrounds (PUBG), I’d choose the latter just about 100% of the time. Now that we know that both games are primed to hit our mobile screens... | Read more »
Siege of Dragonspear (Games)
Siege of Dragonspear 2.5.12 Device: iOS Universal Category: Games Price: \$9.99, Version: 2.5.12 (iTunes) Description: Experience the Siege of Dragonspear, an epic Baldur’s Gate tale, filled with with intrigue, magic, and monsters.... | Read more »

## Price Scanner via MacPrices.net

Sunday Sales: \$200 off 13″ Touch Bar MacBook...
Amazon has new 2017 13″ 3.1GHz Touch Bar MacBook Pros on sale this weekend for \$200 off MSRP, each including free shipping: – 13″ 3.1GHz/256GB Space Gray MacBook Pro (MPXV2LL/A): \$1599.99 \$200 off... Read more
B&H drops prices on 15″ MacBook Pros up t...
B&H Photo has dropped prices on new 2017 15″ MacBook Pros, now up to \$300 off MSRP and matching Adorama’s price drop yesterday. Shipping is free, and B&H charges sales tax for NY & NJ... Read more
Apple restocks Certified Refurbished 2017 13″...
Apple has restocked Certified Refurbished 2017 13″ 2.3GHz MacBook Pros for \$200-\$230 off MSRP. A standard Apple one-year warranty is included with each MacBook, models receive new outer cases, and... Read more
13″ Space Gray Touch Bar MacBook Pros on sale...
Adorama has new 2017 13″ Space Gray Touch Bar MacBook Pros on sale for \$150 off MSRP. Shipping is free, and Adorama charges sales tax in NY & NJ only: – 13″ 3.1GHz/256GB Space Gray MacBook Pro (... Read more
Best deal of the year on 15″ Apple MacBook Pr...
Adorama has New 2017 15″ MacBook Pros on sale for up to \$300 off MSRP. Shipping is free, and Adorama charges sales tax in NJ and NY only: – 15″ 2.8GHz Touch Bar MacBook Pro Space Gray (MPTR2LL/A): \$... Read more
Save \$100-\$150+ on 13″ Touch Bar MacBook Pros...
B&H Photo has 13″ Touch Bar MacBook Pros on sale for \$100-\$150 off MSRP. Shipping is free, and B&H charges sales tax for NY & NJ residents only: – 13″ 3.1GHz/256GB Space Gray MacBook Pro... Read more
Current deals on 27″ Apple iMacs, models up t...
B&H Photo has 27″ iMacs on sale for up to \$150 off MSRP. Shipping is free, and B&H charges sales tax for NY & NJ residents only: – 27″ 3.8GHz iMac (MNED2LL/A): \$2149 \$150 off MSRP – 27″ 3... Read more
Thursday Deal: 13″ 2.3GHz MacBook Pro for \$11...
B&H Photo has the 13″ 2.3GHz/128GB Space Gray MacBook Pro on sale for \$100 off MSRP. Shipping is free, and B&H charges sales tax for NY & NJ residents only: – 13-inch 2.3GHz/128GB Space... Read more
How to save \$100-\$190 on 10″ & 12″ iPad P...
Apple is now offering Certified Refurbished 2017 10″ and 12″ iPad Pros for \$100-\$190 off MSRP, depending on the model. An Apple one-year warranty is included with each model, and shipping is free: –... Read more
Silver 12″ 1.3GHz MacBook on sale at B&H...
B&H Photo has the 2017 12″ 1.3GHz Silver MacBook on sale for \$1399.99 including free shipping plus sales tax for NY & NJ residents only. Their price is \$200 off MSRP, and it’s the lowest... Read more

## Jobs Board

Firmware Engineer - *Apple* Accessories - A...
# Firmware Engineer - Apple Accessories Job Number: 113452350 Santa Clara Valley, California, United States Posted: 28-Feb-2018 Weekly Hours: 40.00 **Job Summary** Read more
Automation and Performance Engineer, *Apple*...
# Automation and Performance Engineer, Apple Pay Job Number: 113557967 Santa Clara Valley, California, United States Posted: 09-Mar-2018 Weekly Hours: 40.00 **Job Read more
Hardware Systems Architect - *Apple* Watch...
# Hardware Systems Architect - Apple Watch Job Number: 113565323 Santa Clara Valley, California, United States Posted: 05-Mar-2018 Weekly Hours: 40.00 **Job Read more
Lead *Apple* Solution Consultant - Apple (U...
# Lead Apple Solution Consultant Chicago IL Job Number: 113560644 Chicago, Illinois, United States Posted: 10-Mar-2018 Weekly Hours: 40.00 **Job Summary** As a Lead Read more
Art Director, *Apple* Music + Beats1 Market...
# Art Director, Apple Music + Beats1 Marketing Design Job Number: 113258081 Culver City, California, United States Posted: 07-Mar-2018 Weekly Hours: 40.00 **Job Read more