Point-wise error estimate in polynomial regression In our application we wish to estimate the actual path of objects. We have a set of samples of object locations $(t_i, x_i, y_i, P_i)$ where $t_i$ is the sample time, $(x_i, y_i)$ is the 2D location, and $P_i$ is the error covariance matrix for measurement $i$.
We estimate the path using polynomial regression as a function of $t$. That is, we have two polynomials $x(t) = a_x + b_x t + c_x t^2 + d_x t^3$ and $y(t) = a_y + b_y t + c_y t^2 + d_y t^3$. We find $\{a, b, c, d\}$ using a log-likelihood estimator assuming gaussian errors distribution. That is, we minimize:
$$\Sigma (x(t_i) - x_i, y(t_i) - y_i) P_i^{-1} (x(t_i) - x_i, y(t_i) - y_i)^T$$
I would like to be able to estimate the approximation error for each $t$. After finding our point $(x(t), y(t))$ using the above polynomials, we also want the error covariance matrix $P(t)$. Does anybody know of a method of doing so? (it might be related to prediction/error bands for simple linear regression)
Note that when $P_i$ are diagonal matrices we get the simple linear regression that has known coefficient error estimates that I can use. Does anybody know how can this be done for the non-diagonal case?
Alex.
 A: This is just weighted least squares and here is how I would approach it. To keep my notation simple I'll just have polynomials of order 1. It's trivial to extend the approach to polynomials of any order. Let $a_x$, $b_x$ and $a_y$, $b_y$ be the 'true' coefficients describing the path.  So the path is
$$ (x(t), y(t)) = ( a_x + b_x t, a_y + b_y t). $$
For convenience I'll pack the coefficients into a column vector
$$p_0 = \left[\begin{array}{cccc} a_x & b_x & a_y & b_y \end{array}\right]^\prime $$
where $'$ denotes transpose.
The path is observed at a number of points in time, $t_1, \dots,t_n$ say, and these observations as subject to noise with known covariance matrices $P_1,\dots,P_n$.  That is, the observations are pairs
$$(x_i, y_i) = (x(t_i), y(t_i)) + W_i$$ 
where $W_1,\dots,W_n$ are bivariate random variables $W_i = (X_i, Y_i)$ with individual covariances given by the $2\times 2$ matrices $P_1, \dots, P_n$.  I'll assume that the $W_1, \dots, W_n$ are independent.
We can write the observations in vector form as
$$f = T p_0 + w$$
where
$$ f = \left[ \begin{array}{c} x_1 \newline y_1 \newline x_2 \newline y_2 \newline \vdots \newline x_n \newline y_n \end{array}\right]^\prime \qquad
T = \left[ \begin{array}{cccc} 
1 & t_1 & 0 & 0 \newline
0 & 0 & 1 & t_1 \newline
1 & t_2  & 0 & 0 \newline 
0 & 0 & 1 & t_2 \newline
\vdots & \vdots & \vdots & \vdots \newline
1 & t_n  & 0 & 0 \newline
0 & 0 & 1 & t_n
 \end{array} \right] \qquad
w = \left[ \begin{array}{c} X_1 \newline Y_1 \newline X_2 \newline Y_2 \newline \vdots \newline X_n \newline Y_n \end{array}\right]^\prime.$$
Let $P$ be the $2n\times 2n$ covariance of $w$.  So $P$ is block diagonal with diagonals given by the $2\times 2$ matrices $P_1,\dots,P_n$.
You take a weighted least squares approach to estimation, that is, your estimators are given by the minimisers of the quadratic form
$$(f - Tp)^\prime P^{-1} (f - Tp). $$
The minimiser is given by
$$ \begin{array}{ll}
\hat{p} &= (T^\prime D T)^{-1} T^\prime D f \newline 
&= M (Tp_0 + w) \newline 
&= p_0 + Mw \newline 
\end{array}$$
where $D = P^{-1}$ and $M = (T^\prime D T)^{-1} T^\prime D$.  So the error in your coefficients is given by $\hat{p} - p_0 = Mw$ and the covariance of the error is
$$ C = \operatorname{cov}(Mw) = M \operatorname{cov}(w) M^\prime = (T^\prime D T)^{-1}. $$
You want to know the covariance of the error at time $t$, that is you want the covariance of
$$\left[\begin{array}{c} \hat{a}_x + \hat{b}_x t \newline \hat{a}_y + \hat{b}_y t \end{array}\right] - \left[\begin{array}{c} a_x + b_x t \newline a_y + b_y t \end{array}\right] \qquad \text{where} \qquad \hat{p} = \left[\begin{array}{cccc} \hat{a}_x & \hat{b}_x & \hat{a}_y & \hat{b}_y \end{array}\right]^\prime. $$
This is given by the covariance of
$$ \left[\begin{array}{cccc} 1 & t & 0 & 0 \newline 
0 & 0 & 1 & t \end{array}\right] (\hat{p} - p_0) = K(t)(\hat{p} - p_0) = K M w $$
say.  This is $K(t) (T^\prime D T)^{-1} K(t)^\prime$.
