```python
%matplotlib notebook
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.patches import Ellipse
from ipywidgets import interact, FloatSlider
from IPython.display import display
```

# Gaussian Mixture Models


References for today's lecture:

1. Bishop, "Pattern recognition and machine learning", **Chapter 9**
1. Murphy, "Machine Learning: A Probabilistic Perspective", *MIT Press*, 2012, **Chapter 11**

## Introduction

**Unsupervised methods:**

- Contrary to supervised (regression, classification) in this scenario there are no labels/targets
- The objective of unsupervised methods is **to find structure in the data** 
- We do this to better understand the data: find hidden causes, find shared features


### Clustering

Group the data into clusters (sets) such that
- "Similar" data samples go to the same cluster
- Elements of one cluster are "different" from the elements of another cluster

Clustering requires a notion of similarity

### Latent variable models (LVM)

- Model the data using a **cause and effect** structure 
    - Observed variables (effect): Data
    - Latent variables (cause): Hidden variables that explains the observed variables 
- The objective is to infer the latent given the observed
- We assume properties for the latent variable, for example
    - Dimension of the latent is smaller than the observed &rarr; Dimensionality reduction
    - Latent variable is continuous and Gaussian distributed &rarr; **PCA**
    - Latent variables are not correlated (observed data might be)



## Mixture model

In this class we focus on LVM for clustering 

> We assume that there are $K$ clusters (sets)

We represent each data point 

$$
x_i \in \mathbb{R}^D \text{ with } i=1,\ldots, N
$$

as having a **categorical latent variable**

$$
z_i \in {1, 2, \ldots, K}
$$

This is the same as saying that $z_i$ has a categorical prior

$$
p(z_i) = \text{Cat}(\vec \pi),
$$

where $\pi_k \in [0, 1]$ are the mixing coefficients and $\sum_{k=1}^K \pi_k=1$ 

> **Summary:** the latent variable associates each data point with one of the $K$ clusters

Another use for mixture models: **Density estimation**

### Likelihood of the mixture model
    
We can write the likelihood as

$$
\begin{align}
p(x_i | \theta) &= \sum_{k=1}^K p(x_i|z_i = k, \theta) p(z_i=k| \theta) \nonumber \\
&= \sum_{k=1}^K p_k(x_i| \theta_k) \pi_k ,  \nonumber
\end{align}
$$

To continue we need to specify the family $p_k(x_i|\theta_k)$

## Gaussian mixture model (GMM)

In the GMM we further assume that each component of the mixture is a Gaussian/Normal distribution.

The we can write the likelihood as 

$$
p(x_i | \theta) = \sum_{k=1}^K \pi_k \mathcal{N}(x_i|\mu_k, \Sigma_k),
$$

where each cluster has three parameters
- The **centroid/prototype** of the cluster is the mean of the clusters: $\mu_k$
- The size and orientation of the cluster, given by the covariance: $\Sigma_k$
- The contribution of the cluster to the mixture (importance): $\pi_k$

> Fitting the GMM: **Finding the best value of these parameters** for all clusters given the data


**Generative model:** Given a value for the parameters we can sample from the GMM


Let's explore the influence of these parameters on the "shape" of the mixture


```python
def create_mixture(pis, mus, Sigmas, N=5000, rseed=None):
    """
    pis (list of floats): mixing coefficients : cluster importance
    mus (list of arrays): mean: cluster centers
    Sigmas (list of arrays): covariance: cluster shape
    N (int): number of points to create
    """
    if rseed is not None:
        np.random.seed(rseed)
    points = list()
    labels = list()
    for k in range(len(pis)):
        Nc = int(pis[k]*N)
        if len(mus[k]) > 1:
            x = np.random.multivariate_normal(mus[k], Sigmas[k], Nc)
        else:
            x = mus[k] + np.sqrt(Sigmas[k])*np.random.randn(Nc)
        points.append(x)
        labels.append([k]*Nc)
    return np.concatenate(points), np.concatenate(labels)
```


```python
data, labels = create_mixture(pis=[0.4, 0.5, 0.1], 
                              mus=[[-1], [2], [2.5]], 
                              Sigmas=[0.1, 0.5, 0.1])

fig, ax = plt.subplots(1, figsize=(6, 3), tight_layout=True)
ax.hist(data, bins=50);
```


