Linear regression on 3D position measurements

Question

I have 3D position measurements $x_i,y_i,z_i$ and their corresponding timestamps $t_i$ in a buffer. The time intervals are not equal between all timestamps. I would like to carry out linear regression, to obtain an estimate of the velocity (= slope of the linear fit) "during" the buffer. My idea is to perform linear regression with least squares.

I have multiple ideas but cannot connect the dots. Main question is under point 2, since I have some ideas but don't know if I am wrong or if I am right, how to exactly do it.

1. I would prefer a solution in vector and matrix notation, that solves the slopes/velocities in all 3 dimensions at the same time, if possible.

2. I know some linear regression with least squares from machine learning courses I took and was wondering how to apply that to this scenario. In linear regression we have matrix $\textbf{X}$, which contains our observations $\vec{x_1},...,\vec{x_N}$, we have $\vec{y}$, containing the target we want to predict, we have a weight vector $\vec{w}$ (or should I use a weight matrix $\textbf{W}$?)) and we have a bias $\vec{b}$ if we don't absorb it into the weights. Then the equation looks like this:
   $$\vec{y} = \vec{w}^T \textbf{X}+\vec{b}$$ But I have position and time measurements and I don't know which to treat as observations and which as targets, or if all of the four should be treated as measurements. Also since I have a time series/multiple measurements, how do I put my position and time measurements into matrices and vectors and adjust this equation?
   I also noticed that this resembles the motion equation: $$ \vec{s} = \begin{pmatrix} x \\ y \\z \end{pmatrix} = \vec{v} \cdot t + \vec{v_0}$$
   I guess if you consider the time series here, the time would become a vector and the velocities would become a Matrix, but so would the positions. Does this resemblance between the two formulae have some meaning / is it related to the solution of the linear regression or is it a coincidence and maybe even wrong?
   To summarize this part: How do I model the linear regression problem and how do I solve it with least squares? Is it somehow related to the motion equation?

3. My intuition tells me that I could simply obtain the linear fit by calculating the position differences $x_i - x_{i-1}$ and time intervals $t_i-t_{i-1}$, obtain the velocities between two measurements by dividing and then taking the average of these velocities between intervals. Then repeat for $y$ and $z$. Would this be equivalent to a specific solution method of linear regression?

Answer 1

You can probably start assuming multi-linear normal distribution for your measured position, and check whether resulting residuals agree with this assumption.

So let the observed position be $\mathbf{r}^{obs}_i=\left(x^{obs}_i,y^{obs}_i,\,z^{obs}_i\right)^T$, at time $t_i$ and the true position be:

$$ \mathbf{r}=\mathbf{v}\cdot\left(t-t_0\right)+\mathbf{r}_0 $$

With:

$$ \mathbf{r}^{obs}\sim N\left(\mathbf{r},\,\boldsymbol{\Sigma}\right) $$

Where constant covariance matrix $\mathbf{\Sigma}$.

We can approximate derivative with:

$$ \frac{d\mathbf{r}^{obs}}{dt}\Bigg|_{t=t_i}\approx \left[\frac{\mathbf{r}_{i+1}^{obs}-\mathbf{r}_{i}^{obs}}{t_{i+1}-t_{i}}+\frac{\mathbf{r}_{i}^{obs}-\mathbf{r}_{i-1}^{obs}}{t_{i}-t_{i-1}}\right]/2 $$

Which can be precise up to order $O\left(\delta t^2\right)$ if your time-steps are all equal, otherwise it will only be up to first order. I expect it should be possible to show that quantity on RHS will be multilinear normal, with equivalent mean, but the variance may now end up depending on time as well as $\boldsymbol{\Sigma}$.

The expectation value of the derivative will be the derivative of the expectation value (since expectation is linear):

$$ E\left[\frac{d\mathbf{r}^{obs}}{dt}\Bigg|_{t=t_i}\right]=\frac{d\mathbf{v}_i}{dt}\cdot\left(t_i-t_0\right)+\mathbf{v}_i $$

Components of the above quantity don't mix, so if the variance of $\frac{d\mathbf{r}^{obs}}{dt}$ did not depend on time too much you could run linear regression on individual components of $\frac{d\mathbf{r}^{obs}}{dt}\Bigg|_{t=t_i}$ the intercept of your result will be the components of velocity and the gradient will be the components of acceleration. With some extra work, and there you will have to go multilinear, you will also be able to work out the covariance matrix $\boldsymbol{\Sigma}$ and how it translates into the uncertainty in estimating your velocity.

If the time dependence of your covariance matrix is not small enough, you may need to go with Max Likelihood estimation, or some other method that can accommodate heteroscedasticity. I would probably try that first.