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

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$?

• I believe the answer to this question requires a clearer explanation of what exactly the "$dy$" values represent. (I suspect their squares might be proportional to the diagonal elements of $M$, but there's not enough information here to be sure.) Could you tell us about the $dy$? – whuber Dec 23 '16 at 15:26
• @whuber, The $dy$ values are the error estimates for the dependent variables, $y$. I just updated the question to reflect this. – Steven C. Howell Dec 23 '16 at 15:29

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.

• The trick of "divide observations and predictors by error" (i.e. multiply by weight) works very generally as far as I know. For example the Wikipedia link by the OP shows the problem being solved is $\min_\beta\|W^{1/2}(y-X\beta)\|^2$, which shown the $\beta$ will be unchanged by multiplying the $W^{1/2}$ through. However note that in standard regression, not just $x$ but also "$1$" are predictors, i.e. $X=[x,1]$. I am not sure if your formulas account for the $1/dy_i$ "predictor". – GeoMatt22 Dec 23 '16 at 17:20
• Although this appears to be a correct description of one way to carry out weighted regression (modulo the caveat by @GeoMatt), it doesn't seem likely to solve your problem. The issue is that the $dy$ are only part of the error variance: there is another component reflecting other random deviations of the data from the idealized linear function you are fitting. Weighting by $1/dy$ in this fashion can grossly inflate the influence of data with small $dy$. To see why, consider an extreme case where $y$ is measured perfectly precisely, so $dy=0$. Your solution will force the line through $(x,y)$! – whuber Dec 23 '16 at 22:48
• You could perhaps introduce a "minimum uncertainty", $\sigma_0$, and then use $\sigma_i^2=dy_i^2+\sigma_0^2$ to compute the weights (a simple shrinkage estimate). – GeoMatt22 Dec 23 '16 at 23:03
• @whuber I see your point but the error on the type of data i am considering is actually quite large. I am actually trying to fully propagate this large error, and the variance from a straight line, doing both of these in an accurate manner so I can make a case that the fit values are not as certain as people typically treat them. – Steven C. Howell Dec 24 '16 at 1:42

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.