Re: [Numpy-discussion] Polyfit

[Charles R H. <cha...@gm...>]

On 10/12/06, Charles R Harris <cha...@gm...> wrote:
>
>
>
> On 10/12/06, Greg Willden <gre...@gm...> wrote:
> >
> > On 10/12/06, Charles R Harris <cha...@gm... > wrote:
> > >
> > > I'm guessing that the rcond number in the lstsq version (default
> > > 1e-10) is the difference. Generally the lstsq version should work better
> > > than the MPL version because at*a is not as well conditioned and vandermonde
> > > matrices are notoriously ill conditioned anyway for higher degree fits. It
> > > would help if you attached the data files in ascii form unless they happen
> > > to contain thousands of data items. Just the x will suffice, zip it up if
> > > you have to.
> > >
> >
> >
> > Here are the files.
> >
> > Since the two algorithms behave differently and each has it place then
> > can both be included in numpy?
> > i.e. numpy.polyfit(x,y,N, mode='alg1')
> > numpy.polyfit (x,y,N, mode='alg2')
> >
> > replacing alg1 and alg2 with meaningful names.
> >
>
> The polyfit function looks seriously busted. If I do the fits by hand I
> get the same results using the (not so hot) MPL version or lstsq. I don't
> know what the problem is. The docstring is also incorrect for the method.
> Hmmm...
>

Polyfit seems overly conservative in its choice of rcond.

In [101]: lin.lstsq(v,y,1e-10)[0]
Out[101]:
array([  5.84304475e-07,  -5.51513630e-03,   1.51465472e+01,
         3.05631361e-02])
In [107]: polyfit(x,y,3)
Out[108]:
array([  5.84304475e-07,  -5.51513630e-03,   1.51465472e+01,
         3.05631361e-02])

Compare

In [112]: lin.lstsq(v,y,1e-12)[0]
Out[112]:
array([ -5.42970700e-07,   1.61425067e-03,   1.99260667e+00,
         6.51889107e+03])

In [113]: dot(lin.inv(vtv),dot(v.T,y))
Out[113]:
array([ -5.42970700e-07,   1.61425067e-03,   1.99260667e+00,
         6.51889107e+03])

So the default needs to be changed somewhere. Probably polyfit shoud accept
rcond as a keyword. Where the actual problem lies is a bit obscure as the
normal rcond default for lin.lstsq is 1e-12. Maybe some sort of import error
somewhere down the line.

Chuck