## Appendix A Derivation of Approximate Inference Algorithm  
The derivation of the algorithm for Bayesian estimation of a BMM model (Algorithm 2) from Equations (5)-(10) is outlined in this appendix. We combine Equations (5)-(10) using Bayesian rules to obtain the joint density function of the observation \(X\) and all the i.i.d. latent variables \(\mathcal{Z}=\{\mathbf{A},\mathbf{B},\mathbf{P},\mathbf{Z}\}\). We assume that \(X\) is conditionally independent of \(P\), given \(\mathbf{Z}\). The joint density function is given by:  
\[f(\mathbf{X},\mathcal{Z})=f(\mathbf{X},\mathbf{A},\mathbf{B}, \mathbf{P},\mathbf{Z})\] \[= f(\mathbf{X}|\mathbf{A},\mathbf{B},\mathbf{Z})f(\mathbf{Z}| \mathbf{P})f(\mathbf{P})f(\mathbf{A})f(\mathbf{B})\] \[= \prod_{n=1}^{N}\prod_{m=1}^{M}\left[p_{m}\mathbf{Beta}(\mathbf{x }_{n}|\mathbf{a}_{m},\mathbf{b}_{m})\right]^{z_{nm}}\cdot C(\mathbf{c})\prod_{ m=1}^{M}p_{m}^{c_{m}-1} \tag{11}\] \[\cdot \prod_{l=1}^{L}\prod_{m=1}^{M}\left[\frac{\alpha_{lm}^{\mu_{lm}} }{\Gamma(\mu_{lm})}a_{lm}^{\mu_{lm}-1}e^{-\alpha_{lm}a_{lm}}\right]\cdot\left[ \frac{\beta_{lm}^{\mu_{lm}}}{\Gamma(\
u_{lm})}b_{lm}^{\mu_{lm}-1}e^{-\beta_{ lm}b_{lm}}\right]\]  
and the logarithm of Equation (11) is  
\[\mathcal{L}(\mathbf{X},\mathcal{Z})=\ln f(\mathbf{X},\mathbf{Z},\mathbf{A},\mathbf{B},\mathbf{P})\] \[= \text{con.}+\sum_{n=1}^{N}\sum_{m=1}^{M}z_{nm}\left\{\ln p_{m}+ \sum_{l=1}^{L}\ln\frac{\Gamma(a_{lm}+b_{lm})}{\Gamma(a_{lm})\Gamma(b_{lm})}\right.\] \[\left.+\sum_{l=1}^{L}\left[(a_{lm}-1)\ln x_{ln}+(b_{lm}-1)\ln(1- x_{ln})\right]\right\} \tag{12}\] \[+ \sum_{l=1}^{L}\sum_{m=1}^{M}\left[(\mu_{lm}-1)\ln a_{lm}-\alpha_ {lm}a_{lm}\right]\] \[+ \sum_{l=1}^{L}\sum_{m=1}^{M}\left[(\
u_{lm}-1)\ln b_{lm}-\beta_{ lm}b_{lm}\right]+\sum_{m=1}^{M}(c_{m}-1)\ln p_{m}\]  
The i.i.d. latent variables we have now are \(\mathbf{A}\), \(\mathbf{B}\), \(\mathbf{P}\) and \(\mathbf{Z}\) with the hyperparameters \(\boldsymbol{\alpha}\), \(\boldsymbol{\beta}\), \(\boldsymbol{\mu}\), \(\boldsymbol{\
u}\) and \(\mathbf{c}\). The updating scheme of the VI can be used to estimate these hyperparameters of the latent variables. Following the principles of the VI framework [32; 33; 34; 15; 35], we have:  
\[c_{m}^{\star}=c_{m_{0}}+\sum_{n=1}^{N}\mathbf{E}\left[z_{nm}\right]=c_{m_{0}}+ \sum_{n=1}^{N}r_{nm} \tag{13}\]  
where \(c_{m_{0}}\) denotes the initial value of \(c_{m}\). We can also obtain the closed-form updating scheme for the hyperparameters \(\alpha_{lm}\) and \(\mu_{lm}\):  
\[= f(\mathbf{X}|\mathbf{A},\mathbf{B},\mathbf{Z})f(\mathbf{Z}| \mathbf{P})f(\mathbf{P})f(\mathbf{A})f(\mathbf{B})\]0 \[= f(\mathbf{X}|\mathbf{A},\mathbf{B},\mathbf{Z})f(\mathbf{Z}| \mathbf{P})f(\mathbf{P})f(\mathbf{A})f(\mathbf{B})\]1 \[= f(\mathbf{X}|\mathbf{A},\mathbf{B},\mathbf{Z})f(\mathbf{Z}| \mathbf{P})f(\mathbf{P})f(\mathbf{A})f(\mathbf{B})\]2  
For the same reasons, the update schemes for the hyperparameters \(\beta_{lm}\) and \(\
u_{lm}\) are:  
\[= f(\mathbf{X}|\mathbf{A},\mathbf{B},\mathbf{Z})f(\mathbf{Z}| \mathbf{P})f(\mathbf{P})f(\mathbf{A})f(\mathbf{B})\]3 \[= f(\mathbf{X}|\mathbf{A},\mathbf{B},\mathbf{Z})f(\mathbf{Z}| \mathbf{P})f(\mathbf{P})f(\mathbf{A})f(\mathbf{B})\]4 \[= f(\mathbf{X}|\mathbf{A},\mathbf{B},\mathbf{Z})f(\mathbf{Z}| \mathbf{P})f(\mathbf{P})f(\mathbf{A})f(\mathbf{B})\]5Further details of the derivations are given in [26, 27]. The update Equations (A3)-(A7) are calculated through the following expectations:  
\[ \begin{split}&\overline{a}=\frac{\mu}{\alpha},\;\overline{b}=\frac{ \
u}{\beta},\;\mathbf{E}\left[z_{nm}\right]=r_{nm},\;\mathbf{E}\left[\ln p_{m} \right]=\psi(c_{m})-\psi(\widehat{c}),\\ &\mathbf{E}_{a}[\ln a]=\psi(\mu)-\ln\alpha,\;\mathbf{E}_{b}[\ln b ]=\psi(\
u)-\ln\beta,\\ &\mathbf{E}_{a}\left[(\ln a-\ln\overline{a})^{2}\right]=[\psi( \mu)-(\ln\mu)]^{2}+\psi^{{}^{\prime}}(\mu)\\ &\mathbf{E}_{b}\left[(\ln b-\ln\overline{b})^{2}\right]=[\psi( \
u)-(\ln\
u)]^{2}+\psi^{{}^{\prime}}(\
u)\end{split}  \tag{10}\]