Point-wise error estimate in polynomial regression

Question

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.

Answer 1

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$.