this is my first post here. I tried my best to follow all the guidelines and to state the problem as clear as possible but should something be sloppy or unclear, I will be happy to update the question with further information.
Problem:
I'm trying to solve a 2D localization of a robot based on uncertain measurements.
Given a set of $n$ uncertain measurements of the robot's pose $\{\textbf{m}_1,\dots, \textbf{m}_n\}$, $\textbf{m}_i = \begin{bmatrix}x_i\\y_i\end{bmatrix} \in \mathbb{R}^2$, each associated with a covariance matrix $\mathbf{\Sigma}_i \in \mathbb{R}^{2\times2}$, my goal is to combine/fuse these measurements into one final estimate of robot position $\textbf{p} = \begin{bmatrix}x\\y \end{bmatrix}$ and an associated covariance matrix $\mathbf{\Sigma}$.
Example
Hopefully, the following figure will clarify what I mean by combine/fuse. It visualizes four measurements each having a mean $\textbf{m}_i$ (cross inside the ellipse) and a covariance matrix $\mathbf{\Sigma}_i$ visualized as the covariance ellipse.
The four measurements with covariance ellipses (red, green, brown and magenta) all depict an estimate of the same point $\textbf{p}$ only from different sources and I would like to create an estimate (mean and covariance) based on these measurements (something like the light blue ellipse).
What I'm doing so far
So far I've been trying to combine the measurements iteratively pair-wise using these Kalman-ish update equations $$ \mathbf{K} = \mathbf{\Sigma}_{i-1}\left(\mathbf{\Sigma}_{i-1}+\mathbf{\Sigma}_i\right)^{-1},\\ \textbf{p}= \textbf{m}_{i-1} + \textbf{K}\left(\textbf{m}_i - \textbf{m}_{i-1}\right),\\ \mathbf{\Sigma} = \mathbf{\Sigma}_{i-1} - \mathbf{K\Sigma}_{i-1}, $$ where $\textbf{m}_{i-1,i}$ are a pair of measurements, $\mathbf{\Sigma}_{i-1,i}$ are the associated covariance matrices and $\textbf{p}$ and $\mathbf{\Sigma}$ are the updated final position based on these two measurements. This yields some reasonable results, however I fear that the result depends on the order of combined measurements.
Question
I guess this would lead to performing something like a weighted average of the measurements, as per the answer here but I'm not entirely sure if this is the right approach and if it can be directly extended to the multivariate case. What is the correct way (mathematical operation) to go about it and merge the Gaussians?
If you haven't made an algebra mistake, the "Kalmanish" approach should work. This is actually a Kalman filter for an immobile object.
The way to think about these kind of problems is by using "Maximum Likelihood" estimation (MLE). You assume that there is a "true" position $p_{true}$ and work backwards: the probability to observe $m_i$ given $p_{true}$ is: $$P(m_i|p_{true})=G(m_i-p_{true},\Sigma_i) $$ where $G(x,\Sigma)$ is the Gaussian with covariance $\Sigma$. The probability to observe all of the $m_k$, assuming independence is the product $$P(m_1,,, m_n|p_{true})=\prod_i G(m_i-p_{true}, \Sigma_i)$$ Of course we do not know $p_{true}$, but MLE asks the question: what value of $p$ can we substitute in $p_{true}$ to maximize this product? As a sketch, we maximize the log of the product (same maximum argument): $p_{max likelihood}=\underset{p}{\arg\!\min} \sum_i (m_i-p)^T\Sigma_i (m_i-p) $ (selected the $\min$ to get rid of the minus in the exponent). If you do it recursively you will recover the Kalman formula, alternatively, you can get a closed formula for the position. As a bonus you can get the estimated covariance of the error of $p_{max likelihood}$ (exercises for the reader :)).