From octave-sources-request at bevo dot che dot wisc dot edu Thu Aug 28 16:01:10 2003
Subject: polyfit.m patch
From: "John W. Eaton" <jwe at bevo dot che dot wisc dot edu>
To: "Pascal A. Dupuis" <Pascal dot Dupuis at esat dot kuleuven dot ac dot be>
Cc: octave-sources at bevo dot che dot wisc dot edu
Date: Thu, 28 Aug 2003 16:00:40 -0500

On 13-Jun-2003, Pascal A. Dupuis <Pascal dot Dupuis at esat dot kuleuven dot ac dot be> wrote:

| please find herewith enclosed a patch for m/polynomial/polyfit.m.
| The goal is to also return the design matrix, that is, the matrix such
| that yf = X*p. This permits to estimate the standard deviation on
| the coefficients of p: let sy be the square of the standard deviation
| between y and yf, then sigma_p = diag(sqrt(sy * inv(X.'*X)))

Instead of doing exactly what you proposed, I've made the following
changes so Octave will be more compatible with Matlab.

This change is not backward-compatible with previous versions of
Octave because now the second output value is a structure.  But the
structure does contain the information that used to be returned in the
second output value, so maybe the change will not be too painful for
people who have used that feature.

jwe


scripts/ChangeLog:

2003-08-28  John W. Eaton  <jwe at bevo dot che dot wisc dot edu>

	* polynomial/polyfit.m: Avoid calling flipud.
	From Pascal A. Dupuis <Pascal dot Dupuis at esat dot kuleuven dot ac dot be>.
	Return structure as second output value for improved Matlab
	compatibility.


Index: scripts/polynomial/polyfit.m
===================================================================
RCS file: /usr/local/cvsroot/octave/scripts/polynomial/polyfit.m,v
retrieving revision 1.19
diff -u -r1.19 polyfit.m
--- scripts/polynomial/polyfit.m	9 Aug 2002 18:58:15 -0000	1.19
+++ scripts/polynomial/polyfit.m	28 Aug 2003 20:58:38 -0000
 at @ -18,7 +18,7 @@
 ## 02111-1307, USA.
 
 ## -*- texinfo -*-
-##  at deftypefn {Function File} {[@var{p}, @var{yf}] =} polyfit (@var{x}, @var{y}, @var{n})
+##  at deftypefn {Function File} {[@var{p}, @var{s}] =} polyfit (@var{x}, @var{y}, @var{n})
 ## Return the coefficients of a polynomial  at var{p}(@var{x}) of degree
 ##  at var{n} that minimizes
 ##  at iftex
 at @ -33,19 +33,30 @@
 ##  at end ifinfo
 ##  to best fit the data in the least squares sense.
 ##
-## The polynomial coefficients are returned in a row vector if  at var{x}
-## and  at var{y} are both row vectors; otherwise, they are returned in a
-## column vector.
+## The polynomial coefficients are returned in a row vector.
 ##
-## If two output arguments are requested, the second contains the values of
-## the polynomial for each value of  at var{x} dot 
+## If two output arguments are requested, the second is a structure
+## containing the following fields:
+##  at table @code
+##  at item R
+## The Cholesky factor of the Vandermonde matrix used to compute the
+## polynomial coefficients.
+##  at item X
+## The Vandermonde matrix used to compute the polynomial coefficients.
+##  at item df
+## The degrees of freedom.
+##  at item normr
+## The norm of the residuals.
+##  at item yf
+## The values of the polynomial for each value of  at var{x} dot 
+##  at end table
 ##  at end deftypefn
 
 ## 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)
+function [p, s, mu] = polyfit (x, y, n)
 
 
   if (nargin != 3)
 at @ -66,22 +77,25 @@
   x = reshape (x, l, 1);
   y = reshape (y, l, 1);
 
-  X = (x * ones (1, n+1)) .^ (ones (l, 1) * (0 : n));
+  X = (x * ones (1, n+1)) .^ (ones (l, 1) * (n : -1 : 0));
 
   p = X \ y;
 
-  if (nargout == 2)
-    yf = X * p;
+  if (nargout > 1)
+
+    yf = X*p;
 
     if (y_is_row_vector)
-      yf = yf.';
+      s.yf = yf.';
+    else
+      s.yf = yf;
     endif
-  endif
 
-  p = flipud (p);
+    [s.R, dummy] = chol (X'*X);
+    s.X = X;
+    s.df = l - n - 1;
+    s.normr = norm (yf - y);
 
-  if (y_is_row_vector && rows (x) == 1)
-    p = p';
   endif
 
 endfunction