While this does not answer the question asked (and will therefore not be accepted), I provide this response in hopes it will help others and promote worthwhile discussion.

In his book [Experimentation: An Introduction to Measurement Theory and Experiment Design](https://books.google.com/books?id=LicvAAAAIAAJ&dq=editions%3AD4pkF1UZ6JsC&q=least+squares+fit#v=onepage&q&f=false), David Baird provides a simple explanation for doing a linear least squares fit using the diffences between the measured and fit values to estimate the error of the fit parameters.  The best fit for the parameters $m$ and $b$ in $$y=mx+b$$ will is determined using his eqn. (6.3):

$$m = \frac{n\sum x_i y_i - \sum x_i \sum y_i}{n\sum(x^2_i) - (\sum x_i)^2} $$

$$b = \frac{\sum(x_i^2)\sum y_i - \sum x_i \sum x_iy_i}{n\sum(x^2_i) - (\sum x_i)^2} $$

After obtaining $m$ and $b$, a standard deviation for the fit parameters can be obtained by calculating the differences of each $y_i$ value from the fit, $\delta y_i = y_i - (m x_i +b)$.  From these differences, one calculates the standard deviation of the data from the fit line using:

$$\sigma_y = \sqrt{\frac{\sum(\delta y_i)^2}{n-2}}$$

and then the standard deviation of the  parameters using
$$\sigma_m = \sigma_y \sqrt{\frac{n}{n\sum{x_i^2}-\left(\sum{x_i}\right)^2}}$$
and
$$\sigma_b = \sigma_y \sqrt{\frac{\sum{x_i^2}}{n\sum{x_i^2}-\left(\sum{x_i}\right)^2}}$$


To include the measurement error, $dy_i$, in the fit one would divide the initial system of equations by $dy_i$ giving $$\frac{y_i}{dy_i} = m\frac{x_i}{y_i} + b \frac{1}{dy_i}$$ then repeat Baird's derivation to get the weighted fit parameters

$$m = \frac{\sum \frac{1}{dy_i}\sum \frac{x_i y_i}{dy_i^2} - \sum \frac{x_i}{dy_i^2} \sum \frac{y_i}{dy_i}}{\sum \frac{1}{dy_i}\sum \frac{x^2_i}{dy_i^2} - \sum \frac{x_i}{dy_i}\sum \frac{x_i}{dy_i^2}} $$

$$b = \frac{\sum \frac{x_i^2}{dy_i^2}\sum \frac{y_i}{dy_i} - \sum \frac{x_i}{dy_i} \sum \frac{x_iy_i}{dy_i^2}}{\sum \frac{1}{dy_i}\sum \frac{x^2_i}{dy_i^2} - \sum \frac{x_i}{dy_i}\sum \frac{x_i}{dy_i^2}} $$

Notice that the $b\frac{1}{dy_i}$ term makes it so you cannot simply divide $x_i$ and $y_i$ by $dy_i$ (as pointed out in the comments below).

Unfortunately, this does not propagate the measurement error into an error in the fit parameters though.