Error propagation methods for orbit parameters

Question

Recently I've encountered an article in our local astronomical popularization magazine. It is about a well known task of Sgr A*-black hole's mass estimation. The article is written as a step by step guide for young researchers. For example, the orbital parameters of the S2 star should be found using a drawing on a sheet of paper.

That is why I believed there is no need for very precise calculations for this task.

However, in addition to the $ x $, $ y $ coordinates of the S2 star, I was also given the corresponding uncertainty values $ \Delta x $, $ \Delta y $ (19 coordinate measurements for the time period from 1992 to 2003):

\begin{array}{|c|c|c|c|} \hline time & x & \Delta x & y & \Delta y \\ \hline 1992.226 & 0.104 & 0.003 & -0.166 & 0.004 \\ 1994.321 & 0.097 & 0.003 & -0.189 & 0.004 \\ ... & .. & .. & .. & \\ 2003.353 & 0.077 & 0.002 & -0.030 & 0.002 \\ 2003.454 & 0.081 & 0.002 & -0.036 & 0.002\\ \hline \end{array}

Using direct least squares method for S2 orbital fitting by an ellipse I've found parameters $a,\ b,\ c,\ d,\ e,\ f$ of the conic equation: $$ax^2+bxy+cy^2+ dx+ey+f=0$$ I've also decided to simulate the error propagation in a Monte Carlo fashion to find the uncertainty in the estimates of $a,\ b,\ c,\ d,\ e,\ f$. And I succeeded in that.

But, here is the question: is there any other appropriate, "analytical" way (not a Monte Carlo method) to find the uncertainty in the estimates of $a,\ b,\ c,\ d,\ e,\ f$ which may be used in usual astronomers' practice?

Answer 1

ProfRob's comment is actually already the answer:

Unless you can linearise the model, there is no analytic error estimate.

I consider error estimation very important, so I would like to show a bit more details here: A way of estimating the error of a model is by Gaussian error propagation, see in particular the part of the Wikipedia entry on non-linear combinations. Let me briefly summarize the idea behind as I find the Wikipedia article not so easy to follow. I will also use different variable names so that they do not overlap with those in your question.

We start with a function $\varphi$ which depends on different variables $x_1, x_2, \ldots$, meaning $\varphi = \varphi(x_1, x_2, \ldots)$. We assume that we can linearize around a certain point $\tilde{\bf x} = (\tilde{x}_1, \tilde{x}_2, \ldots)$, which mainly means writing the function at that point as a Taylor series.

Then, we can determine the Gaussian (maximal) error in the vicinity of ${\bf \tilde{x}}$.

$$\left. \Delta \varphi \right|_{\bf \tilde{x}} = \left. \frac{\partial \varphi}{\partial x_1} \right|_{(\tilde{x}_1,\, \tilde{x}_2,\, \ldots )} \!\!\!\!\cdot \Delta x_1 + \left. \frac{\partial \varphi}{\partial x_2} \right|_{(\tilde{x}_1,\, \tilde{x}_2,\, \ldots )} \!\!\!\!\cdot \Delta x_2+\cdots$$

In this formula you calculate partial derivatives $\frac{\partial \varphi}{\partial x_i}$ with respect to each variable $x_i$ and plug in the values of the variables at the central point ${\bf \tilde{x}}$. The $\Delta x_i$ is the error estimate for each individual variable $x_i$.

How to we apply that to your problem where you have fitted an implicit equation? Since you excluded the Monte-Carlo approach (which is the usual method for this case), I suggest the following recipe:

1. Formulate explicit functions $a(x,y), b(x,y), \ldots f(x,y)$, e.g. $$ a = - \frac{bxy + cy^2+ dx+ ey + f}{x^2}$$
2. Linearize the functions around $x_0,y_0$. The $x_0$ is called $x$ in your table, and respectively $y_0$ is in the $y$-column.
3. Determine $\Delta a, \ldots, \Delta f$ by partial differentiation.
4. Plugin in $\Delta x$ and $\Delta y$ as also given in your table.