From help-octave-request at bevo dot che dot wisc dot edu  Tue Dec 15 20:37:43 1998
Subject: Re: fit data - trendline through zero
From: Joao Cardoso <jcardoso at inescn dot pt>
To: Eric Ortega <eortega at seas dot upenn dot edu>
CC: help-octave at bevo dot che dot wisc dot edu
Date: Wed, 16 Dec 1998 02:36:50 +0000

Eric Ortega wrote:
> 
> 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
> 
> 

why not use: 

  b = pinv(x)*y

where x and y are column vectors? try "help ols"

Joao
-- 
Joao Cardoso, INESC  |  e-mail:
R. Jose Falcao 110   |  tel:    + 351 2 2094345
4050 Porto, Portugal |  fax:    + 351 2 2008487