Re: polyfit problem?

[Haisam K. Ido]

Here's what I use instead of polyfit

numobs = rows( y );

i = 1;
j = 0;

while ( i <= numobs )

  while( j <= n ) 

    H(i,j+1) = x(i)**(j);
    j++;

  endwhile

  j = 0;
  i++;

endwhile

# a = the polynomial coefficients

a = inverse( H' * H ) * H' * y;

#--- End Octave ---

Regards,
+-------------------------------------------------------------------+
                 Haisam K. Ido <address@hidden>
           INTELSAT, The Flight Dynamics Section,  Box 60
        3400  International Drive; Washington, DC 20008  USA
         +1 202.944.7654 (voice), +1 202.944.7032 (telefax)
+-------------------------------------------------------------------+           
          

On Fri, 19 Mar 1999 address@hidden wrote:

> > 
> > Hi -
> >  
> > > When I use octave to get a 9th degree polynomial for xy below, via
> > > polyfit I get very strange coefficients, while when I use gnuplot's fit
> > > capability on the same data I get reasonable coefficients.
> > 
> > octave's polyfit does no pre-scaling before setting up the matrix.
> > I find that necessary for any kind of real work with polynomial fits.
> > When I ran into this, I prescaled the data before calling polyfit,
> > rather than "fixing" polyfit.  Sorry.
> > 
> >      - Larry
> > 
> > 
> > 
> 
> Taking Larry's advice, I made some changes to polyfit to do the scaling
> and translating, doing the fit, and then transforming the polynomial 
> back to the original data.  I called it polyfitw (for lack of a better
> name).  It gives the same results as gnuplot for the stated problem, and
> seems to give reasonable results on other test data.  It is attatched
> below. Comments are welcome.
> 
> Bill Lash
> address@hidden
> 
> ## Copyright (C) 1996 John W. Eaton
> ##
> ## This file is part of Octave.
> ##
> ## Octave is free software; you can redistribute it and/or modify it
> ## under the terms of the GNU General Public License as published by
> ## the Free Software Foundation; either version 2, or (at your option)
> ## any later version.
> ##
> ## Octave is distributed in the hope that it will be useful, but
> ## WITHOUT ANY WARRANTY; without even the implied warranty of
> ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
> ## General Public License for more details.
> ##
> ## You should have received a copy of the GNU General Public License
> ## along with Octave; see the file COPYING.  If not, write to the Free
> ## Software Foundation, 59 Temple Place - Suite 330, Boston, MA
> ## 02111-1307, USA.
> 
> ## usage:  [p, yf] = polyfitw (x, y, n)
> ##
> ## Returns the coefficients of a polynomial p(x) of degree n that
> ## minimizes sumsq (p(x(i)) - y(i)), i.e., that best fits the data
> ## in the least squares sense.
> ##
> ## In this version, the x values are scaled and translated so that 
> ## they fit into the range [-1,1] before calculating the polynomial
> ## to avoid bad conditioning of the inverse.  The resulting polynomial
> ## is then transformed back to the original range.
> ##
> ## If two outputs are requested, also return the values of the
> ## polynomial for each value of x.
> 
> ## Author: KH <address@hidden>
> ## Created: 13 December 1994
> ## Adapted-By: jwe
> ## Warping added by: Bill Lash <address@hidden>
> 
> function [p, yf] = polyfitw (x, y, n)
> 
> 
>   if (nargin != 3)
>     usage ("polyfit (x, y, n)");
>   endif
> 
>   if (! (is_vector (x) && is_vector (y) && size (x) == size (y)))
>     error ("polyfit: x and y must be vectors of the same size");
>   endif
> 
>   if (! (is_scalar (n) && n >= 0 && ! isinf (n) && n == round (n)))
>     error ("polyfit: n must be a nonnegative integer");
>   endif
> 
>   y_is_row_vector = (rows (y) == 1);
> 
>   l = length (x);
>   x = reshape (x, l, 1);
>   y = reshape (y, l, 1);
> 
>   ## Map x into the range [-1:1] so that x_warped = x_scale * x + x_xlat
>   x_scale = 2/ (max(x) - min(x));
>   x_xlat  = -(min(x*x_scale) + 1);
>   x_warped = x_scale * x + x_xlat;
> 
>   ## Unfortunately, the economy QR factorization doesn't really save
>   ## memory doing the computation -- the returned values are just
>   ## smaller.
> 
>   ## [Q, R] = qr (X, 0);
>   ## p = flipud (R \ (Q' * y));
> 
>   ## XXX FIXME XXX -- this is probably not so good for extreme values of
>   ## N or X...
> 
>   X = (x * ones (1, n+1)) .^ (ones (l, 1) * (0 : n));
>   X_warped = (x_warped * ones (1, n+1)) .^ (ones (l, 1) * (0 : n));
> 
>   p_warped = (X_warped' * X_warped) \ (X_warped' * y);
> 
>   ## Ok p_warped solved the problem for p_warped[n]*x_warped^n ~= y
>   ## need to determine the polynomial p for p[n]*x^n ~= y.  
> 
>   ## build matrix containing columns of x_warped^n in terms of x
>   if n == 0
>       A = 1;
>   elseif n == 1
>         A=[1,x_xlat;0,x_scale];
>   else
>       A=[1,x_xlat*ones(1,n);0,x_scale*ones(1,n);zeros(n-1,n+1)];
>       for i=2:n+1
>               temp=conv(A(:,i-1),A(:,i));
>               A(:,i)=temp(1:n+1)';
>       endfor
>   endif
> 
>   ## then a simple multiply to get p from p_warped
>   p=A*p_warped;
> 
>   if (nargout == 2)
>     yf = X * p;
> 
>     if (y_is_row_vector)
>       yf = yf';
>     endif
>   endif
> 
>   p = flipud (p);
> 
>   if (! prefer_column_vectors)
>     p = p';
>   endif
> 
> endfunction
> 
>