```python
data, labels = create_mixture(pis=[0.9, 0.1], 
                              mus=[[-1, -1], [1, 1]], 
                              Sigmas=[0.5*np.eye(2), [[1, -0.9],[-0.9, 1]]])

fig, ax = plt.subplots(1, 2, figsize=(6, 3), tight_layout=True)

ax[0].scatter(data[:, 0], data[:, 1], c=labels, s=10, cmap=plt.cm.coolwarm,
              linewidths=0.1, edgecolors='k', vmin=0, vmax=1);
ax[1].hist2d(x=data[:, 0], y=data[:, 1], bins=30, cmap=plt.cm.Spectral_r);
```

## Posterior probability a.k.a. picking a cluster

Let's assume that we know the parameters of the mixture

Typically we want to infer the most probable cluster for a given data point

> This is given by the posterior probability of $z$ given $x$


Using the Bayes rule and our previous assumptions (likelihood) we have

$$
\begin{align}
r_{ik} = p(z_i = k|x_i, \theta) &= \frac{ p(x_i|z_i = k, \theta) p(z_i=k| \theta)}{p(x_i)} \nonumber \\
&= \frac{ p(x_i|z_i = k, \theta) p(z_i=k| \theta)}{ \sum_{k=1}^K p(x_i|z_i = k, \theta) p(z_i=k| \theta) } \nonumber \\
&= \frac{ \mathcal{N}(x_i|\mu_k, \Sigma_k)  \pi_k}{ \sum_{k=1}^K \mathcal{N}(x_i|\mu_k, \Sigma_k) \pi_k }, \nonumber\end{align}
$$

this posterior is also called the **responsibility** $r_{ik}$ of cluster $k$ for point $i$

Note that it is a **soft-cluster assignment**

In some cases we only need a **hard-cluster assignment**, *i.e.* the single most probable cluster, which can be obtained as

$$
c_{i} = \text{arg}\max_k r_{ik} = \text{arg}\max_k \log \mathcal{N}(x_i|\mu_k, \Sigma_k) + \log \pi_k,
$$

because of the maximum operator we can omit the evidence and use the log likelihood


```python
def jointGMM(data, pi, mu, Sigma):
    # Assume 2 dimensional data
    data_centered = data - mu
    norm = pi/(2*np.pi*np.sqrt(np.linalg.det(Sigma)))
    return norm*np.exp(-0.5*np.sum(np.multiply(data_centered, np.linalg.solve(Sigma, data_centered.T).T), axis=1))

fig, ax = plt.subplots(1, 2, figsize=(6, 3), tight_layout=True)

pis = [0.6, 0.4]; mus=[[-1, -1], [1, 1]]; 
Sigmas=[0.5*np.eye(2), [[1, -0.3],[-0.3, 1]]]
data, labels = create_mixture(pis, mus, Sigmas, rseed=0)

posterior = np.zeros(shape=(len(data), len(pis)))
for k in range(len(pis)):
    posterior[:, k] = jointGMM(data, pis[k], mus[k], Sigmas[k])
posterior = posterior/np.sum(posterior, axis=1)[:, np.newaxis]

ax[0].scatter(data[:, 0], data[:, 1], c=posterior[:, 0], s=10, 
              cmap=plt.cm.coolwarm,
              linewidths=0.1, edgecolors='k', vmin=0, vmax=1); ax[0].set_title("Soft label")
ax[1].scatter(data[:, 0], data[:, 1], c=np.argmax(posterior, axis=1), s=10, 
              cmap=plt.cm.coolwarm,
              linewidths=0.1, edgecolors='k', vmin=0, vmax=1); ax[1].set_title("Hard label");
```

## How do we fit the GMM?

We want to find the optimal point estimates of the parameters of a statistical model 

Can use **MLE**? The log likelihood in this case is 

$$
\log L(\pi,\mu,\Sigma) = \sum_{i=1}^N \log \sum_{k=1}^K \pi_k \mathcal{N}(x_i|\mu_k, \Sigma_k)
$$

**Problem** With the exception of $K=1$, we won't be able to find analytic expressions for the parameters


> In general when we don't have complete data, i.e. there are missing or hidden variables, the posterior does not factorize

In this case we have the latent/hidden variables $z_i$, we need to marginalize it (sum over k) to compute the likelihood

Other problems of the MLE/MAP

