How to propagate uncertainties in weighted linear regression? 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?
 A: Perhaps it's my own biases for Bayesian methods speaking, but I think this is a prime example of how one might use Bayesian analysis to simultaneously account for uncertainties at different levels of your analysis. 
library(rstan)
noise_model <- '
    data{
        int         N;
        int       D;
        matrix[N,D] x;
        vector[N] y_obs;
        real delta_y;
        vector[N] weights;
    }
    parameters{
        vector[N] z;
        vector[D] beta;
    }
    transformed parameters{
        vector[N] y_hat;
        vector[N] y_star;
        y_hat <- x*beta;
        y_star <- 2*delta_y*inv_logit(z)+y_obs-delta_y;
    }
    model{
        y_star ~ normal(y_obs, something);
        y_star ~ normal(y_hat, weights);
    }
'

N <- 100
beta_true <- c(-2, 1.5)
x <- cbind(1, rnorm(N, mean=-0.5, sd=1.5))
y_obs <- x%*%beta_true+runif(N, -0.5, 0.5)

plot(x[,2],y_obs)

stan_data <- list(D=2, N=N, y_obs=as.vector(y_obs), x=x, delta_y=0.5)

fit <- stan(model_code=noise_model, data=stan_data)

str(extract(fit, "y_star"))

This model postis that the y_obs has been corrutped by some (known) noise delta_y. So we sample possible "true" values y_star to be centered around y_obs with the appropriate amount of noise. Now we use these values y_star as the dependent variable in our regression analysis.
A previous version of this answer required a kludge on the back-end, but I realized that we can just enforce the hard constraints on $y^*$ by means of a transformation. This model won't work in its current state, because we need to define what something is in the model statement, and I'm not quite sure how to go about it. Perhaps a smarter person could point out what, precisely, we meed. I'm not entirely sure there's sufficient information in the question for me to answer this. I leave this answer here as a partial solution which others might build on.
Stated another way, the specification of the errors in $y^*$ will permit one to replace this line
 y_star ~ normal(y_obs, something);

with something appropriate to OP's needs, and the model will be complete.
If you wish, you may place appropriate priors on beta as well.
A: 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
A: FIRST POST:
R's default variance-covariance of the fit takes account of weighting. If you need to convince yourself, run the unweighted model and see that the results of summary, vcov, and predict with se.fit=TRUE are all different.
The weighted least squares estimator is given by:
$$\hat{\beta}_{WLS} = \left( X^T W^{-1} X \right) ^{-1} X^T W^{-1} Y$$
and, when the weighting is properly specified,
$$\mbox{vcov} \left( \hat{\beta}_{WLS} \right)  = \hat{\sigma^2} \left( X^T W^{-1} X \right) ^{-1} $$
which is the efficient estimator by the Gauss-Markov theorem.
EDIT:
I see that you have constructed an interesting variance structure for $y$ that is nonnormal. Now, the important bit is that, regardless of the variance structure of $y$, the usual least squares estimate is still consistent for the right slope, although the error estimate may be conservative or anticonservative. Least squares happens to be the maximum likelihood estimate when the residuals are normal. Weighted estimates give different variance estimates. Due to the nature of the nonnormal residuals, we can't be certain of whether the weighted version is more efficient or not regardless of whether the weighting is correct.
The question remains: can you use ML for the known nonnormal variance structure? The answer is yes. For an unweighted estimate, the EM algorithm can do very well since the objective is to "squeeze" the LS line within a certain range of values. So in the max step, penalize the LARGEST observation since the ML estimate of a unif distribution is based on the maximum. Same logic applies to the weighted estimates... though (interestingly) I can't seem to get the hessian to come up with anything that makes a lick of sense. Tsk, the problem with boundary estimates! I would bootstrap these to get the real SEs!
delta_y <- 1
x <- 1:10
y <- 1:10 + runif(10, -delta_y, delta_y)
w <- 1:10
X <- cbind(1, x)

## unweighted:
logLik <- function(beta, y, X) {
  resid <- y - X %*% beta
  ## estep:
  dyest <- max(abs(resid))
  ## mstep:
  -10*dunif(dyest, -dyest, dyest, log=T)
}

unweighted <- nlm(logLik, p=c(0,1), X=X, y=y, hessian=T)

## weighted:
wLogLik <- function(beta, y, X, w) {
  resid <- y - X %*% beta
  ## estep:
  dyest <- max(abs(resid))
  ## mstep:
  -sum(w*dunif(dyest, -dyest, dyest, log=T))
}

weighted <- nlm(wLogLik, p=c(0,1), X=X, y=y, w=w, hessian=T)

## same betas:
weighted$estimate
    unweighted$estimate

## diff ses::
diag(solve(weighted$hessian))
    diag(solve(unweighted$hessian))

## bootstrap SEs:
val <- replicate(1000, {
  i <- sample(1:10, replace=T)
  y <- y[i]
  x <- x[i]
  w <- w[i]
  X <- cbind(1, x)
  weighted <- nlm(wLogLik, p=c(0,0), y=y, X=X, w=w)$estimate
      unweighted <- nlm(logLik, p=c(0,0), y=y, X=X)$estimate
  return(c(weighted, unweighted))
})

## ests
rowMeans(val)
apply(val, 1, sd) ## slightly different SEs

The problem here is that likelihoods are highly irregular, with tons of local maxima to capture you. They depend where you "start" and the variance estimates are horribly miscalibrated.
