[SciPy-User] Least-square fitting of a 3rd degree polynomial

[Charles R Harris]

On Tue, May 1, 2012 at 11:31 AM, camille chambon <camillechambon at yahoo.fr>wrote:

> Hello,
>
> I would like to fit 'a' such as y = a * x**3 + b * x**2 + c * x + d,
> where x and y are measured data, and b = 0.0, c = -0.00458844157413 and d
> = 5.8 are fixed.
>
> According to
> http://www.scipy.org/Cookbook/FittingData#head-27373a786baa162a2e8a910ee0b8a48838082262,
> I try to use scipy.optimize.leastsq fit routine like that:
>
> #### CODE #####
>
> from numpy import *
> from pylab import *
> from scipy import optimize
>
> # My data points
> x = array([1078.0, 1117.0, 1212.1, 1368.7, 1686.8, 1880.0])
> y = array([5.8, 5.6, 5.4, 4.9, 2.4, 0.0])
>
> b, c, d =  0.0, -0.00458844157413, 5.8 # Fixed parameters
>
> fitfunc = lambda p, x: poly1d([p[0], b, c, d])(x) # Target function
> errfunc = lambda p, x, y: fitfunc(p, x) - y # Distance to the target
> function
>
> p0 = [-4.0E-09] # Initial guess for the parameters
> p1, success = optimize.leastsq(errfunc, p0[:], args=(x, y))
>
> time = linspace(x.min(), x.max(), 100)
> plot(x, y, "ro", time, fitfunc(p1, time), "r-") # Plot of the data and the
> fit
>
> title("a fitting")
> xlabel("time [day]")
> ylabel("number []")
> legend(('x position', 'x fit', 'y position', 'y fit'))
>
> show()
>
> ################################
>
> But the fit does not work (as one can see on the attached image).
>
> Does someone have any idea of what I'm doing wrong?
>
> Thanks in advance for your help.
>
> Cheers,
>
> Camille
>
>
You can use numpy for this.

In [1]: from numpy.polynomial import Polynomial as P

In [2]: x = array([1078.0, 1117.0, 1212.1, 1368.7, 1686.8, 1880.0])

In [3]: y = array([5.8, 5.6, 5.4, 4.9, 2.4, 0.0])

In [4]: p = P.fit(x,y,3)

In [5]: plot(*p.linspace())
Out[5]: [<matplotlib.lines.Line2D at 0x2e51ad0>]

In [6]: plot(x, y, '.')
Out[6]: [<matplotlib.lines.Line2D at 0x2e62c90>]

In [7]: p.mapparms()
Out[7]: (-3.6882793017456361, 0.0024937655860349127)

Note two things, the coefficients go from degree zero upwards and you need
to make the substitution x <- -3.6882793017456361 + 0024937655860349127*x
if you want to use them in a publication.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20120501/01341426/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: example.png
Type: image/png
Size: 19989 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20120501/01341426/attachment.png>