Re: Problem of polyfit

[Paul Kienzle]

On Oct 6, 2005, at 6:34 AM, Tetsuro KURITA wrote:

> >  p = inv(X'*X)*X'*y
> >    = inv(X)*inv(X')*X'*y
> >    = inv(X)*I*y
> >    = inv(X)*y
> >    = X \ y
> >
> > wpolyfit does the same with weighting on y, except that it uses QR 
> > decomposition to solve X \ y.
>
> In generaly,
>
> inv(a*b) \= inv(b)*inv(a)
>
> If both of a and b are nonsingular matrixes, the following equation is 
> right.
> inv(a*b) = inv(b)*inv(a)
>
> Acutually, in some case, polyfit.m give bad result.

Please check if wpolyfit works better for you.  I tested it using the 
NIST Statistical Reference Datasets 
(http://www.itl.nist.gov/div898/strd/).  See 
octave-forge/main/optim/test_polyfit.m for details.

wpolyfit gives a similar fit to polyfit but does a better job of 
estimating the uncertainty.  In any case, the results are generally 
accurate to 10^-9 or better relative error.

I tried your formula in wsolve (from octave-forge):
        x = inv(A'*A)*(A'*y);
and the results were worse, particularly on the Filippelli test.

If you have a dataset for which polyfit performs particularly poorly 
please post it to the list.

Thanks,

        - Paul



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------