Let $\mathcal{Y} = (\mathbf{y}_1, \dots, \mathbf{y}_N)$ be data observed, such that each $\mathbf{y}_i \in \mathbb{R}^2$. Now conditional on unobserved cluster centres (means) $\mathcal{X} = (\mathbf{x}_1, \dots, \mathbf{x}_K)$, where $K >>0, K \in [1,\infty)$ is also unknown, we have
$$L(\mathcal{Y}|\mathcal{X}, \mathcal{C}, \Sigma,K) = \prod_{i =1}^N \prod_{j=1}^K [\mathcal{N}_2(\mathbf{y}_i;\mathbf{x}_j,\Sigma)]^{\mathcal{C}_{ij}},$$ where $\Sigma$ is a known covariance matrix and $\mathcal{C}_{ij} = 1$ if observation $i$ belongs to cluster $j$ and is zero otherwise.
Now, denote the number of observations per cluster as $S_1, \dots S_K$. We have that $S_j = \sum_{i =1}^{N} \mathcal{C}_{ij}$ for each $j = 1, \dots, K.$ The distribution of each $S_j$ follows $S_j \sim p(S)$ independently and identically, and while it is not multinomial, it is known.
I would like to infer $K$ and $\mathcal{X}$ given the above set-up. Since $K$ in this application is typically quite large (and can be infinite), I was hoping to use some theory in Bayesian nonparametrics, specifically from Dirichlet mixture models. However, since I do not have a typical Gaussian Mixture Model, I am unsure how to set this up.
Any help or pointers to specific papers which outline the theory that I need would be very kindly appreciated. Thanks!