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.



Point-wise error estimate in polynomial regression
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$.


