# Linear Fit when Data has Uncertainty

I am attempting to find the slope and y-intercept (along with their uncertainty) from a set of data. In this case, I am graphing Gamma Energy (MeV) vs. Peak Centroid (Channel). Here is my data:

Gamma Energy (MeV): 1.17, 1.33, 0.032, 0.662, 0.511, 1.275, 0.088

Peak Centroid (Channel): 622.65, 712, 21, 360.38, 280.64, 676.85, 18.68

Peak Centroid Uncertainty: 0.0342, 0.0347, 2, 0.0155, 0.0231, 0.0288, 0.1346

As can be seen, one Centroid Uncertainty value is significantly higher than the rest, so I cannot use a normal non-weighted least squares fit. I researched this a bit and found that I can use the formula $\left ( X^TWX \right )\hat{\beta }=X^TWy$ to solve for $\hat{\beta}$. Now, I know $y_{i}$ is simply each of my Gamma Energies. I also know $W_{ii}=\frac{1}{\sigma_{i}^2}$ where $\sigma_{i}$ corresponds to each of my Centroid Uncertainties. All other elements in $W$ are $0$. My issue is, I do not know how to contruct $X$. I have looked at various sources, and none of them really made sense to me. What exactly does each element correspond to, and how do I implement that in my case. When I find $\hat{\beta}$, won't this simply correspond to the $m$ and $b$ values in $y=mx+b$? How can I find the uncertainty in $m$ and $b$? Also, in my case, the uncertainties are in the x-values where (I think) the formula I am using assumes them to be in the y-values. Will this affect my results? As an aside, if I had errors in both x and y-values, how would I solve that?

• the $\bf X$ and $\bf Y$ to use in the matrix formula are the respective deviations from the mean: the line passes through $(X_m, Y_m)$. Oct 20, 2017 at 17:04
• and, If the uncertainties are on the $X$, but not on the $Y$, just reverse your line as $(x-x_m)=q(y-y_m)$ Oct 20, 2017 at 17:12
• @GCab I'm sorry, but I don't follow. So you are saying I take the mean of y along with mean of x (with uncertainties figured in), and those are my elements? Because from what I read, y is a vector while X is a nxn+1 matrix. So what are all of the elements in X? Oct 20, 2017 at 18:54
• For what I know, both $\bf X$ and $\bf Y$ are vectors (while $\bf W$ is a diagonal matrix), otherwise the matrix equation you cite would be inconsistent. The vector $\bf X$ would be given by $(x_k-x_m)\; |\,k=1, \cdots, n$. Similar for $\bf Y$. Note that the least square method assumes that the $x$ be "certain", and the the $y$ subject to error, following the same centered normal distribution.If your situation is the reverse, then you have to reverse the role of $x$ and $y$, and look for the regression $$(x-x_m)=p \cdot (y-y_m)$$. Oct 20, 2017 at 22:29
• To be more helpful please cite the reference you found. Oct 20, 2017 at 22:29

1) The least square method, in its simplest standard form, assumes that the $x$ be "certain", and the the $y$ subject to error, following the same centered normal distribution, i.e. with the same $\sigma$ (homoscedasticity).

You instead want to find the reverse relation: then, once you have got $y=mx+b$, you can always reverse it into $x=(1/m)y-b/m$.

Trying to avoid confusion let's follow the standard way, understanding that the $x$ refer to the error-free variable, while $y$ to that subject to error.
So we are looking for the regression line $y=mx+b$.

2) We assume to have $n$ independent random variables $Y_k$, whose PDF (probability density function) is $$\bbox[lightyellow] { p\left( {\eta _k } \right) = {1 \over {\sqrt {2\pi \sigma _k ^2 } }}\exp \left( { - {1 \over 2}{{\left( {\eta _k - \hat m\,x_k - \hat b} \right)^2 } \over {\sigma _k ^2 }}} \right) } \tag{1}$$ where $\hat y = \hat m\,x_k + \hat b$ is the supposed "exact" value of $y$, defined by the the "exact" values of the parameters $m$ and $b$ indicated with the hat.

