Problem of polyfit

[Tetsuro KURITA]

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: ***@mac.com
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
-------------------------------------------------------------

[Paul Kienzle]

Post by Tetsuro KURITA
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

p = inv(X'*X)*X'*y
= inv(X)*inv(X')*X'*y
= inv(X)*I*y
= inv(X)*y
= X \ y

wpolyfit does the same with weighting on y, except that it uses QR
decomposition to solve X \ y.


- Paul



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

[Tetsuro KURITA]

Post by Paul Kienzle
p = inv(X'*X)*X'*y
= inv(X)*inv(X')*X'*y
= inv(X)*I*y
= inv(X)*y
= X \ y
wpolyfit does the same with weighting on y, except that it uses QR
decomposition to solve X \ y.

In generaly,

inv(a*b) \= inv(b)*inv(a)

If both of a and b are nonsingular matrixes, the following equation is
right.
inv(a*b) = inv(b)*inv(a)

Acutually, in some case, polyfit.m give bad result.

=======================================================
Tetsuro KURITA
E-mail: ***@mac.com
=======================================================



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

[Paul Kienzle]

Post by Tetsuro KURITA

Post by Paul Kienzle
p = inv(X'*X)*X'*y
= inv(X)*inv(X')*X'*y
= inv(X)*I*y
= inv(X)*y
= X \ y
wpolyfit does the same with weighting on y, except that it uses QR
decomposition to solve X \ y.

In generaly,
inv(a*b) \= inv(b)*inv(a)
If both of a and b are nonsingular matrixes, the following equation is
right.
inv(a*b) = inv(b)*inv(a)
Acutually, in some case, polyfit.m give bad result.

Please check if wpolyfit works better for you. I tested it using the
NIST Statistical Reference Datasets
(http://www.itl.nist.gov/div898/strd/). See
octave-forge/main/optim/test_polyfit.m for details.

wpolyfit gives a similar fit to polyfit but does a better job of
estimating the uncertainty. In any case, the results are generally
accurate to 10^-9 or better relative error.

I tried your formula in wsolve (from octave-forge):
x = inv(A'*A)*(A'*y);
and the results were worse, particularly on the Filippelli test.

If you have a dataset for which polyfit performs particularly poorly
please post it to the list.

Thanks,

- Paul



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

[Tetsuro KURITA]

Post by Paul Kienzle
wpolyfit gives a similar fit to polyfit but does a better job of
estimating the uncertainty. In any case, the results are generally
accurate to 10^-9 or better relative error.

wpolyfit show following warnning, and result is same to the on of
polyfit.m .

warning: in /sw/share/octave/2.1.57/site/m/octave-forge/optim/wsolve.m
near line 76, column 5:
warning: matrix singular to machine precision, rcond = 1.73942e-56
warning: attempting to find minimum norm solution

Post by Paul Kienzle
If you have a dataset for which polyfit performs particularly poorly
please post it to the list.

Following is an example which show that polyfit.m does not works well.

y = [2.720;
2.811;
2.899;
2.946;
2.940;
2.956;
2.914;
2.949;
2.897;
2.821;
2.736;
2.657;
2.565;
2.363;
2.251;
2.109;
2.024;
1.981;
1.996;
2.068;
2.132;
2.263;
2.408;
2.544;
2.689;
2.825;
2.901;
2.939;
2.932;
2.966;
3.027;
3.129]

x = [900058.30 ;
1000061.71;
1100054.15;
1200085.73;
1300026.64;
1400014.73;
1500007.74;
1600050.93;
1700013.05;
1800024.58;
1900039.38;
2000011.13;
2100062.61;
2200059.35;
2300035.86;
2400023.69;
2500051.17;
2600062.55;
2700001.99;
2800050.76;
2900066.95;
3000066.57;
3100061.96;
3200061.67;
3300070.35;
3400015.58;
3500041.39;
3600066.28;
3700009.04;
3800063.72;
3900001.42;
4000001.16]

p1 = polyfit(x,y,8);
y1 = polyval(p1,x);

gplot ([x,y]),([x,y1])

p2 = wpolyfit(x,y,8);
y2 = polyval(p2,x);
gplot ([x,y]),([x,y2])

=======================================================
Tetsuro KURITA
E-mail: ***@mac.com
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
-------------------------------------------------------------

[Dmitri A. Sergatskov]

Post by Tetsuro KURITA
Following is an example which show that polyfit.m does not works well.

You would get a very good fit if you normalize your data

x = (x-mean(x))/(max(x)-min(x))
and the same with y.

Regards,

Dmitri

[Tetsuro KURITA]

Post by Dmitri A. Sergatskov
You would get a very good fit if you normalize your data

Thanks.
It seems that my opinion was wrong.
I feel embarrassed about making much ado.

Why normailing data is requred.

But there was a person who meets same trouble to mine.

http://www.octave.org/octave-lists/archive/help-octave.1999/
msg01007.html
http://www.octave.org/octave-lists/archive/help-octave.1999/
msg01008.html

I hope polyfit will be improved for less hustle.

Also following two equations do not return same results, when
calculating data I posted.

p = inv(X'*X)*(X'*y);

p = (X'*X)\(X'*y);

In the case of data I posted, first equation give better results. It
looks that the operator "\" have a certain limitation. What is that ?

=======================================================
Tetsuro KURITA
E-mail: ***@mac.com
=======================================================



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

[Fredrik Lingvall]

Post by Tetsuro KURITA
Also following two equations do not return same results, when
calculating data I posted.
p = inv(X'*X)*(X'*y);
p = (X'*X)\(X'*y);
In the case of data I posted, first equation give better results. It
looks that the operator "\" have a certain limitation. What is that ?

Can you post the output of norm( inv(X'*X)*X'*y - (X'*X)\X'*y )?

A simple test on my machine:

octave:26> X=randn(100,100);
octave:27> y=randn(100,1);
octave:28> norm( inv(X'*X)*X'*y - (X'*X)\X'*y )
ans = 3.7190e-13

Fredrik



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

[Tetsuro KURITA]

Post by Fredrik Lingvall
Can you post the output of norm( inv(X'*X)*X'*y - (X'*X)\X'*y )?

The result is as follows.

octave> norm(inv(X'*X)*(X'*y) - (X'*X)\(X'*y))
warning: matrix singular to machine precision, rcond = 3.02246e-112
warning: attempting to find minimum norm solution
ans = 50.504

When calculating (X'*X)\(X'*y), warning message raise.


=======================================================
Tetsuro KURITA
E-mail: ***@mac.com
=======================================================



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

[Fredrik Lingvall]

Post by Tetsuro KURITA

Post by Fredrik Lingvall
Can you post the output of norm( inv(X'*X)*X'*y - (X'*X)\X'*y )?

The result is as follows.
octave> norm(inv(X'*X)*(X'*y) - (X'*X)\(X'*y))
warning: matrix singular to machine precision, rcond = 3.02246e-112
warning: attempting to find minimum norm solution
ans = 50.504
When calculating (X'*X)\(X'*y), warning message raise.
=======================================================
Tetsuro KURITA
=======================================================

OK. Then X'*X do not have full rank (i.e., som eigenvalues are very
small) and the inverse do not exist.

Try for example,

[U,D,V] = svd(X'*X);
semilogy(diag(D));

You will probably see that some singular values are very close to zero
and a direct inverse of
X'*X will result in strong noise amplification.

Fredrik





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

[Paul Kienzle]

Post by Tetsuro KURITA

Post by Fredrik Lingvall
Can you post the output of norm( inv(X'*X)*X'*y - (X'*X)\X'*y )?

The result is as follows.
octave> norm(inv(X'*X)*(X'*y) - (X'*X)\(X'*y))
warning: matrix singular to machine precision, rcond = 3.02246e-112
warning: attempting to find minimum norm solution
ans = 50.504
When calculating (X'*X)\(X'*y), warning message raise.

I get warnings for both:

warning: inverse: matrix singular to machine precision, rcond =
3.02246e-112
warning: matrix singular to machine precision, rcond = 3.02246e-112

The norm is the same.

OS X 10.3 octave 2.1.71 HPC build.

- Paul



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

[Paul Kienzle]

Post by Dmitri A. Sergatskov

Post by Tetsuro KURITA
Following is an example which show that polyfit.m does not works well.

You would get a very good fit if you normalize your data
x = (x-mean(x))/(max(x)-min(x))
and the same with y.

I've modified wpolyfit so that

[p,s,mu] = wpolyfit(x,y,dy,n)

first centers and scales x by mu = [mean(x),std(x)]. I've used std(x)
rather than max(x)-min(x) for compatibility.

The fit is a whole lot nicer when you do this :-)

- Paul



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

[Fredrik Lingvall]

Post by Paul Kienzle

Post by Tetsuro KURITA
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

p = inv(X'*X)*X'*y
= inv(X)*inv(X')*X'*y
= inv(X)*I*y
= inv(X)*y
= X \ y
wpolyfit does the same with weighting on y, except that it uses QR
decomposition to solve X \ y.

Try,

[U,D,V] = svd(X);
semilogy(diag(D));

If you have singular values close to zero your in trouble.

Fredrik



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

[Dmitri A. Sergatskov]

Post by Paul Kienzle
wpolyfit does the same with weighting on y, except that it uses QR
decomposition to solve X \ y.

BTW, I just noticed that polifit returns p as 1xn vector and wpolifit
as nx1.
I do not know which one is "right", but I think it should be consistent.

Post by Paul Kienzle
- Paul

Sincerely,

Dmitri.
--



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

[John W. Eaton]

On 6-Oct-2005, Dmitri A. Sergatskov wrote:

| On 10/6/05, Paul Kienzle <***@users.sourceforge.net> wrote:
|
| > wpolyfit does the same with weighting on y, except that it uses QR
| > decomposition to solve X \ y.
| >
|
| BTW, I just noticed that polifit returns p as 1xn vector and wpolifit
| as nx1.
| I do not know which one is "right", but I think it should be consistent.

For compatibility, I think polyfit should always return a row vector.

jwe



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

[Paul Kienzle]

Post by John W. Eaton
|
| > wpolyfit does the same with weighting on y, except that it uses QR
| > decomposition to solve X \ y.
| >
|
| BTW, I just noticed that polifit returns p as 1xn vector and wpolifit
| as nx1.
| I do not know which one is "right", but I think it should be
consistent.
For compatibility, I think polyfit should always return a row vector.

I've made the change to wpolyfit.

- Paul



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

[Tom Holroyd]

Post by Tetsuro KURITA
W = X'*X
z = X'*y
p = inv(W)*z

Is that the same as pinv?

Dr. Tom Holroyd
"A man of genius makes no mistakes. His errors are volitional and
are the portals of discovery." -- James Joyce



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

[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

Post by Tetsuro KURITA
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
-------------------------------------------------------------

[Miroslaw Kwasniak]

Post by Henry F. Mollet
octave:36> (yhat-yf)'
ans =
-3.6364e-05 -3.9394e-05 -4.2424e-05 -4.5455e-05 -4.8485e-05
4.8485e-05 4.5455e-05 4.2424e-05 3.9394e-05 3.6364e-05

Hmm, you have single precision aritmetic?

OCTAVE_VERSION = 2.1.50
-3.1086e-15 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00
0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00

OCTAVE_VERSION = 2.9.2
-8.8818e-16 -4.4409e-15 -7.1054e-15 -3.5527e-15 -1.7764e-15
-3.5527e-15 -5.3291e-15 -3.5527e-15 0.0000e+00 3.5527e-15



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

[Henry F. Mollet]

[~] -bash-2.05b 501$ octave
GNU Octave, version 2.1.71 (powerpc-apple-darwin8.1.0).
This time around I was starting out with format long instead of default
format. I did not know that this was required to get higher precision. I
thought that my calculations using default format would be equally precise
but just not show as many digits on the screen.
Thanks, Henry
octave:1> format long
octave:22> (yhat-yf)'
octave:22> (yhat-yf)'
ans =

Columns 1 through 3:

3.10862446895044e-15 5.32907051820075e-15 7.10542735760100e-15

Columns 4 through 6:

7.10542735760100e-15 5.32907051820075e-15 5.32907051820075e-15

Columns 7 through 9:

5.32907051820075e-15 7.10542735760100e-15 3.55271367880050e-15

Column 10:

7.10542735760100e-15

Post by Miroslaw Kwasniak

Post by Henry F. Mollet
octave:36> (yhat-yf)'
ans =
-3.6364e-05 -3.9394e-05 -4.2424e-05 -4.5455e-05 -4.8485e-05
4.8485e-05 4.5455e-05 4.2424e-05 3.9394e-05 3.6364e-05

Hmm, you have single precision aritmetic?
OCTAVE_VERSION = 2.1.50
-3.1086e-15 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00
0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00
OCTAVE_VERSION = 2.9.2
-8.8818e-16 -4.4409e-15 -7.1054e-15 -3.5527e-15 -1.7764e-15
-3.5527e-15 -5.3291e-15 -3.5527e-15 0.0000e+00 3.5527e-15
-------------------------------------------------------------
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
-------------------------------------------------------------

[Miroslaw Kwasniak]

Post by Henry F. Mollet
[~] -bash-2.05b 501$ octave
GNU Octave, version 2.1.71 (powerpc-apple-darwin8.1.0).
This time around I was starting out with format long instead of default
format. I did not know that this was required to get higher precision.

It sounds horrible , none of my versions of octave has such behaviour !



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

[John W. Eaton]

On 7-Oct-2005, Miroslaw Kwasniak wrote:

| On Thu, Oct 06, 2005 at 06:44:20PM -0700, Henry F. Mollet wrote:
| > [~] -bash-2.05b 501$ octave
| > GNU Octave, version 2.1.71 (powerpc-apple-darwin8.1.0).
| > This time around I was starting out with format long instead of default
| > format. I did not know that this was required to get higher precision.
|
| It sounds horrible , none of my versions of octave has such behaviour !

None of mine do either. The format command only controls how the
output is presented. It has nothing to do with the way calculations
are performed. So if you see a difference, there must be some
explanation other than format.

jwe



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

[Henry F. Mollet]

I believe my problem was/is as follows:
[P, S] = polyfit (x, y, 1)
does not give me yf in the workspace unless I amend polyfit.m.
Therefore I copied the result from the screen and I used yf = [paste]. I
lost precision in this step. The rest was done line by line and yhat was in
the workspace with high precision despite using default format.

So my question is how can I get yf (of polyfit) into my workspace without
changing polyfit.m and without going to format long.
Henry

Post by Miroslaw Kwasniak

Post by Henry F. Mollet
[~] -bash-2.05b 501$ octave
GNU Octave, version 2.1.71 (powerpc-apple-darwin8.1.0).
This time around I was starting out with format long instead of default
format. I did not know that this was required to get higher precision.

It sounds horrible , none of my versions of octave has such behaviour !
-------------------------------------------------------------
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
-------------------------------------------------------------

[Miroslaw Kwasniak]

Post by Henry F. Mollet
[P, S] = polyfit (x, y, 1)
does not give me yf in the workspace unless I amend polyfit.m.
Therefore I copied the result from the screen and I used yf = [paste]. I
lost precision in this step. The rest was done line by line and yhat was in
the workspace with high precision despite using default format.
So my question is how can I get yf (of polyfit) into my workspace without
changing polyfit.m and without going to format long.
Henry

See help for polyfit ;)

2.1.50: Function File: [P, YF] = polyfit (X, Y, N)
2.9.2: Function File: [P, S] = polyfit (X, Y, N)
^^^

So:

octave:2> x=1:2;[p yf]=polyfit (x,x,1); yf
yf = ^^^^^^^^

1 2

octave2.9:1> x=1:2;[p s]=polyfit (x,x,1); s.yf
ans = ^^^^^^^^^

1 2



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

[Henry F. Mollet]

Output S of polyfit is a structure with 5 elements. 'S.fy' is what I should
have used instead of copy/paste from the output on the screen. I have set
silent_functions = 0, therefore both P and S were in the workspace.
Henry

Post by Henry F. Mollet
[P, S] = polyfit (x, y, 1)
does not give me yf in the workspace unless I amend polyfit.m.
Therefore I copied the result from the screen and I used yf = [paste]. I
lost precision in this step. The rest was done line by line and yhat was in
the workspace with high precision despite using default format.
So my question is how can I get yf (of polyfit) into my workspace without
changing polyfit.m and without going to format long.
Henry

Post by Miroslaw Kwasniak

Post by Henry F. Mollet
[~] -bash-2.05b 501$ octave
GNU Octave, version 2.1.71 (powerpc-apple-darwin8.1.0).
This time around I was starting out with format long instead of default
format. I did not know that this was required to get higher precision.

It sounds horrible , none of my versions of octave has such behaviour !
-------------------------------------------------------------
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
-------------------------------------------------------------

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