Re: Problem of polyfit

[Henry F. Mollet]

Here's a test using N = 1.  Polyfit.m gives the expected results as per
least square fit using the hat matrix H.
Henry
octave:29> x'
ans =

   1   2   3   4   5   6   7   8   9  10

octave:30> y' % y = 2*x and changed 5th element to 15
ans =

   2   4   6   8  15  12  14  16  18  20

octave:31> Onex'
ans =

   1   1   1   1   1   1   1   1   1   1
   1   2   3   4   5   6   7   8   9  10

octave:32> yf' % using [P, S] = polyfit (x, y, 1) and saving yf
ans =

 Columns 1 through 8:

   2.6364   4.6061   6.5758   8.5455  10.5152  12.4848  14.4545  16.4242

 Columns 9 and 10:

  18.3939  20.3636

octave:33> H = Onex*(Onex'*Onex)^-1*Onex';
octave:34> yhat = H*y;
octave:35> yhat'
ans =

 Columns 1 through 8:

   2.6364   4.6061   6.5758   8.5455  10.5152  12.4848  14.4545  16.4242

 Columns 9 and 10:

  18.3939  20.3636

octave:36> (yhat-yf)'
ans =

 Columns 1 through 5:

   -3.6364e-05   -3.9394e-05   -4.2424e-05   -4.5455e-05   -4.8485e-05

 Columns 6 through 10:

    4.8485e-05    4.5455e-05    4.2424e-05    3.9394e-05    3.6364e-05




on 10/6/05 2:46 AM, Tetsuro KURITA at address@hidden wrote:

> The following code in "polyfit.m" never give the best fit data in the
> least squares sense as described in document.
> 
>    X = (x * ones (1, n+1)) .^ (ones (l, 1) * (n : -1 : 0));
>    p = X \ y;
> 
> This code do nothing to mimimize `sumsq (p(x(i)) - y(i))'.
> 
> I think this code should be modified as follows.
> 
> X = (x * ones (1, n+1)) .^ (ones (l, 1) * (n: -1 : 0))
> W = X'*X
> z = X'*y
> p = inv(W)*z
> 
> I am using Octave 2.1.57 on MacOS X.
> 
> =======================================================
>   Tetsuro KURITA
>    E-mail: address@hidden
>    http://homepage.mac.com/tkurita/scriptfactory/
> =======================================================
> 
> 
> 
> -------------------------------------------------------------
> Octave is freely available under the terms of the GNU GPL.
> 
> Octave's home on the web:  http://www.octave.org
> How to fund new projects:  http://www.octave.org/funding.html
> Subscription information:  http://www.octave.org/archive.html
> -------------------------------------------------------------
> 




-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------