5
$\begingroup$

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?

$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$ Jun 16, 2013 at 21:47
  • $\begingroup$ @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. $\endgroup$ Jun 17, 2013 at 0:39

1 Answer 1

3
$\begingroup$

In your linked post DrBubbles seemed to have the same issue. I'm not sure about the underlying theory and the link he posted to seems to be proprietary.

Your error propagation seemed fine and I was getting the same thing you were.

I guess what I would do would be to take the root sum square of the two errors and use that as an effective y error. Then find the weights with 1/yErrEffective^2.

I don't see why this method wouldn't be fine for a linear regression, however for regressions with more complex curvature I'm sure it could have some issues. I'm not sure how he was finding his confidence bands either, but they are non-ideal because they do not expand when outside the data range, as the better model does.

Determine the weights:

errs = {{10, 200000/841}, {4, 3125/49}, {3, 2000/81}};
wts = 1/(Part[Transpose[errs], 1]^2 + Part[Transpose[errs], 1]^2)

Do the regression:

data = {{327, 100000/29}, {153, 12500/7}, {88, 10000/9}};
Needs["HypothesisTesting`"];
Needs["LinearRegression`"];
markerSize = 13;
title = "Linear Regression Title";
fontSize = 20;
imSize = 800;
xLabel = "T [s]";
yLabel = "1/\[Rho]";

xMax = Max[
   Part[Transpose[data], 
    1]] + .1*(Max[Part[Transpose[data], 1]] - 
     Min[Part[Transpose[data], 1]]); xMin = 
 Min[Part[Transpose[data], 
    1]] - .1*(Max[Part[Transpose[data], 1]] - 
     Min[Part[Transpose[data], 1]]); fit = 
 NonlinearModelFit[data, a + b*x, {a, b}, x, Weights -> wts];

Framed[Column[{Row[{Show[
      Plot[fit["SinglePredictionBands"], {x, xMin, xMax}, 
       PlotStyle -> Dotted, PlotRange -> {{xMin, xMax}, Automatic}, 
       PlotLegends -> {"95% Single Prediction Bands"}, 
       ImageSize -> imSize, PlotLabel -> title, 
       FrameLabel -> {xLabel, yLabel}, 
       BaseStyle -> {FontSize -> fontSize}, Frame -> True, 
       AxesOrigin -> {0, 0}, LabelStyle -> Black], 
      Plot[fit["MeanPredictionBands"], {x, xMin, xMax}, 
       PlotStyle -> Dashed, 
       PlotLegends -> {"95% Mean Prediction Bands"}], 
      ListPlot[data, 
       PlotMarkers -> {"\[FilledSmallCircle]", markerSize}, 
       PlotStyle -> Black], 
      Plot[fit["BestFit"], {x, xMin, xMax}, 
       PlotLegends -> {"Best Fit"}]]}], 
   Row[{Grid[{{"Fit Model: ", fit["BestFit"], "RSquared:", 
        fit["RSquared"]}, {"Parameters", 
        fit["BestFitParameters"] // TableForm, 
        "Parameter 95% CI", {\[PlusMinus](Abs[
               Part[Part[fit["ParameterConfidenceIntervals"], 1], 1] -
                 Part[Part[fit["ParameterConfidenceIntervals"], 1], 
                 2]]/2) // 
           N, \[PlusMinus](Abs[
                Part[Part[fit["ParameterConfidenceIntervals"], 2], 
                  1] - Part[
                  Part[fit["ParameterConfidenceIntervals"], 2], 2]]/
               2) // N // TableForm} // TableForm}}, 
      Alignment -> Left, Spacings -> {2, 1}, Frame -> All, 
      ItemStyle -> "Text", 
      Background -> {{LightGray, None, LightGray, None}}]}]}], 
 ImageMargins -> 10, FrameMargins -> 10]

I'm assuming you're obviously trying to get parameters from the regression and this has the added bonus of giving a 95% CI.

Hope this helps.

Edit: Added image. enter image description here

$\endgroup$
3
  • $\begingroup$ sorry, it's been to long since I worked on this project and therefore I do not remember much of it, so I cannot answer if it's ok to do so or not. thank you for the answer though. $\endgroup$
    – TheoryX
    Aug 16, 2016 at 15:14
  • $\begingroup$ Shouldn't the wts be wts = 1/(Part[...., 1]^2 + Part[...., 2]^2), instead of taking Part[...., 1] twice? $\endgroup$
    – corey979
    May 2, 2018 at 14:27
  • $\begingroup$ I'm wondering if WFR: OrthogonalLineFit might be useful. $\endgroup$
    – Sterling
    Mar 26, 2021 at 18:05

Your Answer

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

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