I have measurements of a trajectory that is parameterized by time. The data consists of points with two spacial coordinates $(\tilde{x}_i, \tilde{y}_i)$ and a time stamp $(t_i)$. I'm using numpy.polyfit to fit a 2nd order polynomial to the data:

enter image description here

Since the data describe the motion of a physical object in 2D space, the parameters from the polynomial represent the velocity and acceleration of the object (and initial position).


My data points (of course) do not perfectly fit the model, thus there are residuals in $x$ and $y$ (i. e. the measured coordinates have some offset from where they should lie according to the fit at the corresponding time):

enter image description here

Because the derived velocity and acceleration are important quantities for me, I came up with the following system to test how reliable they are:

  • I assume that the curve from the fit truely represents the actual trajectory, and that my measured data points $(\tilde{x}_i, \tilde{y}_i)$ only deviate from the true locations $(x_i, y_i)$ by a gaussian offset:

$\begin{pmatrix} \tilde{x}_i\\ \tilde{y}_i \end{pmatrix} = \begin{pmatrix} x_i\\ y_i \end{pmatrix} + \begin{pmatrix} \delta_{x, i}\\ \delta_{y, i} \end{pmatrix}$

  • Based on this assumption, I create fake data sets $(\hat{x}_i, \hat{y}_i)$ by adding some random offset drawn from gaussian distributions that have the same standard deviations ($\sigma_x$, $\sigma_y$) as the offset distribution of the measured data:

$\begin{pmatrix} \hat{x}_i\\ \hat{y}_i \end{pmatrix} = \begin{pmatrix} x_i\\ y_i \end{pmatrix} + \begin{pmatrix} \hat{\delta}_{x, i}\\ \hat{\delta}_{y, i} \end{pmatrix}$

  • If I'm not mistaken, this should mean that the real dataset is basically indistinguishable from the genrated ones.

  • I then fit polynomials to the generated data sets and compare the velocity and acceleration vectors to those of the measured data. Note that in the plots below, the black outline indicates the vectors of the measured data, while the semi-transparent blue ones are from the generated data. The gray dashed lines indicate the mean values. (The vectors are warped due to non-square aspect ratios.):

enter image description here

For these velocities and accelerations I would then calculate the standard deviations and use them as error estimates for the measured values. Is that reasonable? Or am I commiting some kind of fallacy? I'm suspicious because the vectors from the measured data are awefully close to the mean values.

Also, I guess this counts as a Monte Carlo simulation. Is that correct?

Additional Note

The numpy.polyfit function can also provide diagnostic information like the covariance matrix and singular values of the scaled Vandermonde coefficient matrix, and I would guess that they can also tell me similar things, however I'm not exactly sure because I don't really know what they mean (see also link, link).

  • 1
    $\begingroup$ The parameters of a "3rd order polynomial" also include the rate of change of the acceleration. Are you truly fitting 3rd order polynomials, or are you just fitting 2nd order polynomials? As far as the substance of your question goes, it is natural to analyze the 2D distribution of the residuals. You can map them, for instance, and summarize them using bivariate statistics such as their first and second moments. A plot that connects 2D residuals according to their time sequence can be especially revealing. You can also use Complex numbers, as in stats.stackexchange.com/a/66268/919. $\endgroup$
    – whuber
    Jan 31, 2022 at 17:25
  • 1
    $\begingroup$ Thanks @whuber, you're correct, I meant 2nd order (though in reality it actually depends; I also fit 3rd order polynomials). But I guess that shouldn't change the outcome. And thanks for the other information. I'll check it out! $\endgroup$
    – mapf
    Jan 31, 2022 at 17:51

1 Answer 1


Given a starting position $x_0$, starting velocity $v_{x,0}$ and fixed acceleration $a_x$, you can describe the path of the x-coordinate in time as $$x(t) = x_0 + v_{x,0} \cdot t + \frac{1}{2} a_x t^2$$

Is that the polynomial that you are fitting? $x(t)$ as function of $t$ and $y(t)$ as function of $t$?

(The coordinates $x(t)$ as function of $y(t)$ does not generally follow a polynomial function for this situation)

In that case you can estimate the distribution/error of the parameter estimates in the typical way for ordinary least squares regression (using $(X^tX)^{-1}$ as described here on Wikipedia).

You don't need to use simulations for this (the simulations will just approximate the exact computation).

  • $\begingroup$ It might be that your time measurement does also have some error. If you assume that the error in $y$ and $x$ are independent then you could estimate the error in $t$ (since an error in $t$ will cause correlation in the residuals, and from that correlation you can estimate the error in $t$). $\endgroup$ Feb 2, 2022 at 20:20
  • $\begingroup$ Although, for these measurements I wonder whether you need to be very precise. Are you trying to estimate/predict whether a 10km size comet is gonna hit earth? $\endgroup$ Feb 2, 2022 at 20:21
  • $\begingroup$ Hi, yes, that's the kind of polynomial that I am fitting. In theory the time measurement also has an error but compared to the spacial errors it's totally irrelevant. But thanks for the tip! It's pretty funny how specific your question about the comet is! I'm actually doing something quite similar. I'm trying to estimate the source region of cometary debris by extrapolating particle trajectories back in time. But the dynamics are already very interesting on their own, so I'd like to have an idea about how reliable they are. $\endgroup$
    – mapf
    Feb 3, 2022 at 7:53
  • $\begingroup$ Can you please explain what you mean by the typical way for ordinary least squares regression though? I don't have a very solid math background, so I don't really know what's going on in that wiki article. I would bet that a lot of that information is already provided by the diagnostics from numpy.polyfit. I just don't really know what they mean either. $\endgroup$
    – mapf
    Feb 3, 2022 at 7:58
  • 1
    $\begingroup$ What your are doing is indeed a Monte Carlo simulation in order to estimate a distribution of the estimate derived from the measurement (which tells you something about the error). But that distribution can also be computed directly. The least squares regression is namely a linear estimator (it's a weighted sum of your measurements $x$ and $y$). If the measurements $x$ and $y$ are normal distributed then a sum of the $x$ and $y$ is as well. So you can compute the variance of the estimated parameters directly based on an estimates of the variance in the distribution of the $x$ and $y$. $\endgroup$ Feb 3, 2022 at 8:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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