incorporating error in best fit line equation I have two sets of data {X} and {Y} each element in those two sets has its own unique +/- error. I graphed x vs y in a scatter plot along with every point's error bar (both in x and y (x1 +/- a1 , y1 +/- b1)) and added a best fit line. the best fit line is linear and therefore has an equation of the form mx+b. My question is how do I calculate the error in the constants of the equation y=mx+b.
edit: the elements of each set contain unique uncertainty in both sides for example its x1=10+0.5 or x1=10-0.25 and not x1=+/-0.5. similarly with Y.
update:
Eiso
(10^52 erg s−1)
35.00 ± 7,
1.70 ± 0.9,
114.00 ± 9,
14.10 ± 2.8,
24.00 ± 2,
188.00 ± 10
. Epeak
(KeV)
464.10 ± 26,
26.01 ± 2,
1333.53 ± 28,
816.62 ± 88,
1125.98 ± 48,
1090.06 ± 54.
A sample of the actual data, keep in mind that these errors will be transformed in the log space. so the plotted data will actually be
y points = Log(Eiso) and the x points = log (Epeak/avg(Epeak)) which changes the error. Then a best fit line would be calculated and my main concern is with the error in the constants of y=mx+b.

Let $\mathbf X := (x_1, \ldots x_n)^t$ be the column vector containing all your $x$ values, similarly $\mathbf Y := (y_1,\ldots,y_n)^t$, and finally let $\beta := (m, b)^t$.
The least-squares (LS) method computes your values $m$ and $b$ as follows: Create the design matrix $\mathbf X_d$:
$$
\mathbf X_d := \left(\begin{matrix}
                      x_1 & 1 \\
                      x_2 & 1 \\
                      \vdots & \vdots \\
                      x_n & 1 \\
                     \end{matrix}\right).
$$
Then, our goal can be expressed as trying to find the solution $\beta$ for the equation:
$$
\mathbf X_d \beta = \mathbf Y.
$$
Next, LS computes the Moore-Penrose pseudo-inverse (MPPI) $\mathbf X_d^+$ of $\mathbf X_d$, defined as:
$$
\mathbf X^+ := (X_d^t X_d)^{-1} X_d^t
$$
(where we presume that $(X_d^t X_d)^{-1}$ has full rank, i.e. the $x_i$ are not all the same). Then, $\beta$ is obtained as:
$$
\beta = \mathbf X^+ \mathbf Y.
$$
So this is the very formula that LS uses to compute the coefficients $\beta = (m, b)^t$ from your data. I tell you all this, because, if you presume nonnegligible errors in your data $\mathbf X$, you would have to compute the condition number of $\mathbf X_d$:
$$
cond(\mathbf X_d) := \|\mathbf X_d\| \|\mathbf X_d^+\|.
$$
If this number is large, your LS is ill-conditioned and not reliable. I.e. the error propagated from your $\mathbf X$ to your solution vector $\beta$ could blow up.
Now, presume that the condition number is small and LS can be used. Then, to actually get some error bars for $\beta$ you could apply methods of error propagation, to see how the errors are amplified by the application of the MPPI.
A more general approach, that can also be used if the formulae are not as simple as for LS, is the application of the bootstrap method. Bootstrap computes, for a given sample, approximations of complex statistics. And the variance of your vector $\beta$ is such a statistic. Python provides a ready to use implementation, and for comprehensive R packages see e.g. here or here.

Edit: As @whuber points out in the comments, there is another problem with errors in the explanatory $\mathbf X$ data: If the error of the $x$ values is comparable with the overall standard deviation of the $x$, then one could get a reduced slope estimate.
