How to calculate uncertainty of linear regression slope based on data uncertainty (possibly in Excel/Mathematica)?

Example: Example plot Let's have data points (0,0), (1,2), (2,4), (3,6), (4,8), ... (8, 16), but each y value has an uncertainty of 4. Most functions I found would calculate the uncertainty as 0, as the points perfectly match the function y=2x. But, as shown on the picture, y=x/2 match the points as well. It's an exaggerated example, but I hope it shows what I need.

EDIT: If I try to explain a bit more, while every point in example has a certain value of y, we pretend we don't know if it's true. For example the first point (0,0) could actually be (0,6) or (0,-6) or anything in between. I'm asking if there is an algorithm in any of the popular problems that takes this in account. In the example the points (0,6), (1,6.5), (2,7), (3,7.5), (4,8), ... (8, 10) still fall in the uncertainty range, so they might be the right points and the line that connects those points has an equation: y = x/2 + 6, while the equation we get from not factoring in the uncertainties has equation: y=2x + 0. So uncertainty of k is 1,5 and of n is 6.

TL;DR: In the picture, there is a line y=2x that's calculated using least square fit and it fits the data perfectly. I'm trying to find how much k and n in y=kx + n can change but still fit the data if we know uncertainty in y values. In my example, uncertainty of k is 1.5 and in n it's 6. In the image there is the 'best' fit line and a line that just barely fits the points.

  • 1
    $\begingroup$ If you have uncertainty in your $x$'s, you generally shouldn't use ordinary linear regression because it's biased (though small uncertainties will result in small bias; maybe you don't care so much). Are the uncertainties in your y's always constant or do they vary? How are you fitting your line? $\endgroup$ – Glen_b Sep 21 '13 at 0:38
  • $\begingroup$ I'm mainly asking about uncertainty in y's. But I'd be happy if the solution would consider uncertainty in x's as well. Generally they aren't constant, but a solution that requires them to be constant would be fine as well, I tried fitting with linearfit, fit, findfit in mathematica and linest (and a custom function that I found that weights values by uncertainties) in excel. $\endgroup$ – bedanec Sep 21 '13 at 0:42
  • $\begingroup$ Can you explain as clearly as possible how the "uncertainty" values would be related to what would happen under say a repeat of the experiment? i.e. what, precisely, do these uncertainties represent? $\endgroup$ – Glen_b Sep 21 '13 at 10:56
  • $\begingroup$ Err for example if y is weight, but the digital scale is only accurate to +-6. (Not the best example with these values, but for example when scale shows 255g it could be 255.0 or 255.9) $\endgroup$ – bedanec Sep 21 '13 at 11:36
  • $\begingroup$ Could you explain how that would relate to the first sentence of my question please? $\endgroup$ – Glen_b Sep 22 '13 at 1:51

Responding to "I'm trying to find how much $k$ and $n$ in $y = k x + n$ can change but still fit the data if we know uncertainty in $y$ values."

If the true relation is linear and the errors in $y$ are independent normal random variables with zero means and known standard deviations then the $100(1-\alpha)$% confidence region for $(k,n)$ is the ellipse for which $\sum (k x_i + n - y_i)^2/\sigma_i^2 < \chi_{d,\alpha}^2$, where $\sigma_i$ is the standard deviation of the error in $y_i$, $d$ is the number of $(x,y)$ pairs, and $\chi_{d,\alpha}^2$ is the upper $\alpha$ fractile of the chi-square distribution with $d$ degrees of freedom.

EDIT - Taking the standard error of each $y_i$ to be 3 -- i.e., taking the error bars in the plot to represent approximate 95% confidence intervals for each $y_i$ separately -- the equation for the boundary of the 95% confidence region for $(k,n)$ is $204 (k-2)^2 + 72n(k-2) + 9n^2 = 152.271$.

enter image description here

| cite | improve this answer | |

I did a naive direct sampling with this simple code in Python:

import random
import numpy as np
import pylab
def uncreg(x, y, xu, yu, N=100000):
    out = np.zeros((N, 2))
    for n in xrange(N):
        tx = [s+random.uniform(-xu, xu) for s in x]
        ty = [s+random.uniform(-yu, yu) for s in y]
        a, b = np.linalg.lstsq(np.vstack([tx, np.ones(len(x))]).T, ty)[0]
        out[n, 0:2] = [a, b]
    return out
if __name__ == "__main__":
    P = uncreg(np.arange(0, 8.01), np.arange(0, 16.01, 2), 0.1, 6.)
    H, xedges, yedges = np.histogram2d(P[:, 0], P[:, 1], bins=(50, 50))
    pylab.imshow(H, interpolation='nearest', origin='low', aspect='auto',
                 extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]])

and got this: Histogram of direct sampling results

Of course you can mine the P for the data you want, or change the uncertainty distributions.

| cite | improve this answer | |

I was on the same hunt before and I think this may be a useful place to start. The excel macro function gives linear fit terms and their uncertainties based on tabular points and uncertainty for each point in both ordinates. Maybe look up the paper it is based on to decide if you want to implement it in a different environment, modify, etc. (There is some legwork done for Mathematica.) It seems to have good walk-through documentation on the surface but haven't opened up the macro to see how well annotated it is.

| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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