3) The probability that we get the $n$-uple of values $(y_1,y_2,\cdots,y_n)$, considering the $x$ and $\sigma$ $n$-tuples as fixed, is $$\bbox[lightyellow] { p\left( {y_1 } \right) \cdots p\left( {y_n } \right)\;dy_1 \cdots dy_n } \tag{2}$$ which is maximum when $$\bbox[lightyellow] { E^2 (m,b) = \sum\limits_k {{{\left( {y_k - \hat m\,x_k - \hat b} \right)^2 } \over {\sigma _k ^2 }}} = \sum\limits_k {\left( {{{y_k - \hat y_k } \over {\sigma _k }}} \right)^2 } } \tag{2}$$ is minimum, which means $$\bbox[lightyellow] { \left\{ \matrix{ 0 = {\partial \over {\partial b}}E^2 (m,b)\quad \Rightarrow \quad \sum\limits_k {{{\left( {y_k - \hat m\,x_k - \hat b} \right)} \over {\sigma _k ^2 }}} = 0 \hfill \cr 0 = {\partial \over {\partial m}}E^2 (m,b)\quad \Rightarrow \quad \sum\limits_k {{{\left( {y_k - \hat m\,x_k - \hat b} \right)x_k } \over {\sigma _k ^2 }}} = 0 \hfill \cr} \right. }$$ The first tells us that the estimated value of $b$ (indicated with the tilde) is equal to this expression involving the weighted average of the $y_k$ and $x_k$ values, with weights $(1/\sigma_k^2)$ \bbox[lightyellow] { \eqalign{ & \tilde b = \sum\limits_k {\left( {{{\left( {1/\sigma _k ^2 } \right)} \over {\sum\limits_k {\left( {1/\sigma _k ^2 } \right)} }}} \right)y_k } - \tilde m\sum\limits_k {\left( {{{\left( {1/\sigma _k ^2 } \right)} \over {\sum\limits_k {\left( {1/\sigma _k ^2 } \right)} }}} \right)x_k } = \cr & = \bar y - \tilde m\,\bar x \cr} } \tag{3} So we can rewrite the second as $$\bbox[lightyellow] { 0 = \sum\limits_k {{{\left( {y_k - \hat m\,x_k - \hat b} \right)x_k } \over {\sigma _k ^2 }}} = \sum\limits_k {{{\left( {\left( {y_k - \bar y} \right) - \hat m\,\left( {x_k - \bar x} \right)} \right)\left( {x_k - \bar x} \right)} \over {\sigma _k ^2 }}} }$$ and get $$\bbox[lightyellow] { \tilde m = \sum\limits_k {{{\left( {y_k - \bar y} \right)\left( {x_k - \bar x} \right)} \over {\sigma _k ^2 }}} \;\;\mathop /\limits_{} \;\; \sum\limits_k {{{\left( {x_k - \bar x} \right)^2 } \over {\sigma _k ^2 }}} } \tag{4}$$

4) Passing to the more synthetic vectorial notation, indicating by $$\bbox[lightyellow] { {\bf Y} = \left( {\matrix{ {y_1 } \cr {y_2 } \cr \vdots \cr {y_n } \cr } } \right)\quad {\bf X} = \left( {\matrix{ 1 & {x_1 } \cr 1 & {x_2 } \cr \vdots & \vdots \cr 1 & {x_n } \cr } } \right)\quad {\bf \beta } = \left( {\matrix{ b \cr m \cr } } \right)\quad {\bf W} = \left( {\matrix{ {1/\sigma _1 ^2 } & 0 & \cdots & 0 \cr 0 & {1/\sigma _2 ^2 } & \cdots & 0 \cr \vdots & \vdots & \ddots & \vdots \cr 0 & 0 & \cdots & {1/\sigma _n ^2 } \cr } } \right) }$$ we are going to minimize $$\bbox[lightyellow] { E^2 ({\bf \beta }) = \left( {{\bf W}^{\,{\bf 1/2}} \,\left( {{\bf Y} - {\bf X\beta }} \right)} \right)^{\,{\bf T}} \left( {{\bf W}^{\,{\bf 1/2}} \,\left( {{\bf Y} - {\bf X\beta }} \right)} \right) }$$ which returns $$\bbox[lightyellow] { {\bf \tilde \beta } = \left( {{\bf X}^{\,{\bf T}} {\bf WX}} \right)^{ - \,{\bf 1}} {\bf X}^{\,{\bf T}} \,{\bf W}\;{\bf Y} } \tag{4}$$

Refer also to this Wikipedia article.

Applied to your data this gives $$\bbox[lightyellow] { b=13.27,\;m=522.76 }$$ .

5) Concerning the confidence abour the estimated values of $m$ and $b$, you can find various considerations in the already cited article as well as in this other related one, together with the necessary cautions as per the statistical soundness and applicability.

