Definitions
I have a set of data. Let's assume the true underlying model to be $\eta_i=\beta_0+\beta_1\zeta_i$. Where $\beta_i$ are the true coefficients of the model, $\eta_i$ are the true independent variables and $\zeta_i$ are the dependent variables. These cannot be observed.
The measurements (the data I have) are assumed to both have random errors ($x_i = \zeta_i + \delta_i$ and $y_i = \eta_i + \gamma_i$), which we assume to be normally distributed: $\delta_i \sim N(0,\sigma_{x_i}^2)$ and $\gamma_i \sim N(0,\sigma_{y_i}^2)$.
Through values from datasheets and error propagation, I have estimates of this distrubution: $s_{x_i} \approx \sigma_{x_i}$ and $s_{y_i} \approx \sigma_{y_i}$.
Problem
Because my measurements have errors on both axis, I use Total Least Squares (TLS) to fit a linear model to my data: $\hat y_i = b_0 + b_1 x_0$.
Confusingly, I found people calling this also Orthogonal Distance Regression (ODR) or Bivariate Least Squares (BLS). I use the
scipy.odr(which is using ODRPACK) implementation to be exact.
Now, let's assume I have a new independent data point which I assume follows the same model ($\eta_0=\beta_0+\beta_1\zeta_0$). Again, this true new data point cannot be observed, so I actually have $x_0$. Again, through looking at data sheets and propagating the error, I estimated $s_{x_0} \approx \sigma_{x_0}$.
I use the fitted model to predict $\hat y_0 = b_0 + b_1 x_0$. Now I want to calculate the variance of the prediction $s_{\hat y_0}^2$. How do I do that?
Additional info
I don't think this is relevant to the problem, but I include it to avoid an XY Problem.
I want to use $s_{y_0}$ to construct a prediction interval for $x_0$: $\hat y_0 \pm t_{\alpha,n-2}s_{\hat y_0}$.
Then I check if the measured $y_0$ falls outside of the prediction interval to detect an anomaly in the data. I need the probability, as one needs to know the probability that this new data point is an anomaly. I asked about this method in this question.
My journey
I already found a bunch of implementations, all of which give different results!
Approach 1
This question leads to these lecture notes and this website. The website references Wolberg, J., Data Analysis Using the Method of Least Squares, 2006, Springer for their approach.
I don't fully understand how they end up at their solution, I think they propagate the error of the prediction $\hat y_0 = b_0 + b_1 x_0$ using the Taylor expansion method, but I am unsure. These sources don't use different $s_{y_i}$. They assume one $s_y$ for all data points, including the data point to be predicted. The code calculates: $$ s_{\hat y_0}^2 = s_y^2 + s_{b_0}^2 + x_0^2 s_{b_1}^2 + 2x_0\text{cov}(b_0, b_1) $$ If I read Wolberg, Section 2.6 correctly, this method does not take into account the errors in $x$, but I am not sure.
Approach 2
I also came across this question, which found this paper.
Here, they somehow expand the implementation of the normal Least Squares method to arrive at: $$ s_{\hat y_0}^2 = [1 + \pmb X_0^T(\pmb X^T \pmb W^{-1} \pmb X)^{-1} \pmb X_0 + s_{x_0}^2b_1^2]s^2 $$ where $\pmb W$ is calculated from the weights for the TLS, and $s$ is the experimental error.
Approach 3
The same paper from Approach 2 claims that the same value can be also derived by another method, and references Meloun M, Militky J, Forina M. Chemometrics for Analytical Chemistry. Volume 1: PC-aided Statistical Data Analysis. Ellis Horwood: Chichester, 1992. To my knowledge they also propagate the error of the prediction using the Taylor expansion method. But they arrive at: $$ s_{\hat y_0}^2 = s^2 + s_{b_0}^2 + x_0^2s_{b_1}^2 + b_1^2s_{x_0}^2+ 2x_0\text{cov}(b_0, b_1) $$ Interestingly, this approach matches the approach 2, if I set $s=1$.
Side notes
The paper from Approach 2 and 3 actually calculate the variance of the response mean value (which in my eyes can be used for the confidence interval) with $s_{\hat y_0}^2 = [\pmb X_0^T(\pmb X^T \pmb W^{-1} \pmb X)^{-1} \pmb X_0 + s_{x_0}^2b_1^2]s^2$ or $s_{\hat y_0}^2 = s_{b_0}^2 + x_0^2s_{b_1}^2 + b_1^2s_{x_0}^2+ 2x_0\text{cov}(b_0, b_1)$. They then calculate the variance of the response variable (used for prediction interval) as the mean of $q$ observations with $s_{\hat y_0}^2 = [\frac1q+\pmb X_0^T(\pmb X^T \pmb W^{-1} \pmb X)^{-1} \pmb X_0 + s_{x_0}^2b_1^2]s^2$. As I only have one observation, $q=1$. And as approach 2 should match approach 3, I used the same logic for approach 3 to get the prediction interval.
Additionally, when researching what could have gone wrong, I looked at Melound et. al., Equation (1.40). And if I am not mistaken, the second term has a second partial derivative, which might be neglected in the paper? But I am very sure, that all of these sources have much more knowledge in statistics, so I am pretty sure I simply don't understand what's going on.
