How do I calculate the variance-covariance matrix for a set of 2-D data points with errors: (x, y, dy)

Question

I am doing a weighted linear least squares fit to N measured values of the form $(x, y, dy)$, where $x$ is the independent variable, $y$ is the independent variable, and $dy$ is the error estimates for the independent variable. Following the linked wikipedia page, I am able to get the two fit parameters but not yet able to get the parameter errors. The article lists that the variance-covariance matrix for the parameters, $M^\beta$, can be calculated using $$M^\beta = (X^T W X)^{-1} W M W^T X(X^T W^T X)^{-1}\,,$$ where $X$ is an $n \times 2$ matrix of the $x$-values and ones, $W$ is an $n \times n$ diagonal matrix with $1/dy_i^2$ values on the diagonal, and $M$ is the variance-covariance matrix for the observations. How would I calculate $M$ from my $(x, y, dy)$ data?

This related question and this post use $$M_{jk} = \frac{1}{N-1} \sum_{i=1}^{N}(x_{ij} - \bar{x_j})(x_{ik} - \bar{x_k})$$ but should $dy$ somehow affect this? This is straightforward to calculate but I am not sure which index the $x$ or $y$ values should take (with this only being a 2$\times$2 symmetric matrix it only affects the diagonal).

The article mentions a major simplification when $W = M^{-1}$ but I am not sure if this is the case for my data. Is there a simple way to know this is the true without calculating $M$?

Answer 1

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, 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.

Answer 2

Looking at the formula for $M^\beta$, you can deduce the dimensions of $M$ should be $N \times N$, where $N$ is the number of data points. Therefore, the variance-covariance matrix is not comparing all the $x$'s to all the $y$'s, which would give a $2 \times 2$ matrix, but instead comparing the $N$ data points.

The wiki article explaining the covariance matrix states that, "Indeed, the entries on the diagonal of the covariance matrix $\Sigma$ are the variances of each element of the vector $\mathbf {X}$." If these are $dy_i^2$, and we can assume the errors between points are not related, then $M = W^{-1}$ and the simplification applies.