I have data $\vec x$ and targets $\vec y$. The targets are each uncertain in either direction up to a number $\delta y$. I also have a vector of weights $\vec w$.
As a simple example in R, let's just say
delta_y <- 1
x <- 1:10
y <- 1:10 + runif(10, -delta_y, delta_y)
w <- 1:10
Like so:
Now, I can do weighted linear regression with
fit <- lm(y ~ x, weights=w)
From this I can obtain the slope $\beta$ of this relation using
beta <- coef(fit)[2]
How can I propagate the uncertainty $\delta y$ into $\beta$?
I know I can do
summary(fit)
and obtain Std. Error for my slope, but is that really taking my uncertainties into account?
Best Answer
The first thing you have to ask is whether the error model is still Gaussian(-like). In other words, is it better for the model to be close to the observed $y_i$? Or does it not matter where in the interval $|y - y_i| \leq \delta y$ the model falls? Consider the following picture:
If you do a traditional least squares fit, which assumes a Gaussian error model, you get the green line. But if the error model is uniform inside the error bars, the red and teal lines are just as good as the green line. Under the uniform error model, there are an infinite number of solutions that are all equally good.
So for now, let's assume that the error model is Gaussian. Then the result is a constrained least squares problem that looks something like this: \begin{equation*} \begin{aligned} & \min_\beta & & \tfrac{1}{2}||W(A \beta - y)||_2^2 \\ & \text{subject to} & & |A \beta - y| \leq \delta y \end{aligned} \end{equation*} where $W = diag(w)$ and $A = [x\:e]$. For convenience of notation, I've assumed $\delta y$ is a vector. This actually adds to the flexibility of the model, allowing for a different $\delta y_i$ for each observation. If you multiply out the objective function and expand the constraints, you'll see this is just a quadratic program, which can be solved by your favorite QP solver: \begin{equation*} \begin{aligned} & \min_\beta & & \tfrac{1}{2}\beta^T A^T W^2 A \beta - y^T W^2 A \beta \\ & \text{subject to} & & \begin{bmatrix} A \\ -A \end{bmatrix} \beta \leq \begin{bmatrix} \delta y + y \\ \delta y - y \end{bmatrix} \end{aligned} \end{equation*} Note that this QP may be infeasable, as demonstrated by this image:
Finally, to obtain a covariance estimate for the fitted parameters... This really deserves its own topic because it's so complicated, but I'll give you a brief overview. First, solve the QP. This will identify the constraints active at the solution. Let the subscript $a$ denote the constraints active at the solution. Then form the Hessian of the Lagrangian using only the active constraints. It should look something like this:
\begin{equation*}
\begin{bmatrix}
A^T W^2 A & \begin{bmatrix} A \\ -A \end{bmatrix}_a^T \\ \begin{bmatrix} A \\ -A \end{bmatrix}_a & 0
\end{bmatrix}
\end{equation*}
Invert the Hessian and the principal submatrix is your covariance estimate. I would also do a Monte Carlo simulation to see if the covariance estimate is reasonable.z