From help-octave-request at bevo dot che dot wisc dot edu  Tue Dec 15 16:59:46 1998
Subject: Re: fit data - trendline through zero
From: Eric Ortega <eortega at seas dot upenn dot edu>
To: help-octave at bevo dot che dot wisc dot edu
Date: Tue, 15 Dec 1998 17:59:44 -0500 (EST)


On Tue, 15 Dec 1998 lash at tellabs dot com wrote:

> > 
> > Marcel -
> >  
> > > I use polyfit to find the trendline in my data. How can I force it to go
> > > through zero ?
> > 
> > Divide your data by x before fitting.  Think about it.
> > 
> > On a related note, I have a hacked polyfit (polyweightfit) that
> > accepts weighting factors for each point.  If that sounds generally
> > interesting, I can post or submit it.
> > 
> >      - Larry Doolittle   ldoolitt at jlab dot org
> > 
> 
> I originally thought this as well, but I don't think that the result is
> then truly least squared error.  The error for larger values of x will
> be weighted less since it is divided by x. For example, if
> 
> x=[1,1000000,1000001,1000002]'
> and 
> y=x+0.5
> 
> Using polyfit(x,y./x,0) gives a result of 1.1250 which has pretty large
> squared error.  The actual answer should be very close to 1.
> 
> After thinking a bit, and looking at polyfit.m, I think that if you
> just want to fit y=mx instead of y=mx+b you would want to use the following:
> 
> If x and y are column vectors,
> 
> m = (x'*x)/(x'*y)
> 
> 
> If you wanted to fit y=p[1]*x^2+p[2]x
> 
> X = [x,x.^2];
> p = (X'*X)/(X'*y)
> 
> Someone correct me if I am wrong.
> 
> Bill Lash
> lash at tellabs dot com
> 
> 


If you want to fit `y=mx' and not `y=mx+b', then the easiest way to go
about this, IMHO, is to move your puff of points to straddle the origin,
accomplished as follows:
 
  tx = (max(size(x))^(-1)*(sum(x))  % x average
  x_h = x - tx;
  ty = (max(size(x))^(-1)*(sum(y))  % y average
  y_h = y - ty;

This works with Matlab 5.0, I don't know about octave.
Notice that all you are doing is subtracting the average value (in
each direction) from your set of points.

At this point, the `m' in `y=mx':

  m = sum(x_h'*y_h)/sum(x_h'*x_h);

And you're done.



					adios,
			 		  eo