- Non-convex, it has local optima (Murphy 11.3.2)
- Very hard to compute (NP-complete hard!)
- **Unidentifiability:** There exist $K!$ configurations that yield the exact same solution (label switch)
- **Singularities:** A component may collapse to a single point ($\Sigma \to 0$ and $\log L \to \infty$)

What can we do? 

- Gradient descent with constraints
- Approximate MLE estimation: **Expectation Maximization (EM)**


## Expectation Maximization (EM)

The MLE solution can be computed if we **"fully-observe"** the data 

> We know $z_i$ for each $x_i$, *i.e.* we know which cluster owns $x_i$

The **fully-observed** log likelihood (FOLL) is

$$
\log L_{\text{FO}}(\theta) = \sum_{i=1}^N \log p(z_i) p(x_i|z_i, \theta), 
$$

but in practice we do not observe $z$ so we marginalize it by computing the **expected value** of the FOLL given the posterior

$$
\begin{align}
Q(\theta, \theta^{\text{old}}) &= \mathbb{E}_{p(z|x, \theta^{\text{old}})} \left[\log L_{\text{FO}}(\theta) \right]   \nonumber \\
&= \sum_{k=1}^K \sum_{i=1}^N p(k|x_i, \theta^{\text{old}})  \log p(k) p(x_i|\theta_k) \nonumber
\end{align}
$$

which is called the **auxiliary function** 



After this we update our parameters by **maximizing** 

$$
\theta^{\text{new}} = \text{arg}\max_\theta Q(\theta, \theta^{\text{old}}),
$$

and finally we do $\theta^{\text{old}} \leftarrow \theta^{\text{new}}$ and repeat the procedure until convergence


This iterative procedure of

1. **E-step:** Estimating the expected value of the likelihood given our current parameters
1. **M-step:** Maximizing the estimator and updating our parameters


is called **Expectation Maximization** (EM)


EM is a general algorithm:

- it can be used for other problems with non-analytical solution/non-totally observed data
- it has many variants: Incremental, Variational, Monte-Carlo (Murphy 11.4.8 and 11.4.9)

### EM updates for the GMM

Let's start by defining the auxiliary function for the GMM

$$
\begin{align}
Q(\theta, \theta^{\text{old}}) 
&= \sum_{k=1}^K \sum_{i=1}^N p(k|x_i, \theta^{\text{old}})  \log p(k) p(x_i|\theta_k)  \nonumber \\
&= \sum_{k=1}^K \sum_{i=1}^N r_{ik}^{\text{old}}  \left(\log \pi_k + \log \mathcal{N}(x_i|\mu_k, \Sigma_k) \right),\nonumber
\end{align}
$$

where $r_{ik}^{\text{old}}$ is the responsibility given the parameters of the previous iteration

Now we maximize $Q$ to get an estimator of the parameters

For the means

