How to Derive Prediction Intervals for Orthogonal Distance Regression using `scipy.odr`? Questions

*

*How can I derive prediction intervals for predictions based on new observations from the output of scipy.odr?

*Is it also possible (or necessary) to take into account uncertainties in the new observations?

Background
I would like to perform a linear regression between two sets of variables, both of which have uncertainties associated with them, and also be able to derive prediction intervals as part of the output.
I have identified Orthogonal Distance Regression (ODR) as a possible method by which this could be done, and the scipy.odr library as an implementation of this method.
What I have been unable to find is how to use the output of this library to derive prediction intervals. I understand from this page that I need to add - in quadrature - the estimate of the standard deviation of the predicted value and the estimate of the residual standard deviation obtained when fitting the model to the data, but I am unsure how to derive these quantities from the output of scipy.odr. I am also almost completely unsure how to take into account any uncertainties in new observations, as surely these will affect the resulting prediction intervals.
 A: [I am working on researching a complete answer, this is as far as I've got so far, I will edit as I work out more of the solution.]
Further research has turned up this paper which provides at least part of the answer.
Briefly, ordinary least squares (OLS) does not consider uncertainties on either axis, and weighted least squares (WLS) only takes into account uncertainties in the predictor variable / dependent variable / on the y-axis. The above publication advocates the use of bivariate least squares (BLS) over orthogonal distance regression (ODR), but goes on to derive formulae for prediction intervals which appear to be compatible with ODR.
In this approach, each 'data point' $\left(x_{i}, y_{i}\right)$ is considered to be the result of multiple experiments or observations under ostensibly the same conditions, thus allowing the calculation of the variance for both the predictor and response variable pair.
The predictor ($x_{i}$) and response ($y_{i}$) variables are related by the following equation
$y_{i}=b_{0}+b_{1} x_{i}+e_{i}$
where $b_{0}$ and $b_{1}$ are the estimates of the intercept and slope of the true linear model, and $e_{i}$ is the $i$th residual error. The variance of $e_{i}$ is referred to as the weighting factor and is denoted as $w_{i}$ or $s_{e_{i}}^{2}$:
$w_{i}=s_{e_{i}}^{2}=s_{y_{i}}^{2}+b_{1}^{2} s_{x_{i}}^{2}-2 b_{1} \operatorname{cov}\left(x_{i}, y_{i}\right)$
where $s_{x_{i}}^{2}$ and $s_{y_{i}}^{2}$ are the experimental variances of point $i$, and $\operatorname{cov}\left(x_{i}, y_{i}\right)$ is the covariance between the measurements for each $\left(x_{i}, y_{i}\right)$ data pair.
The variance of the response variable $y_{0}$, being the mean of the $q$ observations performed at $x_{0}$, is given by:
$s_{y_{0}}^{2}=\left(\frac{1}{q}+ X _{0}^{ T }\left( X ^{ T } W ^{-1} X \right)^{-1} X _{0}+s_{x_{0}}^{2} b_{1}^{2}\right) s^{2}$
where $s^{2}=\frac{\sum_{i=1}^{n}\frac{\left(y_{i}-\hat{y}_{i}\right)^{2}}{w_{i}}}{n-2}$ is the estimate of the true experimental variance, $X_{0}$ is a two element column vector $\left|\begin{array}{l}
1 \\
x_{0}
\end{array}\right|$, and $X$ is an $n \times 2$ matrix for which the first column is a column of ones and the second column is formed by the $n$ values of $x$ corresponding to the experimental points; $W$ is an $n \times n$ diagonal matrix whose $i$th diagonal element is the weighting factor $w_{i}$, defined above.
The variance of the predictor mean value at a given observation $y_{0}$ is
$s_{x_{0}}^{2}=\left( Y _{0}^{ T }\left( Y ^{ T } W ^{\prime -1} Y \right)^{-1} Y _{0}+s_{y_{0}}^{2} \frac{1}{b_{1}^{2}}\right) s^{\prime 2}$
where $Y _{0}$ is $\left|\begin{array}{l}
1 \\
y_{0}
\end{array}\right|$, $Y$ is an $n \times 2$ matrix for which the first column is a column of ones and the second column is the values of $y$, $W^{\prime}$ is an $n \times n$ diagonal matrix whose $i$th diagonal element is the weighting factor  $w_{i}^{\prime}=s_{x_{i}}^{2}+\frac{1}{b_{1}^{2}} s_{y_{i}}^{2}-2 \frac{1}{b_{1}} \operatorname{cov}\left(x_{i}, y_{i}\right)$, and $s^{\prime 2}$ the experimental error associated with the predictor variables, given by
$s^{\prime 2}=\frac{\sum_{i=1}^{n} \frac{\left(y_{i}-\hat{y}_{i}\right)^{2}}{w_{i}^{\prime}}}{n-2}$
The variance of the prediction of the predictor variable of a future sample at $y_{0}$, the mean of $q$ observations, is
$s_{x_{0}}^{2}=\left(\frac{1}{q}+ Y _{0}^{ T }\left( Y ^{ T } W ^{-1} Y \right)^{-1} Y _{0}+s_{y_{0}}^{2} \frac{1}{b_{1}^{2}}\right) s^{\prime 2}$
The prediction intervals of the response and predictor variables are then given by
$y_{0} \pm t_{\alpha, n-2} s_{y_{0}}$
$x_{0} \pm t_{\alpha, n-2} s_{x_{0}}$
where $t_{\alpha, n-2}$ is the $t$-value for the required significance level $\alpha$ and $n-2$ degrees of freedom.
TODO: Interface to scipy.odr output.
