self-consistency, we adopted an iterative procedure, based again on the EM algorithm, whereby we derived likelihood models in iteration $l$ from the membership lists calculated in iteration $l-1$ (maximisation step). Then, we used these likelihood models to infer membership probabilities (expectation step), and proceeded to the next iteration, $l+1$.

The first iteration proceeded as described in the previous sections, with the data set $\mathcal{D}$ defined by the set of sources observed in all the variables in representation space, and with the class labels ${\cal C}$ taken from [@2007ApJS..172..663S]. Let $p_0({\cal D}|{\cal C}_i)$ denote the initial likelihood models derived from the initial reference set. In subsequent iterations, the regularised maximum-likelihood model, together with Eq. \[bayes.th\], were used to infer posterior membership probabilities (expectation step). Since the two models by definition comprise all possibilities, we can compute the posterior probability in iteration $l$ as

$$p_l({\cal C}_i|{\cal D}) = \frac{p_{l-1}({\cal D}|{\cal C}_i)\cdot p_{l-1}({\cal
    C}_i)}{\sum_{i=1}^2p_{l-1}({\cal D}|{\cal C}_i)\cdot p_{l-1}({\cal C}_{i})}.$$

We took into account the measurement uncertainties in the likelihood term as follows:

$$\label{uncertainties}
p_{l-1}(\vec{m}_{\rm measured}|{\cal C}_i)= \int p_{l-1}(\vec{m}|{\cal C}_i)\cdot
p_{l-1}(\vec{m}_{measured}|\vec{m}) {\rm d}\vec{m}$$

where $\vec{m}_{measured}$ represents the measured values. We assumed independent Gaussian uncertainties in all measured variables (proper motions and apparent magnitudes) and thus,

