From octave-sources-request at bevo dot che dot wisc dot edu  Fri Mar 19 16:30:34 1999
Subject: Re: polyfit problem?
From: lash at tellabs dot com
To: octave-sources at bevo dot che dot wisc dot edu
Cc: h dot ido at intelsat dot int, help-octave@bevo.che.wisc.edu, ldoolitt@jlab.org
Date: Fri, 19 Mar 1999 16:29:53 -0600 (CST)

> 
> 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
lash at tellabs dot com

## 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 <Kurt dot Hornik at ci dot tuwien dot ac dot at>
## Created: 13 December 1994
## Adapted-By: jwe
## Warping added by: Bill Lash <lash at tellabs dot com>

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