$$
\begin{align}
\frac{d Q}{d\mu_{k'}} &= \frac{d}{d\mu_{k'}} \sum_{k=1}^K \sum_{i=1}^N r_{ik}^{\text{old}}  \left(\log \pi_k + \log \mathcal{N}(x_i|\mu_k, \Sigma_k) \right) \nonumber \\
&= \sum_{i=1}^N r_{ik'}^{\text{old}}  \frac{d}{d\mu_{k'}}  \log \mathcal{N}(x_i|\mu_{k'}, \Sigma_{k'})\nonumber \\
&= \sum_{i=1}^N r_{ik'}^{\text{old}}  \frac{d}{d\mu_{k'}} \left( -\frac{1}{2} (x_i - \mu_{k'})^T \Sigma_{k'}^{-1} (x_i - \mu_{k'}) -\frac{D}{2} \log(2\pi) -\frac{1}{2} \log(|\Sigma_{k'}|) \right)\nonumber \\
&= -\frac{1}{2} \sum_{i=1}^N r_{ik'}^{\text{old}}  \frac{d}{d\mu_{k'}} (x_i - \mu_{k'})^T \Sigma_{k'}^{-1} (x_i - \mu_{k'})  \nonumber \\
&= \sum_{i=1}^N r_{ik'}^{\text{old}} \Sigma_{k'}^{-1} (x_i - \mu_{k'})  = 0 \nonumber
\end{align}
$$

and from here we can show that the new value of the cluster mean is given by

$$
\mu_{k}^{\text{new}}  = \frac{1}{\sum_{i=1}^N r_{ik}^{\text{old}}} \sum_{i=1}^N r_{ik}^{\text{old}} x_i
$$


Following the same procedure for the covariance we arrive to

$$
\begin{align}
\Sigma_k^{\text{new}} &= \frac{1}{\sum_{i=1}^N r_{ik}^{\text{old}}} \sum_{i=1}^N r_{ik}^{\text{old}} (x_i - \mu_k) (x_i - \mu_k)^T \nonumber \\
&= -\mu_k \mu_k^T  + \frac{1}{\sum_{i=1}^N r_{ik}^{\text{old}}} \sum_{i=1}^N r_{ik}^{\text{old}} x_i x_i^T   \nonumber
\end{align}
$$

Finally, we incorporate the constraint that $\sum_k \pi_k = 1$ to arrive to 

$$
\pi_k^{\text{new}} = \frac{\sum_{i=1}^N r_{ik}^{\text{old}}}{N}
$$


[&larr; Go back to the index](#index)

## Fitting a GMM using EM

The following algorithms summarizes the application of EM to the GMM problem

- **(1)** Set initial conditions
- **(2)** Compute cluster responsibilities for each data point using

$$
r_{ik}^{\text{old}} = \frac{ \mathcal{N}(x_i|\mu_k^{\text{old}}, \Sigma_k^{\text{old}})  \pi_k^{\text{old}}}{ \sum_{k=1}^K \mathcal{N}(x_i|\mu_k^{\text{old}}, \Sigma_k^{\text{old}}) \pi_k^{\text{old}} }
$$

- **(3)** Compute the new values for the parameters using

$$
\mu_{k}^{\text{new}}  = \frac{1}{\sum_{i=1}^N r_{ik}^{\text{old}}} \sum_{i=1}^N r_{ik}^{\text{old}} x_i,
$$

$$
\Sigma_k^{\text{new}} = -\mu_k \mu_k^T  + \frac{1}{\sum_{i=1}^N r_{ik}^{\text{old}}} \sum_{i=1}^N r_{ik}^{\text{old}} x_i x_i^T , 
$$

and

$$
\pi_k^{\text{new}} = \frac{\sum_{i=1}^N r_{ik}^{\text{old}}}{N}
$$

- **(4)** Update the parameters $\theta^{\text{old}} \leftarrow \theta^{\text{new}}$

- **(5)** If convergence criterion is met stop otherwise go back to **2**

### Example using "Old faithful" data

Iteration by iteration  (Bishop Fig. 9.8)


Animation (wikipedia)



## Practical details of EM

**Initialization and restarts** 

The EM procedure is guaranteed to converge to a *stationary point* (local minimum/maximum or saddle point)

- One can obtain a better solution using the following heuristic:
    1. Set a **number of retries**
    1. Run GMM with a random initialization
    1. If the likelihood of this solution is larger than the previous replace best
    1. Repeat until **number of retries** is achieved
- Another option is to initialize with the solution from a simpler algorithms (*e.g.* k-means)
    - This can avoid catastrophic results
    - But it can also limit exploration 
 

**Number of clusters**

The number of clusters $K$ is a parameter set *a priori* by the user

> How do we know what value of $K$ is good?

Option 1: The Elbow method, 
- Plot log likelihood (in the validation set) as a function of $K$
- log likelihood may keep decreasing with $K$ (overfit)
- but at some point the decrease is negligible: find the **elbow** of the curve



Option 2: Use Akaike Information Criterion (AIC) or Bayesian Information Criterion (BIC), or other criterion for model selection that penalizes complexity

Option 3: [The Silhouette score](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.silhouette_score.html)


**Convergence**

EM is an iterative procedure

> How many times do we repeat it?

- Check log likelihood vs iteration
- For example: stop if log likelihood increment is less than 0.1\% wrt to previous


## Special cases of GMM: The shape of the covariance

The covariance $\Sigma \in \mathbb{R}^{D\times D}$, where $D$ is the dimension of the data, defines the **shape** and **size** of the clusters

Instead of a full covariance (maximum flexibility) you could impose constraints such as
1. Spherical covariance: $\Sigma = \sigma^2 \text{I}$, $\sigma \in \mathbb{R}_{+}$
1. Diagonal covariance: $\Sigma = \vec \sigma^2 \text{I}$, $\vec \sigma \in \mathbb{R}_{+}^D$


```python
fig, ax = plt.subplots(1, 2, figsize=(6, 3), sharex=True, sharey=True, tight_layout=True)

eps = np.random.randn(2, 5000)

def update(s1, s2, rho):
    Sigma = [[s1**2, rho*s1*s2], [rho*s1*s2, s2**2]]
    data = np.dot(np.linalg.cholesky(Sigma), eps).T
    for ax_ in ax:
        ax_.cla(); 
    ax[0].scatter(data[:, 0], data[:, 1], s=10, 
                  linewidths=0.1, edgecolors='k', vmin=0, vmax=1);
    ax[1].hist2d(x=data[:, 0], y=data[:, 1], range=[[-6, 6], [-6, 6]], 
                 bins=50, cmap=plt.cm.Blues);
    
interact(update, 
         s1=FloatSlider(min=0.5, max=2.), 
         s2=FloatSlider(min=0.5, max=2.),
         rho=FloatSlider(min=-0.9, max=0.9));
```

## Gaussian mixture with [scikit-learn](https://scikit-learn.org/stable/modules/generated/sklearn.mixture.GaussianMixture.html)

From the `sklearn.mixture` scikit learn module

```python
GaussianMixture(n_components=1, # Number of clusters
                covariance_type='full', # Shape of the covariance
                tol=0.001, # termination tolerance of EM
                reg_covar=1e-06, # small number added to the diagonal of the covariance
                max_iter=100, # Max number of EM iterations
                n_init=1, # Number of random initializations
                init_params='kmeans', # Initialization for the responsabilities (kmeans or random)
                ...
               )
```

Create data to test the model


```python
data, labels = create_mixture(pis=[0.33, 0.33, 0.34], 
                              mus=[[-3, 0], [0, 0], [3 ,0]], 
                              Sigmas=[[[1, 0.9],[0.9, 1]], 
                                      [[1, -0.9],[-0.9, 1]], 
                                      [[1, 0.9],[0.9, 1]]])

fig, ax = plt.subplots(figsize=(6, 3), tight_layout=True)
ax.scatter(data[:, 0], data[:, 1], s=10, cmap=plt.cm.coolwarm,
           linewidths=0.1, edgecolors='k');
```

Find best number of clusters using the Bayesian information criterion (BIC). A model with lowest BIC is preferred

$$
BIC = K \log N - 2 \log L
$$

The BIC grows with $k$ the number of parameters and number of samples $n$ and decreases with the maximum value of the log likelihood

Try spherical, diagonal and full covariance


```python
from sklearn.mixture import GaussianMixture

maxK = 20; 
cov = 'full' 
metrics = np.zeros(shape=(maxK-1, ))
for i, K in enumerate(range(2, maxK+1)):
    gmm = GaussianMixture(n_components=K, covariance_type=cov, n_init=3)
    gmm.fit(data)
    metrics[i] = gmm.bic(data)
    
fig, ax = plt.subplots(figsize=(6, 3))
ax.plot(range(2, maxK+1), metrics, linestyle='--', marker='o');
ax.set_ylabel('BIC');
```

Inspect results

- Cluster shapes and cluster assignments
- uncertainty in the responsibilities
- Generative results



```python
gmm = GaussianMixture(n_components=10, covariance_type=cov, n_init=3)
gmm.fit(data)
hatr = gmm.predict(data); 
hatdata = gmm.sample(n_samples=10000)

fig, ax = plt.subplots(3, figsize=(6, 8), sharex=True, sharey=True, tight_layout=True)

ax[0].scatter(data[:, 0], data[:, 1], c=hatr, s=5, cmap=plt.cm.tab20,
              linewidths=0.1, edgecolors='k');
ax[0].set_title('Arg max of responsibilities')

p = gmm.predict_proba(data)
H = -np.sum(p*np.log(p), axis=1)
ax[1].scatter(data[:, 0], data[:, 1], c=H, s=10, cmap=plt.cm.Reds,
              linewidths=0.1, edgecolors='k');
ax[1].set_title('Entropy of of responsibilities')

ax[2].scatter(hatdata[0][:, 0], hatdata[0][:, 1], s=5, alpha=0.1)
ax[2].set_title('Generated data')

ax[2].scatter(gmm.means_[:, 0], gmm.means_[:, 1], s=50, c='k', marker='d')

for pos, covar, w in zip(gmm.means_, gmm.covariances_, gmm.weights_):
    if cov == 'full':
        U, s, Vt = np.linalg.svd(covar)
        angle = np.degrees(np.arctan2(U[1, 0], U[0, 0]))
        width, height = 2 * np.sqrt(s)
    elif cov == 'diag':
        angle = 0
        width, height = 2 * np.sqrt(covar)    
    elif cov == 'spherical':
        angle = 0
        width, height = 2 * np.sqrt(covar), 2 * np.sqrt(covar)     
    
    ax[2].add_patch(Ellipse(pos, 2.5*width, 2.5*height, angle, color='k', alpha=w))
    
```

## Appendix: Relation to K-means

The K-means algorithm is a GMM with the following constraints

- Spherical clusters: $\Sigma = \sigma^2 \text{I}$, $\sigma \in \mathbb{R}_{+}$
- All clusters have equal size: $\sigma^2 = \sigma_1^2 = \sigma_2^2 = \ldots = \sigma_K^2$
- Uniform prior for the mixture coefficients: $\pi_k = \frac{1}{K}$
- Hard labels are used instead of responsibilities

The multivariate normal for spherical clusters

$$
\mathcal{N}(x_i|\mu_k, \Sigma_k = \sigma^2 I) = \frac{1}{\sqrt{(2\pi)^D} \sigma} \exp \left(-\frac{1}{2\sigma^2} \|x_i - \mu_k\|^2 \right)
$$

The update rule for the responsibility 

$$
\begin{align}
\hat r_{ik}^{\text{old}} &= \text{arg} \max_k \frac{ \exp \left(-\frac{1}{2\sigma^2} \|x_i - \mu_k\|^2 \right) }{ \sum_{k=1}^K \exp \left(-\frac{1}{2\sigma^2} \|x_i - \mu_k\|^2 \right)  } \nonumber \\
&=  \text{arg} \max_k  \log \exp \left(-\frac{1}{2\sigma^2} \|x_i - \mu_k\|^2 \right)  \nonumber \\
&=  \text{arg} \min_k  \frac{1}{2\sigma^2}  \|x_i - \mu_k\|^2   \nonumber \\
&=  \text{arg} \min_k  \|x_i - \mu_k\|^2   \nonumber 
\end{align}
$$

and the means
$$
\mu_{k}^{\text{new}}  = \frac{1}{\sum_{i=1}^N \hat r_{ik}^{\text{old}}} \sum_{i=1}^N \hat r_{ik}^{\text{old}} x_i
$$


Hence the EM for K-means reduces to

1. Update "hard" responsibilities
1. Update the cluster means


## Variational mixture of Gaussians (Bishop 10.2, Murphy 21)

In EM the parameters are point estimates (similar to MLE and MAP)

But in a fully Bayesian LVM the parameters are distributions that we estimate (posterior) given some assumptions (priors)

$$
p(z, \theta |x) = \frac{p(x|z,\theta)p(z, \theta)}{\int p(x|z,\theta)p(z, \theta) dz d\theta}
$$

but computing this integral is intractable (sum over all combinations of z)

We have two choices: stochastic approximation (MCMC) or deterministic approximation

**Variational Bayes** is an example of the latter

VB approximates intractable posteriors with analytical factorized distributions



- $p(z|x)$ is the intractable posterior
- $q_\nu(z)$ is a family of "easier" factorized posterior 
- We look for $\nu$ that maximizes a lower bound on the KL divergence between posteriors
- This topic is fully explored in INFO320!


### Priors for the GMM

In GMM we set priors for the weights $\pi$, means $\mu$ and covariance $\Sigma$
- For $\mu$ we use a normal prior
- For $\Sigma$ we use Wishart (covariance) or a mix between LKJ (correlation) and half normal (variance)
- For the weights a Dirichlet prior is used
  
**Recap on Dirichlet distribution**
- Two parameters: K (number of categories) and $\alpha_k$ concentrations
- Its support is the simplex: $x_1, x_2, \ldots, x_K$, $x_k \in [0,1]$ and $\sum_k x_k = 1$
- By setting the concentration (hyperparameter) we can choose between have equally important clusters or a few active clusters


          
[Picture source](http://www.sumsar.net/blog/2015/04/the-non-parametric-bootstrap-as-a-bayesian-model/)

**Example: Variational GMM** 

- Using a large number of clusters $K$ but with a small concentration $\alpha_k = 0.01$
- Clusters with negligible concentration are not plotted



### Variational GMM with [scikit-learn](https://scikit-learn.org/stable/modules/generated/sklearn.mixture.BayesianGaussianMixture.html)

```python
BayesianGaussianMixture(...,
                        weight_concentration_prior_type='dirichlet_process',
                        weight_concentration_prior=None,
                        mean_precision_prior=None,
                        mean_prior=None,
                        degrees_of_freedom_prior=None,
                        covariance_prior=None,
                        ...
                       )
```


```python
BayesianGaussianMixture?
```


```python
from sklearn.mixture import BayesianGaussianMixture

data, labels = create_mixture(pis=[0.33, 0.33, 0.34], 
                              mus=[[-3, 0], [0, 0], [3 ,0]], 
                              Sigmas=[[[1, 0.9],[0.9, 1]], 
                                      [[1, -0.9],[-0.9, 1]], 
                                      [[1, 0.9],[0.9, 1]]])


gmm = BayesianGaussianMixture(n_components=10 , 
                              weight_concentration_prior_type='dirichlet_distribution',
                              weight_concentration_prior=10.,
                              covariance_type="full", max_iter=1000, n_init=5)
gmm.fit(data)
hatr = gmm.predict(data)
print(np.unique(hatr))
hatdata = gmm.sample(n_samples=10000)
```


```python
fig, ax = plt.subplots(3, figsize=(6, 8), sharex=True, sharey=True, tight_layout=True)

ax[0].scatter(data[:, 0], data[:, 1], c=hatr, s=5, cmap=plt.cm.tab20,
              linewidths=0.1, edgecolors='k');
ax[0].set_title('Arg max of responsibilities')

p = gmm.predict_proba(data)
H = -np.sum(p*np.log(p), axis=1)
ax[1].scatter(data[:, 0], data[:, 1], c=H, s=10, cmap=plt.cm.Reds,
              linewidths=0.1, edgecolors='k');
ax[1].set_title('Entropy of of responsibilities')

ax[2].scatter(hatdata[0][:, 0], hatdata[0][:, 1], s=5, alpha=0.1)
ax[2].set_title('Generated data')

ax[2].scatter(gmm.means_[:, 0], gmm.means_[:, 1], s=50, c='k', marker='d')

for pos, covar, w in zip(gmm.means_, gmm.covariances_, gmm.weights_):
    U, s, Vt = np.linalg.svd(covar)
    angle = np.degrees(np.arctan2(U[1, 0], U[0, 0]))
    width, height = 2 * np.sqrt(s)    
    ax[2].add_patch(Ellipse(pos, 2.5*width, 2.5*height, angle, color='k', alpha=w))
```

## More material

- EM for mixture of Bernoulli distributions: Bishop 9.3.3
- EM for mixture of student-t distributions: Murphy 11.4.5
- Online EM: Murphy 11.4.8
- GMM with priors (MAP): Murphy 11.4.2.8
- [Scikit-learn article: Setting the concentration prior](https://scikit-learn.org/stable/auto_examples/mixture/plot_concentration_prior.html)
- [Variational GMM in pymc3](https://docs.pymc.io/notebooks/gaussian-mixture-model-advi.html)





```python

```

b2f979ebcb3e20fa937c79eaef5f3312a75393c1
ipynb
Jupyter Notebook
lectures/9_mixture_models/lecture.ipynb
magister-informatica-uach/INFO337
45d7faabbd4ed5b25a575ee065551b87b097f92e
Unlicense
2021-06-12T04:07:26.000Z
2022-03-27T23:22:59.000Z
lectures/9_mixture_models/lecture.ipynb
magister-informatica-uach/INFO337
45d7faabbd4ed5b25a575ee065551b87b097f92e
Unlicense
lectures/9_mixture_models/lecture.ipynb
magister-informatica-uach/INFO337
45d7faabbd4ed5b25a575ee065551b87b097f92e
Unlicense
2019-11-07T14:49:09.000Z
2019-11-07T14:49:09.000Z
Qwen/Qwen-72B
1. YES
2. YES