MathOverflow is a question and answer site for professional mathematicians. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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?


share|cite|improve this question
up vote 1 down vote accepted

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

share|cite|improve this answer
@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? – Alex Aug 15 '11 at 16:54
Sorry, silly typo. Should work just fine for any set of basis functions. – Robby McKilliam Aug 15 '11 at 18:11

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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