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?



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

  • $\begingroup$ @Robby X′DX is T′DT, right? Because I can't find X defined anywhere before its first usage. Anyway, it seems that this is correct for any degree of polynomials (or even any set of basis functions), right? $\endgroup$ – Alex Shtof Aug 15 '11 at 16:54
  • $\begingroup$ Sorry, silly typo. Should work just fine for any set of basis functions. $\endgroup$ – Robby McKilliam Aug 15 '11 at 18:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.