\section{Mixture Models}\label{sec:mixtures}
Mixture models are generative models focusing on how data could have created. Each point $x_n$ is supposed to be generated from one of the $K$ distributions, that can be easily sampled to generate new points.\\
Given the set of the datapoints $\mathbf{X}$, we define as $\Delta$ the set of all the parameters of the distributions. The model is a reverse-engineering process, aiming to learn $\Delta$, and thus the originating distribution associated to each point.\\
The objective is to maximize the \textit{mixture model likelihood}. Let the $k$-th distribution have probability density function:
$$ p(x_n |z_{nk}=1,\Delta_k)$$
The marginal likelihood associated to data $\mathbf{X}$ is $p(\mathbf{X}|\Delta)$.\\
We leverage the i.i.d. (independent and identically distributed) assumption over the data-points $x_n \in \mathbf{X}, \ n=1,\dots N$, therefore factorizing:
$$ p(X|\Delta) = \prod_{n=1}^N p(x_n|\Delta)$$
Then, we un-marginalize $k$:
$$ p(\mathbf{X}|\Delta) = \prod_{n=1}^N \sum_{k=1}^K p(x_n,z_{nk}=1|\Delta) = \prod_{n=1}^N \sum_{k=1}^K p(x_n|z_{nk}=1,\Delta_k) p(z_{nk}=1|\Delta)$$
Hence, since we want to find $\Delta$ maximizing it:
$$ \underset{\Delta}{\mathrm{argmax} \ }\prod_{n=1}^N \sum_{k=1}^K p(x_n|z_{nk}=1,\Delta_k) p(z_{nk}=1|\Delta)$$
We switch to log-probabilities to simplify the optimization:
$$ \underset{\Delta}{\mathrm{argmax} \ }\sum_{n=1}^N \log \sum_{k=1}^K p(x_n|z_{nk}=1,\Delta_k) p(z_{nk}=1|\Delta)$$
However, applying a logarithm of a sum could be tricky to be optimized, then we resort to the \textit{Jensen's Inequality}:
$$ \log \mathbb{E}_{p(x)} [f(x)] \geq \mathbb{E}_{p(x)}[\log f(x)]$$
So, given our likelihood $\mathcal{L} = \sum_{n=1}^N \log \sum_{k=1}^K p(x_n|z_{nk}=1,\Delta_k) p(z_{nk}=1|\Delta)$, we add a (arbitrary looking) distribution $q(z_{nk} = 1)$, subject to $\sum_k q(z_{nk}=1) = 1$:
$$ \mathcal{L} = \sum_{n=1}^N \log \sum_{k=1}^K \frac{q(z_{nk}=1)}{q(z_{nk}=1)} p(x_n|z_{nk}=1,\Delta_k) p(z_{nk}=1|\Delta)$$
We now have an expectation:
$$  \mathcal{L} = \sum_{n=1}^N \log \mathbb{E}_{q(z_{nk}=1)} \Bigg[\frac{1}{q(z_{nk}=1)}p(x_n|  z_{nk}=1,\ \Delta_k) p(z_{nk}=1|\Delta)\Bigg]$$
telling us how on average the parameters of $q(z_{nk})$ are able to model the likelihood of data, i.e. how on average the produced assignments explain the data disposition in the space. If we look at $q(z_{nk})$ as a posterior distribution (i.e. $q(z_{nk}|x_n)$), the equation could be interpreted as:
$$ \mathbb{E}_{\mathrm{posterior}} \Big[\frac{\mathrm{likelihood}\cdot \mathrm{prior}}{\mathrm{posterior}}\Big] = \mathbb{E}_{\mathrm{posterior}} [\mathrm{marginal~likelihood}]$$
So, applying the Jensen's Inequality:
$$ \mathcal{L} \geq \sum_{n=1}^N \mathbb{E}_{q(z_{nk}=1)} \Bigg[\log \frac{1}{q(z_{nk}=1)}p(x_n|  z_{nk}=1,\ \Delta_k) p(z_{nk}=1|\Delta)\Bigg]$$
$$ = \sum_{n=1}^N \sum_{k=1}^K q(z_{nk}=1) \log \Bigg \{\frac 1 {q(z_{nk}=1)}p(x_n|z_{nk}=1,\ \Delta_k)p(z_{nk}=1|\Delta)\Bigg\}$$
$$ = \sum_{n=1}^N \sum_{k=1}^K q(z_{nk}=1) \log  p(z_{nk}=1 |\Delta) + $$
$$  \sum_{n=1}^N \sum_{k=1}^K q(z_{nk}=1) \log  p(x_n|z_{nk}=1,\ \Delta_k) - $$
$$ \sum_{n=1}^N \sum_{k=1}^K q(z_{nk}=1) \log q(z_{nk}=1)$$
We then define $q_{nk} = q(z_{nk}=1),\ \pi_k = p(z_{nk}=1|\Delta)$ (both just scalar quantities):
$$\mathcal{L}\geq \sum_{n=1}^N \sum_{k=1}^K  q_{nk} \log \pi_k + \sum_{n=1}^N \sum_{k=1}^K q_{nk} \log p(x_n |z_{nk}=1,\ \Delta_k) - \sum_{n=1}^N \sum_{k=1}^K q_{nk} \log q_{nk}$$
and we differentiate this lower bound with respect to $q_{nk}, \ \pi_k, \ \Delta_k$, setting it to zero and obtaining an iterative update.\\
The updates for $\Delta_k, \ \pi_k$ will depend on $q_{nk}$, therefore the latter will be firstly updated.\\
This is another form of the Expectation-Maximization algorithm, derived in a different way. Again, the algorithm is very sensitive to initialization, however it is guaranteed to converge to a local maximum of the lower bound.\\
The overall algorithm can be summarized as:

\begin{algorithm}[!htbp]
    \caption{Mixture Model}
    \begin{algorithmic}%[1]
        \State Guess $\mu_k,\ \sigma^2_k, \ \pi_k \ k=1,\dots, K$
        \Loop
            \NoDo
            \State Compute $q_{nk}, \ n=1,\dots, N,\ k=1,\dots, K$
            \State Update $\mu_k, \ \sigma^2_k,\ k=1,\dots, K$
        \EndLoop
        \State \textbf{until} convergence
    \end{algorithmic}
\end{algorithm}

The clustering assignment step is governed by the probabilities $q_{nk}$ of a point $x_n$ coming from the distribution $k$:
$$ q_{nk} = p(z_{nk}=1 | x_n, \mathbf{X})$$
Each point will be then characterized by a set of $K$ probabilities.
\subsection{Gaussian Mixture Model (GMM)}
We assume the component distributions are Gaussians with diagonal covariance, $p(x_n|z_{nk}=1, \ \mu_k, \ \sigma_k^2) = \mathcal{N}(\mu, \ \sigma^2I)$.\\
When updating $\pi_k$ we are subjected to the constraint $\sum_k \pi_k = 1$, so we need to add a Lagrangian to the update term:
$$ \sum_{n,k} q_{nk} \log (\pi_k) - \lambda \Big(\sum_k \pi_k -1\Big)$$
Differentiating it and setting it to zero:
$$ \frac{\partial}{\partial \pi_k} = \frac 1 {\pi_k} \sum_n q_{nk} - \lambda = 0 \longrightarrow \sum_n q_{nk} = \lambda \pi_k$$
Summing both sides over $k$, we find 
$$\lambda = \sum_{n,k} q_{nk}$$
Therefore, substituting: 
$$\pi_k = \frac{\sum_n q_{nk}}{\sum_{n,j} q_{nj}} = \frac 1N \sum_n q_{nk}$$
\textit{Update for ${q_{nk}}$}: now the whole bound is relevant. Adding the Lagrange term $-\lambda \Big(\sum_k q_{nk}-1\Big)$ and differentiating:
$$ \frac{\partial}{\partial q_{nk}} = \log \pi_k + \log p(x_n |z_{nk} = 1, \ \Delta_k) - (\log q_{nk} +1) -\lambda $$
Re-arranging ($\lambda' = f(\lambda)$):
$$ \pi_k p(x_n |z_{nk}=1,\ \Delta_k) = \lambda' q_{nk}$$
Summing over $k$ to find $\lambda'$ and re-arranging:
$$ q_{nk} = \frac{\pi_k p(x_n|z_{nk}=1, \ \Delta_k)}{\sum_{j=1}^K \pi_j p(x_n |z_{nj}=1, \ \Delta_j)}$$
where evaluating $p(x_n|z_{nk}=1, \ \Delta_k)$ requires to plug $x_n$ inside a multi-variate Gaussian parametrized by $\mu_k, \sigma_k^2$.\\
\textit{Updates for ${\mu_k, \ \sigma_k^2}$}: these are easier (they are not subjected to any constraint). Differentiating the following and setting to zero:
$$ \sum_{n,k} q_{nk} \log \Bigg\{ \frac 1{(2\pi \sigma^2_k)^{D/2}}\exp \Bigg( - \frac 1{2\sigma^2_k}(x_n -\mu_k)^T (x_n -\mu_k)\Bigg)\Bigg\}$$
$$ \mu_k = \frac{\sum_n q_{nk}x_n}{\sum_n q_{nk}}, \quad \sigma_k^2 = \frac{\sum_n q_{nk}(x_n -\mu_k)^T(x_n -\mu_k)}{D \sum_n q_{nk}}$$
\subsection{Dirichlet Process Mixture Model (DPMM)}
Dirichlet Process Mixture Models (DPMM) \parencite{10.1214/06-BA104} is another probabilistic model, leveraging the framework up to an infinite number of mixtures. It is a \textit{non-parametric} model, meaning that the number of parameters can grow with higher amounts of data. This is especially useful when the number of clusters is not fixed, or there is the necessity to extrapolate it from data. The algorithm is Bayesian, setting a probabilistic distribution over the number of clusters.\\
To introduce it, we firstly define a \textit{Dirichlet Process} $\mathcal{DP}(\alpha, G)$, as a distribution over distributions, where:
\begin{itemize}
    \item $\alpha$ is the concentration parameter: the smaller it is, the fewer points result with high probabilities
    \item $G$ is the base distribution we want to shape through the process itself
\end{itemize}
Under a practical point of view, we can refer to a DP as a black-box having as input the distribution $G$, and as output a distribution $G'$, whose similarity to $G$ is determined by $\alpha$.
$$ G \rightarrow \boxed{\mathcal{DP}_\alpha} \rightarrow G'$$
Instead, formally, the underlying statistical process estimates the break-sticking weights $V_k$ as a recursive sampling from a  Beta distribution:
$$ \pi_k = V_k \prod_{j=1}^{k-1} (1-V_k),\ \text{ where } \ V_1, V_2, \dots \sim_{\text{iid}} \mathrm{Beta}(1,\alpha)$$
Given these weights, we can therefore approximate the original distribution $G$ by drawing some elements from it:
$$ \mu_1, \mu_2, \dots \sim_{\text{iid}} G$$
and merging them into a complete distribution:
$$ G = \sum_{k \in N} \pi_k \delta_{\mu_k}$$
where $\delta_{\mu_k}$ is a Dirac's delta (i.e. an indicator function) associating each drawn element $\mu_k$ with its corresponding weight $\pi_k$.\\
Some important properties associated to a Dirichlet Process are then:
\begin{itemize}
    \item Invariance of the expected value of the distribution $G$:
    $$ \mathbb{E}_{\mathcal{DP}(\alpha, G)}[x] = \mathbb{E}_G[x]$$
    \item Convergence at the original distribution $G$ at the bounds, i.e. as $\alpha \rightarrow \infty, \ \mathcal{DP}(\alpha, G) = G$
    \item Composition of DPs:
    $$ H \sim \mathcal{DP}(\alpha, G)$$
    $$ J \sim \mathcal{DP}(\alpha, H)$$
\end{itemize}
Finally, we can define DPs as Mixture Models. We consider a Gaussian as a base distribution $G$, and we draw some elements, that will become the cluster centers. Since we can draw an unbounded number of elements, the number of clusters is possibly unbounded. The model will be then in charge to weigh each of the infinite clusters, eventually excluding the majority of them and adapting to the data $\mathbf{X}$.\\
To explain how inference works, we resort to a typical metaphor called the ``The Chinese Restaurant''. According to this, given a certain number of tables, to generate an observation we need to sit at a table, with a probability proportional to the number of people sitting at that table. Each table will be then identified by a value (i.e. a statistical moment, which is the cluster center/mean vector in our case).\\
As we can easily see, this concept is funded on a Maximum Likelihood approach. However, in order to stick with a Bayesian approach, the framework introduces the concept a posterior (itself, a Dirichlet Process) by specifiying a prior of the same nature, allowing to reach the conjugacy.\\
A new observation will be then assigned to a new cluster proportionally to $\alpha$, however being firstly balanced against the probability of an observation given a cluster. To do this, we typically consider the data-points as Normally distributed and belonging to $K$ different clusters. At the beginning we assume there is no leakage of information about the separability of the clusters, and that we have unknown prior probabilities $\pi_k$ of a particular data point belonging to the cluster $k$. Instead of imagining that each data point is firstly assigned a cluster and then drawn from the distribution associated to that cluster, we now think of each observation being associated with parameter $\mu_k$ drawn from some discrete distribution $G$. That is, we are now treating the $\mu_k$ as being drawn from the random distribution $G$, and our prior information is incorporated into the model by the distribution over distributions $G$. Extending this reasoning to an infinite number of clusters means selecting a random prior distribution: 
$$ G(\mu_k) = \sum_{k \in N} \pi_k \delta_{\mu_k}$$
More practically, we take a random guess initially, providing a mean for each of the $K$ clusters. We incrementally change all of our assignments, unassigning and ``questioning'' each observation (in what is known as Gibbs Sampling mechanism). We do so by computing the following probability:
$$ p(z_i | z_1, \dots, z_{i-1}, z_{i+1}, z_N, \mathbf{X}, \alpha, G) = p(z_i | \{z_j\}_{j=1,\dots, N \cap j \neq i}, \mathbf{X}, \alpha)$$
representing the probability of a clustering assignment $z_i$ of a point $x_i$, given all the clustering assignments of the neighbouring points, and the parameters of the DPMM. We keep doing this for each data point, and it is proved to converge to the true posterior distribution.\\
Developing a bit the math, and defining as $\mu_k$ the means of the generating distributions, we can write the former probability (by means of the chain rule) as:
$$ p(z_i = k | \{z_j\}_{j=1,\dots, N \cap j \neq i}, \mathbf{X}, \{\mu_q\}_{q=1,\dots,K},, \alpha) $$
$$= p(z_i = k | \{z_j\}_{j=1,\dots, N \cap j \neq i}, x_i, \{x_j\}_{j=1,\dots, N \cap j \neq i}, \mu_k, \alpha)$$
$$ = p(z_i = k | \{z_j\}_{j=1,\dots, N \cap j \neq i}, \alpha) p(x_i | \mu_k,\{x_j\}_{j=1,\dots, N \cap j \neq i}) $$
$$ = \begin{cases}
\Big( \frac{N_k}{N+\alpha}\Big) \int_\mu p(x_i |\mu)p(\mu | G, \mathbf{X}) & \text{existing clusters}\\
\alpha \int _\mu p(x_i | \mu) p(\mu |G) & \text{new cluster}
\end{cases}$$
where:
\begin{itemize}
    \item we assume that all the distributions in the integral operator are Normal-like
    \item we see that the assignment weights $\Big( \frac{N_k}{N+\alpha}\Big)$ is explicitly proportional to the number $N_k$ of data points inside a particular cluster, over the total number of points $N$
\end{itemize}
It turns out that everything can be simplified up to:
$$ = \begin{cases}
\Big( \frac{N_k}{N+\alpha}\Big) \mathcal{N}\Big(x, \frac{N}{N+1}\mathbf{X}, \Sigma\Big) & \text{existing clusters}\\
\alpha\mathcal{N}(x, 0, \mathbf{I})  & \text{new cluster}
\end{cases}
$$
when assuming zero-normal distribution with identity covariance matrix.\\
The Gibbs Sampling mechanism thus represents an approximation for the estimation of the infinite-dimensional posterior distribution $p(\pi, \mu | \mathbf{X})$ (i.e., sticking with an infinite number of clusters, we would need to associate to each data point an infinite-dimensional set of prior probabilities, and cluster means as well).\\
Another solution is to resort to Variational Inference, as we do in the case of GMMs. The underlying mathematical process is the same: an approximating distribution $q$ for the posterior is introduced, and a lower bound on the marginal likelihood is optimized.\\
Despite the performant nature of the Gibbs Samplings, in our study we adopt this second approach, being less sensitive to initialization and typically faster to converge.