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

  • 2
    $\begingroup$ 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$? $\endgroup$ – whuber Dec 23 '16 at 15:26
  • 2
    $\begingroup$ @whuber, The $dy$ values are the error estimates for the dependent variables, $y$. I just updated the question to reflect this. $\endgroup$ – 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.

| cite | improve this answer | |
  • 3
    $\begingroup$ 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". $\endgroup$ – GeoMatt22 Dec 23 '16 at 17:20
  • $\begingroup$ 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)$! $\endgroup$ – whuber Dec 23 '16 at 22:48
  • $\begingroup$ 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). $\endgroup$ – GeoMatt22 Dec 23 '16 at 23:03
  • 1
    $\begingroup$ @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. $\endgroup$ – 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.

| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.