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

josef.pktd at gmail.com josef.pktd at gmail.com
Tue May 1 13:44:01 EDT 2012


On Tue, May 1, 2012 at 1:31 PM, 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.

my guess would be that the scale is awful, the x are very large

if you only need the to estimate "a", then it's just a very simple
linear regression problem dot(x,x)/dot(x,y) or something like this.

as aside numpy.polynomial works very well for fitting a polynomial,
but doesn't allow for missing terms

Josef

>
> Cheers,
>
> Camille
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>



More information about the SciPy-User mailing list