forcing y intercept with polyfit

[Seth Gilchrist]

I want to use polyfit in a script file and force the y-intercept of
the fit line to zero. Is there a way to do this besides setting the
(n+1) value generated by p = polyfit(x,y,n) to zero manually?

Thanks for the help
Seth

[Steven Lord]

Post by Seth Gilchrist
I want to use polyfit in a script file and force the y-intercept of
the fit line to zero. Is there a way to do this besides setting the
(n+1) value generated by p = polyfit(x,y,n) to zero manually?
Thanks for the help
Seth

Well, if the y-intercept of the line is zero (i.e. the line goes through the
origin) then your equation must be:

y = a_n*x^n+a_(n-1)*x^(n-1)...+a_1*x

Now if you divide through by x, you get:

y/x = a_n*x^(n-1)+a_(n-1)*x^(n-2)...+a_2*x+a_1

which is an (n-1)st degree polynomial relating x to y/x. [In the case where
x is zero, you know y should be zero so replace any elements of (y/x) where
x==0 by 0.] Now you can polyfit x and (y/x) with an n-1 degree polynomial
and shift all the coefficients left by 1 place.

[Tom Lane]

Post by Steven Lord

Post by Seth Gilchrist
I want to use polyfit in a script file and force the y-intercept of
the fit line to zero.

...

Post by Steven Lord
y/x = a_n*x^(n-1)+a_(n-1)*x^(n-2)...+a_2*x+a_1

Steve's idea is a good one if you are willing to minimize the sum of squared
residuals on this scale:

sum((y/x - estimated y/x)^2)

(using a notation that is somewhere in between math and MATLAB). If you
want to minimize the sum of squared residuals on the original scale:

sum((y-estimated y)^2),

then you can do something like this to fit a polynomial up to 3rd degree:

b = [x.^3, x.^2, x] \ y;

This gives you the coefficient estimates. If you need the 2nd and 3rd
outputs from polyfit you'd have to do a bit more. You shouldn't use polyfit
the usual way and then set the last element of p to 0 manually -- that would
give the wrong answer. You should, though, add an extra 0 element for the
constant term to b, if you intend to use polyval to evaluate the polynomial
later on.

Here's an example:

% Random data
x = rand(50,1);
y = .2 + x - x.^2 + randn(size(x))/10;
scatter(x,y)
xx = linspace(0,1);

% Unconstrained fit
b1 = polyfit(x,y,3);
line(xx,polyval(b1,xx),'color','r');

% Constrained to pass through the origin
b2 = [x.^3 x.^2 x]\y;
line(xx,polyval([b2; 0],xx),'color','g');


I hope this is clear enough. Let me know if not.

-- Tom

[Scott Seidman]

Post by Tom Lane

Post by Steven Lord

Post by Seth Gilchrist
I want to use polyfit in a script file and force the y-intercept of
the fit line to zero.

...

Post by Steven Lord
y/x = a_n*x^(n-1)+a_(n-1)*x^(n-2)...+a_2*x+a_1

Steve's idea is a good one if you are willing to minimize the sum of
sum((y/x - estimated y/x)^2)
(using a notation that is somewhere in between math and MATLAB). If
you want to minimize the sum of squared residuals on the original
sum((y-estimated y)^2),
b = [x.^3, x.^2, x] \ y;
This gives you the coefficient estimates. If you need the 2nd and 3rd
outputs from polyfit you'd have to do a bit more. You shouldn't use
polyfit the usual way and then set the last element of p to 0 manually
-- that would give the wrong answer. You should, though, add an extra
0 element for the constant term to b, if you intend to use polyval to
evaluate the polynomial later on.
% Random data
x = rand(50,1);
y = .2 + x - x.^2 + randn(size(x))/10;
scatter(x,y)
xx = linspace(0,1);
% Unconstrained fit
b1 = polyfit(x,y,3);
line(xx,polyval(b1,xx),'color','r');
% Constrained to pass through the origin
b2 = [x.^3 x.^2 x]\y;
line(xx,polyval([b2; 0],xx),'color','g');
I hope this is clear enough. Let me know if not.
-- Tom

Don't use polyfit. The multilinear regression can be solved with one of
the division operators (I can't remember which), but the simple model is

Y=AX

When using an intercept, you catenate a column of ones onto X. Just
don't catenate the column of ones, and A=Y/X (or Y\X-- I can never get
this straight).

Scott

[Greg Heath]

"Scott Seidman" <***@mindspring.com> wrote in message news:***@130.133.1.4...
-----SNIP

Post by Scott Seidman
Y=AX
When using an intercept, you catenate a column of ones onto X. Just
don't catenate the column of ones, and A=Y/X (or Y\X-- I can never > get

this straight).

I interpret Y/X as Y post-divided by X (the X is "under" the Y). Therefore

(Y/X)*X = Y.

I interpret Y\X as X pre-divided by Y (the Y is "under" the X). Therefore

Y*(Y\X) = X.

Hope this helps.

Greg

[Peter Spellucci]

Post by Steven Lord

Post by Seth Gilchrist
I want to use polyfit in a script file and force the y-intercept of
the fit line to zero. Is there a way to do this besides setting the
(n+1) value generated by p = polyfit(x,y,n) to zero manually?
Thanks for the help
Seth

Well, if the y-intercept of the line is zero (i.e. the line goes through the
y = a_n*x^n+a_(n-1)*x^(n-1)...+a_1*x
y/x = a_n*x^(n-1)+a_(n-1)*x^(n-2)...+a_2*x+a_1
which is an (n-1)st degree polynomial relating x to y/x. [In the case where
x is zero, you know y should be zero so replace any elements of (y/x) where
x==0 by 0.] Now you can polyfit x and (y/x) with an n-1 degree polynomial
and shift all the coefficients left by 1 place.
--
Steve Lord

this is formally correct but introduces much "bias" in the statistics (for
example the confidence intervals for the coefficients would be completely
different). why not simply
a=[x,x.^2,x.^3,.....,x.^n]\y;
???

hth
peter