help-octave
[Top][All Lists]

## Re: Problem of polyfit

 From: Henry F. Mollet Subject: Re: Problem of polyfit Date: Thu, 06 Oct 2005 12:53:56 -0700 User-agent: Microsoft-Entourage/11.1.0.040913

```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
>    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
-------------------------------------------------------------

```