However, not pretending to be exhaustive nor totally rigorous, but just to outline the basis of a possible approach, let's go back to (1) and (2).
We can consider eq. (2) as proportional to the conditional probability of getting the vector $\bf Y$ , given $\bf \beta$ $$P\left( {{\bf Y}\,\,\left| {\;{\bf \beta }} \right.} \right) \propto p\left( {y_1 } \right) \cdots p\left( {y_n } \right)\;dy_1 \cdots dy_n$$ Assuming $P(\bf \beta)$ to be fairly constant, at least within the "admissible" range of values, then we can put that \bbox[lightyellow] { \eqalign{ & P\left( {{\bf \beta }\,\,\left| {\,{\bf Y}\;} \right.} \right) \propto p\left( {y_1 } \right) \cdots p\left( {y_n } \right)\;db\,dm \propto \cr & \propto \exp \left( { - {1 \over 2}\sum\limits_k {\left( {{{y_k - \hat m\,x_k - \hat b} \over {\sigma _k }}} \right)^2 } } \right) = \exp \left( { - {1 \over 2}\sum\limits_k {\left( {{{\,x_k } \over {\sigma _k }}\hat m + {1 \over {\sigma _k }}\hat b - {{y_k } \over {\sigma _k }}} \right)^2 } } \right) \cr} } \tag{5} with $x_k, y_k, \sigma_k$ known and fixed.
Let's pay attention to the fact that presuming that the $\sigma_k$ are known, and not to be determined by some estimation on the spreadness of the $y_k$, is a crucial point in fixing the statistics of the parameters (Normal instead of Student).

The sum of the quadratic forms in $b,m$ will provide a quadratic form that can be written as \bbox[lightyellow] { \eqalign{ & Q(b,m) = \sum\limits_k {\left( {{1 \over {\sigma _k }}\hat b + {{\,x_k } \over {\sigma _k }}\hat m - {{y_k } \over {\sigma _k }}} \right)^2 } = \cr & = \left( {T^2 \left( {\hat b - \tilde b} \right)^2 + 2R^2 \left( {\hat b - \tilde b} \right)\left( {\hat m - \tilde m} \right) + U^2 \left( {\hat m - \tilde m} \right)^2 + V^2 } \right) \cr} } \tag{6} where it is clear that for $\hat b = \tilde b,\;\hat m = \tilde m$ $Q(b,m)$ attains the minimum value, which is the condition defining the estimated values.

