From octave-sources-request at bevo dot che dot wisc dot edu Tue Mar 23 09:47:59 1999 Subject: Re: polyfit problem? From: "Haisam K. Ido" To: octave-sources at bevo dot che dot wisc dot edu, help-octave@bevo.che.wisc.edu Date: Tue, 23 Mar 1999 15:47:59 +0000 (GMT) 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 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 lash at tellabs dot com 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 > 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 > ## Created: 13 December 1994 > ## Adapted-By: jwe > ## Warping added by: Bill Lash > > 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 > >