polyfit.m patch

On 13-Jun-2003, Pascal A. Dupuis <[email protected]> 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  <[email protected]>

	* polynomial/polyfit.m: Avoid calling flipud.
	From Pascal A. Dupuis <[email protected]>.
	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
@@ -18,7 +18,7 @@
 ## 02111-1307, USA.
 
 ## -*- texinfo -*-
-## @deftypefn {Function File} {[@var{p}, @var{yf}] =} polyfit (@var{x}, @var{y}, @var{n})
+## @deftypefn {Function File} {[@var{p}, @var{s}] =} polyfit (@var{x}, @var{y}, @var{n})
 ## Return the coefficients of a polynomial @var{p}(@var{x}) of degree
 ## @var{n} that minimizes
 ## @iftex
@@ -33,19 +33,30 @@
 ## @end ifinfo
 ##  to best fit the data in the least squares sense.
 ##
-## The polynomial coefficients are returned in a row vector if @var{x}
-## and @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 @var{x}.
+## If two output arguments are requested, the second is a structure
+## containing the following fields:
+## @table @code
+## @item R
+## The Cholesky factor of the Vandermonde matrix used to compute the
+## polynomial coefficients.
+## @item X
+## The Vandermonde matrix used to compute the polynomial coefficients.
+## @item df
+## The degrees of freedom.
+## @item normr
+## The norm of the residuals.
+## @item yf
+## The values of the polynomial for each value of @var{x}.
+## @end table
 ## @end deftypefn
 
 ## Author: KH <[email protected]>
 ## 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)
@@ -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