From bug-octave-request at bevo dot che dot wisc dot edu  Mon Oct 19 17:31:34 1998
Subject: Unidentified subject!
From: "John W. Eaton" <jwe at bevo dot che dot wisc dot edu>
To: Daniel Tourde <ted at ffa dot se>
Cc: bug-octave at bevo dot che dot wisc dot edu, ted@lorient.ffa.se
Date: Mon, 19 Oct 1998 17:31:35 -0500 (CDT)

On 18-Oct-1998, Daniel Tourde <ted at ffa dot se> wrote:

| The polyfit function of Octave does not work correctly. Results
| given are wrong (especially if the order of the polynome is high (>6))
|
| I have send an exemple two days ago. If you don't find it at
| bug-octave at bevo dot che dot wisc dot edu you will find it at
| help-octave at bevo dot che dot wisc dot edu dot  This example contains two files (y =
| f(x) where x is called globind.gz and y is called displ.gz). I
| joined to this exemple the results given by Matlab (The results of
| Matlab have been validated by an other software.

Please try the following patch.  The reason that the solution was not
handled this way all along is probably because the polyfit function
was written before the left division operator could also solve
non-square problems.

Thanks,

jwe



Mon Oct 19 17:26:35 1998  John W. Eaton  <jwe at bevo dot che dot wisc dot edu>

	* polynomial/polyfit.m: Just use the \ operator to handle the
	least-squares solution.


*** octave-2.0.13/scripts/polynomial/polyfit.m	Tue Nov 18 14:57:24 1997
--- octave-2.0/scripts/polynomial/polyfit.m	Mon Oct 19 17:26:15 1998
***************
*** 51,69 ****
    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;
--- 51,59 ----
    x = reshape (x, l, 1);
    y = reshape (y, l, 1);
  
    X = (x * ones (1, n+1)) .^ (ones (l, 1) * (0 : n));
  
!   p = X \ y;
  
    if (nargout == 2)
      yf = X * p;