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]:=
  wi[i_,m_]:=wxi[[i]] wyi[[i]]/(m^2 wyi[[i]]+wxi[[i]]);
  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}])*
  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[
  • $\begingroup$ Without formatting this is rather unreadable; in the present form I would suggest to ask a machine, not a human being... :) $\endgroup$
    – gwr
    Commented Sep 18, 2016 at 11:07

1 Answer 1


"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])];

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@#}& @

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.

  • $\begingroup$ 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. $\endgroup$
    – NicolasW
    Commented Oct 1, 2013 at 19:12
  • $\begingroup$ There is a syntax error in the While statement, which needs test and body separated by a comma. $\endgroup$
    – gwr
    Commented Sep 18, 2016 at 10:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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