MathGroup Archive 2009

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: linear regression with errors in both variables

On Feb 10, 2:56 am, Joerg <scha... at> wrote:
> Hi,
> I want to test the hypothesis that my data
> follows a known simple linear relationship,
> y = a + bx. However, I have (known) measurements
> errors in both the y and the x values.
> I suppose just a simple linear regression
> does not do here.
> Any suggestions how do test this correctly?
> Thanks,
> joerg

You don't say whether the errors are the same for all data points.
If they are, and sx & sy are the (known) standard deviations of the
measurement errors in x & y, then

  ClearAll[a,b]; Minimize[Evaluate[Expand[
  #.#&[a + b*x - y]]/(b^2 sx^2 + sy^2)], {a,b}]

will give what you want. You can save a little time if you eliminate
'a' by centering the data. And if Mathematica honors the parentheses,
centering will also reduce potential roundoff error problems.

  xbar = Mean@x; ybar = Mean@y; ClearAll[a,b];
  {#[[1]], {a -> ybar - #[[2,1,2]]*xbar, #[[2,1]]}}& @
  #.#&[b(x-xbar)-(y-ybar)]]/(b^2 sx^2 + sy^2)], b]

If the errors differ from point to point then sx & sy will be lists.

  Minimize[Tr[ (a + b*x - y)^2 / (b^2 sx^2 + sy^2) ], {a,b}]

will work, but it's slow. Pre-evaluating doesn't help.
The following is faster:

  worals[x_,y_,sx_,sy_] := Block[
  {a,b,f,z, u = 1/sx, v = 1/sy, w = (sy/sx)^2},
  {a,b} = (y*v).PseudoInverse@{v,x*v}; f = #.#&[(a+b*x-y)v];
  While[f > (z = (x*w + (y-a)b)/(b^2 + w);
             {a,b} = (y*v).PseudoInverse@{v,z*v};
             f = #.#&@Join[(z-x)u,(a+b*z-y)v])]; {f,{a,b}}]

('worals' is an acronym for Weighted Orthogonal Regression by
Alternating Least Squares.)

  • Prev by Date: Re: Mathematicas simplifications
  • Next by Date: Re: Custom label for a slider in Manipulate
  • Previous by thread: Re: linear regression with errors in both variables
  • Next by thread: Re: linear regression with errors in both variables