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 have a question about finding a linear regression weighted with x and y errors. I found the answer provided by 0x4A4D in Estimate error on slope of linear regression given data with associated uncertainty, but didn't know how to post it under his post, so I'am putting it here. Sorry in advance if there is anything wrong with this.

If i use the following code provided by 0x4A4D:

ortlinfit[data_?MatrixQ, errs_?MatrixQ] := 
 Module[{n = Length[data], c, ct, dk, dm, k, m, p, s, st, ul, vl, w, wt, xm, ym},
        (* yes, I know I could have used FindFit[] for this... *)
        {ct, st, k} = Flatten[MapAt[Normalize[{1, #}] &, 
           NArgMin[Norm[Function[{x, y}, y - \[FormalM] x - \[FormalK]] @@@ data],
                   {\[FormalM], \[FormalK]}], 1]];
        (* find orthogonal regression coefficients *)
        {c, s, p} = FindArgMin[{
           Total[(data.{-\[FormalS], \[FormalC]} -
                 \[FormalP])^2/((errs^2).{\[FormalS]^2, \[FormalC]^2})],
           \[FormalC]^2 + \[FormalS]^2 == 1},
          {{\[FormalC], ct}, {\[FormalS], st}, {\[FormalP], k/ct}}];
        (* slope and intercept *)
        {m, k} = {s, p}/c;
        wt = 1/errs^2; w = (Times @@@ wt)/(wt.{1, m^2});
        {xm, ym} = w.data/Total[w];
        ul = data[[All, 1]] - xm; vl = data[[All, 2]] - ym;
        (* uncertainties in slope and intercept *)
        dm = w.(m ul - vl)^2/w.ul^2/(n - 2);
        dk = dm (w.data[[All, 1]]^2/Total[w]);
        {Function[\[FormalX], Evaluate[{m, k}.{\[FormalX], 1}]], Sqrt[{dm, dk}]}] /;
     Dimensions[data] === Dimensions[errs]

{lin, {sm, sk}} = ortlinfit[data, errs]

Show[
     Graphics[{AbsolutePointSize[4], Point[data], MapThread[Circle, {data, errs}]},
              Frame -> True], 
     Plot[{lin[x], lin[x] - sm x - sk, lin[x] + sm x + sk}, {x, -1, 9},
          PlotStyle -> {Directive[Thick, Red],
                        Directive[Dashed, Gray], Directive[Dashed, Gray]}]
     ]

there is a problem with some kind of data. I get my data like this: {x, xerr, y, yerr}

info={{327, 10, 29, 2}, {153, 4, 56, 2}, {88, 3, 90, 2}}

and change it to work with 0x4A4D's code for my problem to data = {x, 1/y} and errs = {xerr, yerr/y^2} (absolute errors which I want to plot od graph) with following code

Input:

data = {#[[1]], 1/#[[3]]} & /@ info
errs = {#[[2]], #[[4]]/#[[3]]^2} & /@ info

Output:

data = {{327, 1/29}, {153, 1/56}, {88, 1/90}}
errs = {{10, 2/841}, {4, 1/1568}, {3, 1/4050}}

If i use this data and errs it works perfectly (or at least so it looks at the first glance). Result is:

enter image description here

But for the sake of my report I have to multiply "y" and "yerr" with a 10^(-5) factor, so I get to work with following points:

Input:

info[[All,3]]*=10^(-5);
info[[All,4]]*=10^(-5);

data = {#[[1]], 1/#[[3]]} & /@ info;
errs = {#[[2]], #[[4]]/#[[3]]^2} & /@ info;

Output:

data = {{327, 100000/29}, {153, 12500/7}, {88, 10000/9}};
errs = {{10, 200000/841}, {4, 3125/49}, {3, 2000/81}};

And now his algorithm doesn't work as it should anymore since the result is this:

enter image description here

Does anyone possibly knows why?

share|improve this question
    
Is "weighted with x and y errors" the same thing as (weighted) "total least squares"? If so, there may be useful methods using a singular values decomposition. –  Daniel Lichtblau Jun 16 '13 at 21:47
    
@Daniel, well, there is indeed an SVD-based method in the unweighted TLS case, but I don't think I've seen something using SVD when weights are involved. –  J. M. Jun 17 '13 at 0:39
add comment

Know someone who can answer? Share a link to this question via email, Google+, Twitter, or Facebook.

Your Answer

 
discard

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

Browse other questions tagged or ask your own question.