Error propagation methods for orbit parameters 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    &amp; x      &amp;   \Delta x    &amp; y      &amp; \Delta y \\ \hline
1992.226 &amp; 0.104 &amp; 0.003 &amp; -0.166 &amp; 0.004 \\
1994.321 &amp; 0.097  &amp; 0.003 &amp; -0.189 &amp; 0.004 \\
...     &amp; ..     &amp; ..    &amp; ..     &amp;       \\
2003.353 &amp; 0.077  &amp; 0.002 &amp; -0.030 &amp; 0.002 \\
2003.454 &amp; 0.081  &amp; 0.002 &amp; -0.036 &amp; 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,
&quot;analytical&quot; 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?
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:

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}$$
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.
Determine $\Delta a, \ldots, \Delta f$ by partial differentiation.
Plugin in $\Delta x$ and $\Delta y$ as also given in your table.