From bug-octave-request at bevo dot che dot wisc dot edu  Fri Oct 10 11:27:32 1997
Subject: Error from polyfit -Reply
From: "John W. Eaton" <jwe at bevo dot che dot wisc dot edu>
To: Seung Lee <SJL at nrc dot gov>
cc: bug-octave at bevo dot che dot wisc dot edu
Date: Fri, 10 Oct 1997 11:25:34 -0500

On 10-Oct-1997, Seung Lee <SJL at nrc dot gov> wrote:

| Here is one example:  The obvious soultion is y = 1 + 2*x  + 3*x*x.
| 
|  x          y
| 1.0      6.0
| 2.0    17.0
| 3.0    34.0
| 4.0    57.0
| 5.0    86.0
| 6.0  121.0
| 
| In octave, my command is [p, pf] = polyfit(x, y, 6); and octave
| gives p = [3, 2, 1] and pf = [6, 11, 18, 27, 38, 51] which should be
| same as y.  If I changed statement in polyfit.m as I mentioned
| before, octave gives the same answer as y.  Please check this
| example and let me know your opinion.

I see now that the problem is the the value of yf that is returned.  I
misunderstood and thought you were saying that the polynomial
coefficients were returned incorrectly.  I've fixed this, and also
forced the orientation of the returned yf to match the orientation of
y.  My current version of polyfit is appended below.

Thanks,

jwe


## 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] = polyfit (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.
##
## 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

function [p, yf] = polyfit (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);

  ## 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));

  p = (X' * X) \ (X' * y);

  if (nargout == 2)
    yf = X * p;
  endif

  p = flipud (p);

  if (! prefer_column_vectors)
    p = p';
  endif

  if (y_is_row_vector)
    yf = yf';
  endif

endfunction