Friday, September 28, 2012

Least squares curve fitting with a recursive function [+ python code]


I was recently asked by someone to help him use this recursive function -
$$ f(x+1) = f(x) + A(f(x) - B)^3 $$
to fit his experimental data -
x
f(x)
    1    
    235    
2
130
3
115
4
105
5
100
6
99

Objective - 
My objective was to find values of the parameters $A$ and $B$ which would give a good curve fitting. Note that there can be more than one solutions to this. However, the one which gives the smallest norm ie one which fits most tightly should be preferred.
Note that the function is recursive. (I don't know why did he choose such a function.)


Mathematics -
I decided to fit the data using the linear least squares fit technique, since it is the simplest one. You cannot directly use this technique on the above equation since the parameters do not occur linearly.
A little simplification is required -

$ f(x+1) - f(x) = A(f(x) - B)^3 $
$ (f(x+1) - f(x))^{1/3} = A^{1/3}(f(x) - B) $
$ (f(x+1) - f(x))^{1/3} = af(x) + b $, where $a=A^{1/3}$ and $b=-B$

Now, since the parameters occur linearly, I can use the linear least square technique.


Solution -
Using the linear least squares technique, I get -
$a=-0.0244177$, $b=0.93556103$, or
$A=-1.45584135705 \times 10^{-05}$, $B=38.3148771756$

To understand the mathematics behind the least squares curve fitting, you can refer to this article on wikipedia. In numpy, python's numerical computation package, the function numpy.linalg.lstsq can easily finds these coefficients once they occur in a linear fashion. I have used this function in my python code (posted below) to get the values of these coefficients.

Note -
However, another solution set, $A2=-1.1177 \times 10^{-05}$, $B2=24$, also seems to fit the data good enough. I don't know which fitting technique was used to obtain it. I'll update the post when I find it out. I have compared the plots obtained using these 2 solutions below.


Comparative Plots -
The following plot shows all the 3 curves

Plot with $A$, $B$ is my plot. Plot with $A2$, $B2$ is the plot obtained using the other solution .
Note that all the 3 plots have the same values at $x$ = 1, since it is the starting point.

[Click to view a larger image]
Zoomed view of the region where plots do not match => $x=2, 3, 4, 5, 6$
Note that my parameters fit the curve better at some points, while the other ones fits them better at the remaining points.
[Click to view a larger image]



Observation -
If you observe the $2^{nd}$ plot carefully , my parameters fits the curve well for $2/5$ points, while the other parameters fit the data better for $3/5$ points. Which one is better? Well it depends on the criteria. (I'll soon elaborate this point when I would get the exact way how the other solution - $A2$, $B2$ - was obtained.)


Python Program -
The following python program was used in solving for $A$ and $B$ and plotting the 3 graphs
(Does not YET include the method for obtaining $A2$, $B2$) -