Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I was informing myself about the least square algorithms with errors in x and y. I found this post and the top answer wasn't working for magnitudes around 10^-20. Because I couldn't figure out why, I completed the code from belisarius to also produce errors on slope and intercept. But my code is really, really slow, if I take 10 or more {x, y} pairs with errors it takes way too long. I think it's the Reduce part that slows it down so much, but I don't know how to make it more efficient. Also the function output isn't really elegant yet.

LinFit[xi_List, yi_List, errx_List, erry_List]:=
  Module[{n=Length@xi,wi,ui,vi,wmean,d,g,a,b,set,least,c,wxi,wyi},
  wxi=errx^-2;
  wyi=erry^-2;
  wi[i_,m_]:=wxi[[i]] wyi[[i]]/(m^2 wyi[[i]]+wxi[[i]]);
  ui[i_,m_]:=xi[[i]]-wmean[xi,m];
  vi[i_,m_]:=yi[[i]]-wmean[yi,m];
  wmean[q_List,m_]:=Sum[wi[i,m] q[[i]],{i,n}]/Sum[wi[i,m],{i,n}];
  d[m_]:=Sum[wi[i,m]^2 ui[i,m]^2/wxi[[i]],{i,n}];
  g[m_]:=-Sum[wi[i,m] ui[i,m] vi[i,m],{i,n}]/d[m];
  a[m_]:=2 Sum[wi[i,m]^2 ui[i,m] vi[i,m]/wxi[[i]],{i,n}]/(3 d[m]);
  b[m_]:=(Sum[wi[i,m]^2 vi[i,m]^2/wxi[[i]],{i,n}]-Sum[wi[i,m] ui[i,m]^2{i,n}])/(3d[m]);

  set={ToRules@Reduce[\[FormalM]^3-3 a[\[FormalM]] \[FormalM] \[FormalM]+
        3 b[\[FormalM]] \[FormalM]-g[\[FormalM]]==0&&\[FormalC]==wmean[yi,\[FormalM]]-
        \[FormalM]wmean[xi,\[FormalM]]&&\[FormalA]==Sqrt[1/(n-2) Sum[wi[i,\[FormalM]]
        ( \[FormalM] ui[i,\[FormalM]]-vi[i,\[FormalM]])^2,{i,n}]/
        Sum[wi[i,\[FormalM]] ui[i,\[FormalM]]^2,{i,n}]]&&\[FormalB]==Sqrt[(
        Sum[wi[i,\[FormalM]] xi[[i]]^2,{i,n}]/Sum[wi[i,\[FormalM]],{i,n}])*
        \[FormalA]^2],{\[FormalM],\[FormalC],\[FormalA],\[FormalB]}, 
        Backsubstitution->True]};
  least=Sum[wxi[[i]] (xi[[i]]-(yi[[i]]-\[FormalC])/\[FormalM])^2+wyi[[i]] (yi[[i]]-
           (\[FormalM] xi[[i]]+\[FormalC]))^2,{i,Length@xi}]/.set[[Flatten@Position[
           \[FormalM]/.set,_Real]]];
  c=Flatten@set[[Flatten@Position[\[FormalM]/.set,_Real]]][[Position[
      least,Min@least][[1]]]];
  {Function[(\[FormalM]/.c[[1]])#+\[FormalC]/.c[[2]]][x],
      {\[FormalA]/.c[[3]],\[FormalB]/.c[[4]]}}
]
share|improve this question
add comment

1 Answer

up vote 0 down vote accepted

"worals" is an acronym for Weighted Orthogonal Regression by Alternating Least Squares.
Arguments: x and y are lists of measured values,
sx and sy are lists of the corresponding standard errors of measurement.
The returned values are {chisquare, {intercept, slope}}.

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

If the true relation is linear, and if the errors are independent normal with zero means and standard deviations as given in sx & sy, (or, more to the point, if the foregoing are not too far from truth) then the returned chisquare will have a Chi-Square distribution with n-2 degrees of freedom, where n = Length@x. The corresponding p-value is GammaRegularized[(n-2)/2, chisquare/2].

EDIT, responding to questions.

The model is {x = z + d, y = a + b*z + e}, where d & e are random errors, and a, b, & z are unknowns for which values are to be found that minimize the weighted sum of squares f = #.# & @ Join[(z-x)/sx,(a+b*z-y)/sy]. We start with z = x, then minimize f alternately with respect to either {a,b} or z with the other held constant, until f no longer changes.

To estimate the error in {a,b}, I usually jackknife the solution:

n = Length@x; {f,ab} = worals[x,y,sx,sy];
{jab,jc} = {ab + (n-1)(ab-Mean@#), (n-1)^2/n Covariance@#}& @
           Table[Last[worals@@(Delete[#,i]&/@{x,y,sx,sy})],{i,n}]

jab is the jackknifed estimate of {a,b}, and jc is the jackknifed estimate of its covariance matrix; Sqrt@Diagonal@jc gives the estimates of the two standard errors.

share|improve this answer
    
Thanks, but I have still trouble understanding how this algorithm works. Could you tell me how i can evaluate the error on slope and intercept in a similar way. –  Gebbo Oct 1 '13 at 19:12
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.