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

[josef.pktd at gmail.com]

On Wed, May 2, 2012 at 3:35 AM, camille chambon <camillechambon at yahoo.fr> wrote:
> Thanks for your answer.
>
> I tried to reduce my problem to a linear regression problem:
>
> 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
> z = y - b * x**2 - c * x
> p = numpy.polyfit(x, z, 3)
> time = linspace(x.min(), x.max(), 100)
> plot(x, y, "ro", time, numpy.poly1d([p[0], b, c, d])(time), "r-")
>
> but it doesn't work (see attached image).
>
> Where am I wrong?

Are you sure you got your fixed values b,c,d right?

y2 = y - d - c * x - b * x**2
>>> y2
array([ 4.94634002,  4.92528924,  5.16165003,  5.38019998,  4.33978325,
        2.82627016])

It doesn't look like you can fit y2  as a function of x**3 and get a
good result.

I get a similar plot to your original leastsq solution.


import numpy as np
import matplotlib.pyplot as plt


x = np.array([1078.0, 1117.0, 1212.1, 1368.7, 1686.8, 1880.0])
y = np.array([5.8, 5.6, 5.4, 4.9, 2.4, 0.0])
b, c, d =  0.0, -0.00458844157413, 5.8
y2 = y - d - c * x - b * x**2
x2 = x**3
a = np.dot(x2, y2) / np.dot(x2, x2)

#check with statsmodels
import statsmodels.api as sm
a2 = sm.OLS(y2, x2).fit().params
print a, a2[0]

xx = np.linspace(x.min(), x.max())
yhat = d + c * xx + b * xx**2 + a * xx**3
plt.plot(x, y, 'o')
plt.plot(xx, yhat, '-')

plt.figure()
plt.plot(x, y2, 'o')
plt.show()


Josef

>
> Cheers,
>
> Camille
>
> ________________________________
> De : "josef.pktd at gmail.com" <josef.pktd at gmail.com>
> À : camille chambon <camillechambon at yahoo.fr>; SciPy Users List
> <scipy-user at scipy.org>
> Envoyé le : Mardi 1 mai 2012 19h44
> Objet : Re: [SciPy-User] Least-square fitting of a 3rd degree polynomial
>
> 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
>>
>
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>