Taking the first and second derivatives it is possible to express $T,R,U,V$ wrt the known values. $$\bbox[lightyellow] { \left\{ \matrix{ T^{\,2} = \sum\limits_{1\, \le \,k\, \le \,n} {\left( {1/\sigma _{\,k} } \right)^{\,2} } \hfill \cr U^{\,2} = \sum\limits_{1\, \le \,k\, \le \,n} {\left( {x_{\,k} /\sigma _{\,k} } \right)^{\,2} } \hfill \cr R^{\,2} = \sum\limits_{1\, \le \,k\, \le \,n} {\left( {x_{\,k} /\sigma _{\,k} ^{\,2} } \right)} \hfill \cr V^{\,2} = \sum\limits_k {\left( {{1 \over {\sigma _k }}\tilde b + {{\,x_k } \over {\sigma _k }}\tilde m - {{y_k } \over {\sigma _k }}} \right)^2 } \hfill \cr} \right. }$$

We can translate the above into the equation of a Bivariate Normal Distribution, absorbing the term in $V$ into the normalizing factor. $$\bbox[lightyellow] { {1 \over {2\,\pi \,\sigma _{\,1} \,\sigma _{\,2} \sqrt {1 - \rho ^{\,2} } }}\exp \left( { - {1 \over {2\left( {1 - \rho ^{\,2} } \right)}}\left( {{{\left( {x_{\,1} - \mu _{\,1} } \right)^{\,2} } \over {\sigma _{\,1} ^{\,2} }} - 2\rho {{\left( {x_{\,1} - \mu _{\,1} } \right)\left( {x_{\,2} - \mu _{\,2} } \right)} \over {\sigma _{\,1} \,\sigma _{\,2} }} + {{\left( {x_{\,2} - \mu _{\,2} } \right)^{\,2} } \over {\sigma _{\,2} ^{\,2} }}} \right)} \right) }$$ where $\rho$ is the correlation $$\bbox[lightyellow] { \rho = {\rm corr}\left( {x_{\,1} ,x_{\,2} } \right) = {{V_{12} } \over {\sigma _{\,1} \,\sigma _{\,2} }} }$$

Finally we reach to $$\bbox[lightyellow] { \left\{ \matrix{ \rho = - {{R^{\,2} } \over {T\,U}} \approx - \,0.93396 \hfill \cr \sigma _{\,b} = {1 \over {T\sqrt {\left( {1 - \rho ^{\,2} } \right)} }} \approx {\rm 0}{\rm .029524} \hfill \cr \sigma _{\,m} = {1 \over {U\sqrt {\left( {1 - \rho ^{\,2} } \right)} }} \approx {\rm 0}{\rm .033654} \hfill \cr} \right. } \tag{7}$$ which result I could counter-check with my old CAS program.

From here you can easily find the ellipse which delimites the parameters at the required confidence.

Now you can reverse the relation $x->y$.

• I have verified that we have arrived at equivalent formulas. I have accepted your answer as it clearly is the better answer by explaining the process of arriving at the given equations. Oct 25, 2017 at 0:13
• @JavaNewbie: very good ! so we learnt something about the weighted regression together. Oct 25, 2017 at 11:23
• That we did! Thank you very much for all of the time and effort you put into this. I greatly appreciate it. Oct 25, 2017 at 11:29

After researching and receiving some direction from G Cab, I have found the following: $$m=\frac{\sum \frac{x_{i}}{\sigma_{i}^2}\cdot \sum \frac{y_{i}}{\sigma_{i}^2}-\sum \frac{x_{i}y_{i}}{\sigma_{i}^2}\cdot \sum \frac{1}{\sigma_{i}^2}}{\left ( \sum \frac{x_{i}}{\sigma_{i}^2} \right )^2-\sum \frac{x_{i}^2}{\sigma_{i}^2}\cdot \sum \frac{1}{\sigma_{i}^2}} \; \; \; \; \; \; \; \; \; \; \sigma_{m}=\sqrt{\frac{\sum \frac{1}{\sigma_{i}^2}}{\sum \frac{x_{i}^2}{\sigma_{i}^2}\cdot \sum \frac{1}{\sigma_{i}^2}-\left ( \sum \frac{x_{i}}{\sigma_{i}^2} \right )^2}}$$

$$b=\frac{\sum \frac{x_{i}}{\sigma_{i}^2}\cdot \sum \frac{x_{i}y_{i}}{\sigma_{i}^2}-\sum \frac{y_{i}}{\sigma_{i}^2}\cdot \sum \frac{x_{i}^2}{\sigma_{i}^2}}{\left ( \sum \frac{x_{i}}{\sigma_{i}^2} \right )^2-\sum \frac{x_{i}^2}{\sigma_{i}^2}\cdot \sum \frac{1}{\sigma_{i}^2}} \; \; \; \; \; \; \; \; \; \; \sigma_{b}=\sqrt{\frac{\sum \frac{x_{i}^2}{\sigma_{i}^2}}{\sum \frac{x_{i}^2}{\sigma_{i}^2}\cdot \sum \frac{1}{\sigma_{i}^2}-\left ( \sum \frac{x_{i}}{\sigma_{i}^2} \right )^2}}$$ This is for the case of unequal errors. If we have equal errors, this obviously reduces to its normal form. Please note that this assumes the error is in the y-values (dependent variable) while the x-values are error-free (dependent variable). In order to use this in my application, I had to fit the equation $x=qy+p$ where I used the equations above (switching $x_{i}$ and $y_{i}$ in the formulas). In this case, $p$ corresponds to $m$ and $q$ corresponds to $b$. After this, we rearrange the fit in the way we want it: $y=\frac{1}{q}x-\frac{p}{q}$. Thus, we arrive at

$$m=\frac{1}{q} \; \; \; \; \; \; \; \; \; \; \sigma_{m}=\sqrt{\sigma_{q}^2\left ( \frac{\partial m}{\partial q} \right )^2}=\frac{\sigma_{q}}{q}$$

$$b=-\frac{p}{q} \; \; \; \; \; \; \; \; \; \; \sigma_{b}=\sqrt{\sigma_{q}^2\left ( \frac{\partial b}{\partial q} \right )^2+\sigma_{p}^2\left ( \frac{\partial b}{\partial p} \right )^2}=\sqrt{\frac{\sigma_{p}^2}{q^2}+\frac{\sigma_{q}^2\cdot p^2}{q^4}}$$

• I posted my new answer, and then so your auto-reply: I did not check with mine. Oct 24, 2017 at 22:35