of $k_1$ values. The EM algorithm is not guaranteed to find the global maximum of the likelihood landscape, but repeated runs with subsamples of the complete data set did not result in large variations, either in the number of components, or in the fitted parameters.

We denote the values of $\vec{\theta}_1^k$ that maximise the likelihood function with $k$ Gaussian components as $\vec{\hat\theta_1^k}$. ${\cal L}_{\hat\theta^k_1}$ is the maximum likelihood found by the EM algorithm for the family of models defined as a sum of $k$ multivariate normal distributions with full covariance matrices, and the data set defined by $n$ data points (sources) of dimension $d$. We then selected the likelihood model that maximises the Bayesian information criterion (BIC), given by

$${\cal M}= -2\ln({\cal L}_{\hat\theta_1^k})+k\frac{d(d+3)}{2}\ln(n).$$

The BIC represents a regularised figure-of-merit with a penalty term for complex models derived in the Bayesian framework for multivariate distributions [@Schwarz1978Estimating].

Principal curve models for the cluster members {#PC}
----------------------------------------------

In Sect. \[MG\] we modelled the likelihood term in Eq. \[bayes.th\] for the field sources as a sum of multivariate normal distributions. To avoid the overfitting of the data, we used a regularised maximum-likelihood criterion with a penalty for complex models (i.e., the BIC). While this 