Least-square fitting of a 3rd degree polynomial

classic Classic list List threaded Threaded
9 messages Options
Reply | Threaded
Open this post in threaded view
|

Least-square fitting of a 3rd degree polynomial

CamilleC
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


_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user

a_fit.png (31K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Least-square fitting of a 3rd degree polynomial

josef.pktd
On Tue, May 1, 2012 at 1:31 PM, camille chambon <[hidden email]> 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
> [hidden email]
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: Least-square fitting of a 3rd degree polynomial

Charles R Harris
In reply to this post by CamilleC


On Tue, May 1, 2012 at 11:31 AM, camille chambon <[hidden email]> 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


_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user

example.png (26K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Least-square fitting of a 3rd degree polynomial

Charles R Harris
In reply to this post by josef.pktd


On Tue, May 1, 2012 at 11:44 AM, <[hidden email]> wrote:
On Tue, May 1, 2012 at 1:31 PM, camille chambon <[hidden email]> 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


It does recognize NA in devel, but we will see where that goes ;)

Chuck


_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re : Least-square fitting of a 3rd degree polynomial

CamilleC
In reply to this post by josef.pktd
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?

Cheers,

Camille


De : "[hidden email]" <[hidden email]>
À : camille chambon <[hidden email]>; SciPy Users List <[hidden email]>
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 <[hidden email]> 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
> [hidden email]
> http://mail.scipy.org/mailman/listinfo/scipy-user
>



_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user

a_fit.png (37K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re : Least-square fitting of a 3rd degree polynomial

CamilleC
In reply to this post by Charles R Harris
Thanks for your answer.

But I would like to fit 'a' while 'b', 'c' and 'd' are fixed. Your method gives me another b, c and d. Is there any way to fit a polynomial fixing some coefficients?

Cheers,

Camille


De : Charles R Harris <[hidden email]>
À : camille chambon <[hidden email]>; SciPy Users List <[hidden email]>
Envoyé le : Mardi 1 mai 2012 19h45
Objet : Re: [SciPy-User] Least-square fitting of a 3rd degree polynomial



On Tue, May 1, 2012 at 11:31 AM, camille chambon <[hidden email]> 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




_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re : Least-square fitting of a 3rd degree polynomial

CamilleC
In reply to this post by Charles R Harris
I tried to set a polynomial with a NA coefficient:

from numpy.polynomial import Polynomial as P
import numpy
x = numpy.array([1078.0, 1117.0, 1212.1, 1368.7, 1686.8, 1880.0])
y = numpy.array([5.8, 5.6, 5.4, 4.9, 2.4, 0.0])
b, c, d =  0.0, -0.00458844157413, 5.8
my_poly = P([d, c, b, numpy.nan], [x.min(), x.max()])
print my_poly # print "poly([  5.80000000e+00  -4.58844157e-03   0.00000000e+00              nan], [ 1078.  1880.])"
p = my_poly.fit(x,y,3)
print p # print "poly([ 4.20020036 -2.66837734 -1.33882427 -0.20317739], [ 1078.  1880.])"

but new coefficients are calculated. Does that mean it doesn't work?

Cheers,

Camille


De : Charles R Harris <[hidden email]>
À : SciPy Users List <[hidden email]>
Envoyé le : Mardi 1 mai 2012 19h47
Objet : Re: [SciPy-User] Least-square fitting of a 3rd degree polynomial



On Tue, May 1, 2012 at 11:44 AM, <[hidden email]> wrote:
On Tue, May 1, 2012 at 1:31 PM, camille chambon <[hidden email]> 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


It does recognize NA in devel, but we will see where that goes ;)

Chuck


_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user



_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: Re : Least-square fitting of a 3rd degree polynomial

josef.pktd
In reply to this post by CamilleC
On Wed, May 2, 2012 at 3:35 AM, camille chambon <[hidden email]> 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 : "[hidden email]" <[hidden email]>
> À : camille chambon <[hidden email]>; SciPy Users List
> <[hidden email]>
> 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 <[hidden email]>
> 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
>> [hidden email]
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>
>
>
> _______________________________________________
> SciPy-User mailing list
> [hidden email]
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re : Re : Least-square fitting of a 3rd degree polynomial

CamilleC
Thanks for your answer.

I realized I made a mistake with my data. It works much better now. Sorry, I have my head in the clouds...

Cheers,

Camille


De : "[hidden email]" <[hidden email]>
À : SciPy Users List <[hidden email]>
Envoyé le : Mercredi 2 mai 2012 12h45
Objet : Re: [SciPy-User] Re : Least-square fitting of a 3rd degree polynomial

On Wed, May 2, 2012 at 3:35 AM, camille chambon <[hidden email]> 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 : "[hidden email]" <[hidden email]>
> À : camille chambon <[hidden email]>; SciPy Users List
> <[hidden email]>
> 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 <[hidden email]>
> 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
>> [hidden email]
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>
>
>
> _______________________________________________
> SciPy-User mailing list
> [hidden email]
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user



_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user