I have two datasets, A and B, of weighted (x,y) pairs. I computed the best fit lines, L_A and L_B, respectively, of these datasets, and then computed the intersection of these two lines, (x*,y*).
Now, A and B are drawn from random distributions (of the form y=mx+b+e, where e is a normally distributed error term). I'd like to compute a confidence interval around x*. How can I do this?
A naive approach I tried involved computing confidence intervals for the slope/intercept estimates of L_A and L_B, generating approximate distributions for those 4 parameters, and repeatedly randomly sampling from these distributions. But a given line's slope and intercept estimates are not independent of each other, and I wasn't sure how to salvage this approach to account for that.
Best Answer
One straightforward way is to obtain the maximum likelihood estimator of $x^{*}$ directly. Using the first subscript to designate the line ($1$ or $2$), the model is
$$y_{ij} = m_i(x_{ij} - x^{*}) + y^{*} + \varepsilon_{ij},$$
$1 \le j \le n_i$, $\varepsilon_{ij} \sim \text{Normal}(0, \sigma_i / \omega_{ij})$ and independent. The $n_i$ count the data for each line. The parameters are the point of intersection $(x^{*}, y^{*})$, the slopes $m_1$ and $m_2$, and (as nuisance parameters) the scale factors $\sigma_1$ and $\sigma_2$. The $\omega_{ij}$ are specified weights (not parameters). The standard ML machinery will provide a confidence interval for any of the parameters, including $x^{*}$.
This illustration show two linear fits from $10$ independent samples of two lines, each passing through $(7,10)$, with slopes $1$ and $1/3$. The colored bands are $95\%$ confidence bands around the fits (using least squares). The estimated point of intersection, shown as a large black dot, occurs at $(6.3, 9.4)$. A $95\%$ confidence interval for its x-coordinate is portrayed as a dashed black segment: it extends from $4.4$ to $8.3$.
Some details
(added in response to comments)
This model is identical to performing two separate weighted least squares regressions. It merely combines their 3+3 = 6 parameters (intercept, slope, and scale) in a way that isolates the x-coordinate of the point of intersection. Therefore the parameter estimates will be the same. The point is that in fitting the combined model, any maximum likelihood procedure will report standard errors (equivalently, confidence intervals) for the parameters and this is how we solve the original problem.
I would also like to point out that the there is some freedom in choosing the weights (which can be helpful in persuading an optimization routine to behave well). The foregoing analysis shows it suffices to examine the case of fitting a single line. Writing $\omega_i$ for the weights the negative log likelihood is equal to
$$-\log(\Lambda) = \frac{1}{2}\sum_i{\log(2 \pi (\sigma/\omega_i)^2) + \frac{(m(x_i-x^*) + y^* - y_i)^2}{(\sigma/\omega_i)^2}}.$$
Removing additive constants (which don't affect the ML procedure) this simplifies somewhat to
$$\sum_i{\log(\sigma) + \frac{u_i^2\omega_i^2}{2\sigma^2}}$$
where $u_i^2 = (m(x_i-x^*) + y^* - y_i)^2$. To estimate $\sigma$ we equate the derivative with $0$:
$$0 = \frac{\partial(-\log(\Lambda))}{\partial \sigma} = \sum_i^n{\frac{1}{\sigma} - \frac{u_i^2\omega_i^2}{\sigma^3}},$$
implying
$$(\hat{\sigma})^2 = \frac{1}{n}\sum_i^n{u_i^2 \omega_i^2}.$$
This shows that the estimate of $\sigma^2$ is directly proportional to the scale of the weights. This is clear: if we multiply all weights by a positive value $\lambda$, then we get exactly the same model by dividing $\sigma$ by $\lambda$. Consequently, we are free to normalize the weights as we wish. A good choice is to make $\sum{\omega_i^2}=n$, because this is what we would get for an unweighted model. In particular, note that the actual value of $\sigma$ is meaningless except in comparison to the $L^2$ norm of the weights. Accordingly, the standard error of $\hat{\sigma}$ is also meaningless except in comparison to the weights.