Machine Learning

Machine Learning

---

5. Gaussian PDF and Maximum Likelihood Estimation (MLE)

Prof. Iacopo Masi - Computer Science Department, Sapienza University of Rome

Today's lecture

Unsupervised Learning

Density Estimator

Single Gaussian

Maximum Likelihood Estimation (MLE)

Intro to GMM (Gaussian Mixture Model)

This lecture material is taken from

• Cimi Book - Chapter 15
• Bishop - Chapter 9.2
• MLE Principle - D2L.ai
• Stanford Kmeans
• Stanford Kmeans
• Illustrations
• Code

K-means as a loss minimization problem

Like PCA, we can define a loss or cost function that can better formalize the K-means algorithm.

$$ \mathcal{L}(\mbf{\mu},y;\mbf{D}) = \sum_{i=1}^N \left| \left| \mbf{x}_i - \mbf{\mu}_{y_i} \right|\right|_2^2 $$

• Criterion for separating samples in $K$ groups of equal variance, minimizing the "inertia" or "within-cluster sum-of-squares" (within cluster variance).
• Sum of squared distances from any data point to its assigned center

Homework 2a

Given the scalar function:

$$ \mathcal{L}(\mbf{\mu};\mbf{D}) = \sum_{i=1}^N \left( {x}_i - {\mu} \right)^2 $$

Write a proof sketch that ${\mu}$ as the average of $\mbf{D}\doteq\{{x}_i\}_{i=1}^N$ is a critical point of the function. (It is also a global minimum)

Homework 2b

Given the function vector to scalar function, where the assignments $\bar{y}$ are given somehow (they cannot "move"):

$$ \mathcal{L}(\mbf{\mu},\bar{y};\mbf{D}) = \sum_{i=1}^N \left| \left| \mbf{x}_i - \boldsymbol{\mu}_{\bar{y}_i} \right|\right|_2^2 $$

Write a proof sketch that $\boldsymbol{\mu}$ as the average of the vectors $\mbf{D}\doteq\{\mbf{x}_i\}_{i=1}^N$ is a critical point of the function. (It is also a global minimum)

Hint: set the derivative/gradient of mu to zero and solve; for 2b) pay attention to how the L2 norm is defined

Unsupervised Learning

Objective and Motivation: The goal of unsupervised learning is to find hidden patterns in unlabeled data.

\begin{equation} \underbrace{\{\mathbf{x}_i\}_{i=1}^N}_{\text{known}} \sim \underbrace{\mathcal{D}}_{\text{unknown}} \end{equation}

• Unlike in supervised learning, any data point is not paired with a label.

• As you can see the unsupervised learning problem is ill-posed (which hidden patterns?) and in principle more difficult than supervised learning.

• Unsupervised learning can be thought of as "finding structure" in the data.

Discrete Random Variable

Inverse Transform Sampling

• Step 1. Transform the closed set $D$ of distances into a distribution (probability mass function - pdf) as

$$p(d^{\prime})= \frac{d^{\prime}}{\sum_{d \in D} d} $$

dist = np.array([1.18518298, 1.30917493, 1.10973212, 2.24523519, 1.01625606])

to pmf = np.array([0.17262675 0.19068668 0.16163702 0.32702769 0.14802185]) simply with

pmf = dist/dist.sum()

fig, axs = plt.subplots(1, 2)
fig.set_figheight(3)
fig.set_figwidth(15)
# PDF
axs[0].stem(pmf, linefmt='b-', markerfmt='bo', basefmt='--')
axs[0].set_title('PMF')
axs[0].set_xlabel('Index of distance')
axs[0].set_ylabel('Probability')
axs[0].set_aspect('auto')
# CUMSUM
axs[1].plot(pmf.cumsum(), 'o--')
axs[1].set_title('CDF')
axs[1].set_xlabel('Index of distance')
axs[1].set_ylabel('Cumulative Probability')
axs[1].set_aspect('auto')
plt.show()

Inverse Transform Sampling in Action

Gaussian (Normal) Distribution

• $\mu$ is the mean
• $\sigma^2$ is the variance ($k=1$)

$$ f(x) = \frac{1}{\sqrt{2 \pi \bbox[yellow, 5px]{\sigma^2}}} \exp\left( -\frac{1}{2} \bbox[yellow, 5px]{\frac{(x - \mu)^2}{\sigma^2}} \right) $$

Univariate Gaussian (1-D Gaussian)

# Univariate Gaussian (1-D Gaussian)
mu, sigma = 0, 1
x = np.arange(-3, 3, 0.01)
p = 1 / np.sqrt(2 * np.pi * sigma**2) * np.exp(-(x - mu)**2 / (2 * sigma**2)) 
_ = plt.plot(x, p);
_ = plt.title('Gaussian PDF')

Cumulative Density Function

Multivariate Gaussian (N-D Gaussian)

• $\mu$ is the mean
• $\Sigma$ is the covariance matrix
• $k$ is the dimension of the space where $x$ takes values (k=2 in this case)

$$ f(x) = \frac{1}{\sqrt{(2 \pi)^k \bbox[yellow, 5px]{\det \Sigma}}} \exp\left( -\frac{1}{2} \bbox[yellow, 5px]{(x - \mu)^T \Sigma^{-1} (x - \mu)} \right), $$

from scipy.stats import multivariate_normal
mu = [0, 0]
Sigma = [[2, 0.0],
         [0.0, 1]]
F = multivariate_normal(mu, Sigma )

from scipy.stats import multivariate_normal
X, Y = np.mgrid[-2:2:0.01, -2:2:0.01]
pos = np.dstack((X, Y))
Z = F.pdf(pos)

# plot using subplots
fig = plt.figure(figsize=(20,20))
ax1 = fig.add_subplot(1, 2, 1, projection='3d')

ax1.plot_surface(X, Y, Z, rstride=3, cstride=3, linewidth=1, antialiased=True,
                 cmap='jet')
ax1.view_init(55, -70)
# ax1.set_xticks([])
# ax1.set_yticks([])
# ax1.set_zticks([])
ax1.set_xlabel(r'$x_1$')
ax1.set_ylabel(r'$x_2$')

ax2 = fig.add_subplot(1, 2, 2, projection='3d')
ax2.contourf(X, Y, Z, zdir='z', offset=0, cmap='jet')
ax2.view_init(90, 270)

ax2.grid(False)
ax2.set_xticks([])
ax2.set_yticks([])
ax2.set_zticks([])
ax2.set_xlabel(r'$x_1$')
ax2.set_ylabel(r'$x_2$')

plt.show()

from scipy.stats import multivariate_normal
step = 0.01
X, Y = np.mgrid[-2:2:step, -2:2:step]
pos = np.dstack((X, Y))
##############################
### Gaussian PDF definition
Sigma = [[.001, 0.0],
         [0.0, .001]]
mu = [0, 0]
F = multivariate_normal(mu, Sigma )
##############################
Z = F.pdf(pos)

# plot using subplots
fig = plt.figure(figsize=(20,20))
ax1 = fig.add_subplot(1, 2, 1, projection='3d')

ax1.plot_surface(X, Y, Z, rstride=3, cstride=3, linewidth=1, antialiased=True,
                 cmap='jet')
ax1.view_init(55, -70)
# ax1.set_xticks([])
# ax1.set_yticks([])
# ax1.set_zticks([])
ax1.set_xlabel(r'$x_1$')
ax1.set_ylabel(r'$x_2$')
ax1.set_xlim3d(-2,2)
ax1.set_ylim3d(-2,2)
#ax1.set_zlim3d(0,0.15)

ax2 = fig.add_subplot(1, 2, 2, projection='3d')
ax2.contourf(X, Y, Z, zdir='z', offset=0, cmap='jet')
ax2.view_init(90, 270)

ax2.grid(False)
ax2.set_xticks([])
ax2.set_yticks([])
ax2.set_zticks([])
ax2.set_xlabel(r'$X_1$')
ax2.set_ylabel(r'$X_2$')
ax2.set_xlim(-2,2)
ax2.set_ylim(-2,2)


plt.show()

The Mahalanobis distance

The Mahalanobis distance of a point $\mbf{x}$ from a multivariate Gaussian $\mathcal{N}(\mbf{\mu},\mbf{\Sigma})$

$$ D_M(\vec{x}) = \sqrt{(\mbf{x} - \mbf{\mu})^\mathsf{T} \mathbf{\Sigma}^{-1} (\mbf{x} - \mbf{\mu})} $$

• It is a multi-dimensional generalization of the idea of measuring how many standard deviations away $\mbf{x}$ is from the mean of $\mathcal{N}(\mbf{\mu},\mbf{\Sigma})$

• This distance is zero for $\mbf{x}$ at $\mbf{\mu}$ and grows as $\mbf{x}$ moves away from the mean along each principal component axis.

• If each of these axes is re-scaled to have unit variance, then the Mahalanobis distance corresponds to the standard Euclidean distance in the transformed space. The Mahalanobis distance is thus unitless, scale-invariant, and takes into account the correlations of the data set.

Joint, Marginal, Conditional and Bayes Theorem

Joint

• The distribution above can be seen as a joint distribution $p(X,Y)$ of two random variables $X$ and $Y$
• It defines the density of two variables that co-occur together.

Joint = Marginal $\times$ Conditional

$$\underbrace{p(X,Y)}_{\text{joint}} = \underbrace{p(X)}_{\text{marginal}} \underbrace{p(Y|X)}_{\text{cond.}}$$

Which is also equal to: $$\underbrace{p(X,Y)}_{\text{joint}} = \underbrace{p(Y)}_{\text{marginal}} \underbrace{p(X|Y)}_{\text{cond.}}$$

Bayes Theorem

$$ p(X,Y) = p(X)p(Y|X) = p(Y)p(X|Y)$$$$ p(Y|X) = \frac{p(Y)p(X|Y)}{p(X)} \quad p(X) \text{ non zero} $$

Bayes Theorem (other direction)

$$ p(X,Y) = p(X)p(Y|X) = p(Y)p(X|Y)$$$$ p(X|Y) = \frac{p(X)p(Y|X)}{p(Y)} \quad p(Y) \text{ non zero} $$

You can also write the denominator more explicitly by expanding the marginal distribution (discrete case).

$$ p(Y|X) = \frac{p(Y)p(X|Y)}{p(X)} = \frac{p(Y)p(X|Y)}{\sum_{Y^{\prime}}p(Y^{\prime})p(X|Y^{\prime})} $$

Independent Random Variables

The joint probability simply becomes the product fo the marginals.

$$ p(X,Y) = p(X)p(Y|X) = p(X)p(Y) \text{ because } p(Y) = p(Y|X)$$

when they $X \perp Y $

Terminology for Statistics

Terminology for Machine Learning

Problem: Given data, learn the parameters

Assumptions

1. We have a set of observations or training data $X=\{\mbf{x}_1,\ldots,\mbf{x}_N \}$ in d-dimensional space (for now the dimensionality is not important).
2. We also assume that the generation process of the data is governed by single Gaussian $\mu,\Sigma$.
   • We do not know the values for these parameters $\mu,\Sigma$.
   • We wish to find them.
3. Important assumption: Data points that are drawn independently from the same distribution are said to be independent and identically distributed (i.i.d).

$$ X=\{\mbf{x}_1,\ldots,\mbf{x}_N \} \rightarrow \mu,\Sigma $$

The Maximum Likelihood Principle

This has a Bayesian interpretation which can be helpful to think about. Suppose that we have a model with parameters $\boldsymbol{\theta}\doteq\mu,\Sigma$ and a collection of data examples $X=\{\mbf{x}_1,\ldots,\mbf{x}_N \}$.

If we want to find the most likely value for the parameters of our model, given the data, that means we want to find

$$\mathop{\mathrm{argmax}}_ {\boldsymbol{\theta}} P(\boldsymbol{\theta}\mid X).$$

By Bayes' rule, this is the same thing as

$$ \mathop{\mathrm{argmax}}_ {\boldsymbol{\theta}} P(\boldsymbol{\theta}\mid X) = \mathop{\mathrm{argmax}}_ {\boldsymbol{\theta}} \frac{P(X \mid \boldsymbol{\theta})P(\boldsymbol{\theta})}{P(X)}. $$

• The expression $P(X)$, a parameter agnostic probability of generating the data, does not depend on $\boldsymbol{\theta}$ at all, and so can be dropped without changing the best choice of $\boldsymbol{\theta}$.
• Similarly, we may now posit that we have no prior assumption on which set of parameters is better than any others, so we may declare that $P(\boldsymbol{\theta})$ does not depend on theta either! (uninformative prior)

Thus we see that our application of Bayes' rule shows that our best choice of $\boldsymbol{\theta}$ is the maximum likelihood estimate for $\boldsymbol{\theta}$:

$$ \hat{\boldsymbol{\theta}} = \mathop{\mathrm{argmax}} _ {\boldsymbol{\theta}} P(X \mid \boldsymbol{\theta}) = \mathop{\mathrm{argmax}} _ {\boldsymbol{\theta}} P(\boldsymbol{\theta}\mid X ) $$

This is the reason why likelihood and pdf are the same but we just change the way we interpret the variables inside!

As a matter of common terminology, the probability of the data given the parameters $p(X\mid\boldsymbol{\theta})$ is referred to as the likelihood.

Maximum Likehood Estimator (MLE) for a single Gaussian

Probability Density can be interpreted in two ways:

1. as a function of datapoints, given the parameters $P$ (probability of ${\color{red}data}$, given params)
2. as a function of parameters, given the datapoints $L$ (likehood of ${\color{green}params}$, given data)

$$ p({\color{red}\bx};\boldsymbol{\mu}, \boldsymbol{\Sigma}) = \mathcal{N}({\color{red}\bx}; \boldsymbol{\mu}, \boldsymbol{\Sigma})=\frac{1}{(2 \pi)^{D / 2}} \frac{1}{|\boldsymbol{\Sigma}|^{1 / 2}} \exp \left\{-\frac{1}{2}({\color{red}\bx}-\boldsymbol{\mu})^{\mathrm{T}} \boldsymbol{\Sigma}^{-1}({\color{red}\bx}-\boldsymbol{\mu})\right\} $$

Probability Density can be interpreted in two ways:

1. as a function of datapoints, given the parameters $P$ (probability of ${\color{red}data}$, given params)
2. as a function of parameters, given the datapoints $L$ (likelihood of ${\color{green}params}$, given data)

$$ L(\mbf{\color{green}\mu}, \mbf{\color{green}\Sigma};\bx) = \frac{1}{(2 \pi)^{D / 2}} \frac{1}{|\mbf{\color{green}\Sigma}|^{1 / 2}} \exp \left\{-\frac{1}{2}(\bx - \mbf{\color{green}\mu})^{\mathrm{T}} \mbf{\color{green}\Sigma}^{-1} (\bx - \mbf{\color{green}\mu})\right\} $$

1. We have a set of observations or training data $D=\{\mbf{x}_1,\ldots,\mbf{x}_N \}$ in d-dimensional space (for now the dimensionality is not important).
2. Second equation we apply i.i.d hypothesis.

$$ p({\color{red}D};\mbf{\mu}, \mbf{\Sigma}) = p({\color{red}\{ \bx_1, \ldots, \bx_N \}};\mbf{\mu}, \mbf{\Sigma}) = $$$$ = p({\color{red}\bx_1};\mbf{\mu}, \mbf{\Sigma}) \cdot p({\color{red}\bx_2};\mbf{\mu}, \mbf{\Sigma}) \cdots p({\color{red}\bx_N};\mbf{\mu}, \mbf{\Sigma}) = $$

1. We have a set of observations or training data $D=\{\mbf{x}_1,\ldots,\mbf{x}_N \}$ in d-dimensional space (for now the dimensionality is not important).
2. First equation we apply i.i.d hypothesis.

$$ p({\color{red}D};\mbf{\mu}, \mbf{\Sigma}) = \prod_{i=1}^N p({\color{red}\bx_i};\mbf{\mu}, \mbf{\Sigma}) = \prod_{i=1}^N \mathcal{N}({\color{red}\bx_i}; \mbf{\mu}, \mbf{\Sigma}) $$

This holds also for the likelihood:

$$ L({\color{green}{\mbf{\mu}, \mbf{\Sigma}}};\mbf{x}) = \prod_{i=1}^N \mathcal{N}(\bx_i; {\color{green}{\mbf{\mu}, \mbf{\Sigma}}}) $$

• We can find ${\color{green}{\mbf{\mu}^*, \mbf{\Sigma}^*}}$ by maximizing the likelihood:

$$ {\color{green}{\mbf{\mu}^*, \mbf{\Sigma}^*}} = \arg\max_{\mbf{\mu}, \mbf{\Sigma}} L({\color{green}{\mbf{\mu}, \mbf{\Sigma}}};\mbf{x}) $$

• We can find maximum points by solving:

$$ \frac{\partial L({\color{green}{\mbf{\mu}, \mbf{\Sigma}}};\mbf{x})}{\partial {\color{green}{\mbf{\mu}}}} = \mbf{0} $$

$$ {\color{green}{\mbf{\mu}^*, \mbf{\Sigma}^*}} = \arg\max L({\color{green}{\mbf{\mu}, \mbf{\Sigma}}};\bx) $$

• We can find maximum points by taking the derivative and setting it to zero: $$ \frac{\partial L({\color{green}{\mbf{\mu}, \mbf{\Sigma}}};\bx)}{\partial {\color{green}{\mbf{\mu}}}} = \mbf{0} $$

• We can find maximum points by taking the derivative and setting it to zero: $$ \frac{\partial L({\color{green}{\mbf{\mu}, \mbf{\Sigma}}};\bx)}{\partial {\color{green}{\mbf{\Sigma}}}} = \mbf{0} $$

Minimizing the Negative Log-Likelihood

We use log because:

• It is a strictly monotonic function
• It improves numerical stability for huge values
• There is a connection with the information theory loss function
• Math will simplify because with log and exp and $\prod$ will change to $\sum$ because of log properties

$$ {\color{green}{\mbf{\mu}^*, \mbf{\Sigma}^*}} = \arg\max L({\color{green}{\mbf{\mu}, \mbf{\Sigma}}};\bx) = \arg\max \ln \big( L({\color{green}{\mbf{\mu}, \mbf{\Sigma}}};\bx) \big)$$$$ {\color{green}{\mbf{\mu}^*, \mbf{\Sigma}^*}} = \arg\min \bbox[yellow, 5px]{-\ln \big( L({\color{green}{\mbf{\mu}, \mbf{\Sigma}}};\bx) \big)}$$

Minimizing the Negative Log-Likelihood

• Replacing with log works because we need argmax, not max

$$ {\color{green}{\mbf{\mu}^*, \mbf{\Sigma}^*}} = \arg\min \bbox[yellow, 5px]{-\ln \big( L({\color{green}{\mbf{\mu}, \mbf{\Sigma}}};\bx) \big)} = -\ln \left( \prod_{i=1}^N L({\color{green}{\mbf{\mu}, \mbf{\Sigma}}};\bx_i) \right) = -\sum_{i=1}^N \ln L({\color{green}{\mbf{\mu}, \mbf{\Sigma}}};\bx_i) $$

General recipe for optimizing with MLE

Universal Template

$$ {\color{green}{\mbf{\theta}^*}} = \arg\min_{\color{green}{\mbf{\theta}}} \left[ -\sum_{i=1}^N \ln \big( L({\color{green}{\mbf{\mu}, \mbf{\Sigma}}};\bx_i) \big) \right] $$

or

$$ {\color{green}{\mbf{\theta}^*}} = \arg\max_{\color{green}{\mbf{\theta}}} \left[ \sum_{i=1}^N \ln \big( L({\color{green}{\mbf{\mu}, \mbf{\Sigma}}};\bx_i) \big) \right] $$

MLE with a single Gaussian

This holds also for the likelihood:

$$ L({\color{green}{\mbf{\mu}, \mbf{\Sigma}}};\mbf{X}) = \prod_{i=1}^N \mathcal{N}(\bx_i; {\color{green}{\mbf{\mu}, \mbf{\Sigma}}}) $$

• We can find ${\color{green}{\mbf{\mu}, \mbf{\Sigma}}}$ by maximizing the likelihood:

$$ {\color{green}{\mbf{\mu}^*, \mbf{\Sigma}^*}} = \arg\max L({\color{green}{\mbf{\mu}, \mbf{\Sigma}}};\bx) = \arg\max \ln \big( L({\color{green}{\mbf{\mu}, \mbf{\Sigma}}};\bx) \big) $$

• Using the iid assumption, we turn the product into a sum of logs:

$$ \arg\max \ln \left( \prod_{i=1}^N \mathcal{N}(\bx_i; {\color{green}{\mbf{\mu}, \mbf{\Sigma}}}) \right) = \arg\max \sum_{i=1}^N \ln \mathcal{N}(\bx_i; {\color{green}{\mbf{\mu}, \mbf{\Sigma}}}) $$

• Using log: $$ \arg\max \sum_{i=1}^N \ln \big( \mathcal{N}(\bx_i; {\color{green}{\mbf{\mu}, \mbf{\Sigma}}}) \big) $$

• And replacing Gaussian equation: $$ \arg\max \sum_{i=1}^N \ln \left( \frac{1}{(2 \pi)^{D / 2} |{\color{green}{\mbf{\Sigma}}}|^{1 / 2}} \exp \left\{ -\frac{1}{2}(\bx_i - {\color{green}{\mbf{\mu}}})^{\mathrm{T}} {\color{green}{\mbf{\Sigma}^{-1}}} (\bx_i - {\color{green}{\mbf{\mu}}}) \right\} \right) $$

• Using log property:

$$\arg\max \sum_{i=1}^N \left[ \ln \left( \frac{1}{(2 \pi)^{D / 2}} \right) + \ln \left( \frac{1}{|{\color{green}{\mbf{\Sigma}}}|^{1 / 2}} \right) + \ln \left( \exp \left\{ -\frac{1}{2}(\bx_i - {\color{green}{\mbf{\mu}}})^{\mathrm{T}} {\color{green}{\mbf{\Sigma}^{-1}}} (\bx_i - {\color{green}{\mbf{\mu}}}) \right\} \right) \right]$$

• We can simplify as:

$$ \arg\max \sum_{i=1}^N \left[ \text{const} - \frac{1}{2} \ln |{\color{green}{\mbf{\Sigma}}}| - \frac{1}{2}(\bx_i - {\color{green}{\mbf{\mu}}})^{\mathrm{T}} {\color{green}{\mbf{\Sigma}^{-1}}} (\bx_i - {\color{green}{\mbf{\mu}}}) \right] $$

• Now, we can find ${\color{green}{\mbf{\mu}}}$ so that: $$ \frac{\partial L({\color{green}{\mbf{\mu}, \mbf{\Sigma}}};\bx)}{\partial {\color{green}{\mbf{\mu}}}} = \mbf{0} $$

$$ \frac{\partial}{\partial {\color{green}{\mbf{\mu}}}} \left[ \sum_{i=1}^N \text{const} - \frac{1}{2}\ln |{\color{green}{\mbf{\Sigma}}}| - \frac{1}{2}(\bx_i - {\color{green}{\mbf{\mu}}})^{\mathrm{T}} {\color{green}{\mbf{\Sigma}^{-1}}} (\bx_i - {\color{green}{\mbf{\mu}}}) \right] = \mbf{0} $$

Let's write it as a gradient

$$ \nabla_{\color{green}{\mbf{\mu}}} \sum_{i=1}^N \left[ \text{const} - \frac{1}{2}\ln |{\color{green}{\mbf{\Sigma}}}| - \frac{1}{2}(\bx_i - {\color{green}{\mbf{\mu}}})^{\mathrm{T}} {\color{green}{\mbf{\Sigma}^{-1}}} (\bx_i - {\color{green}{\mbf{\mu}}}) \right] $$$$ = \sum_{i=1}^N \nabla_{\color{green}{\mbf{\mu}}} \left[ -\frac{1}{2}(\bx_i - {\color{green}{\mbf{\mu}}})^{\mathrm{T}} {\color{green}{\mbf{\Sigma}^{-1}}} (\bx_i - {\color{green}{\mbf{\mu}}}) \right] $$

• Change of variable:
• Compute $\nabla_{\mbf{\hat{x}}} (\mbf{\hat{x}}^{\mathrm{T}} \mbf{A} \mbf{\hat{x}})$ where $\mbf{\hat{x}} = \bx_i - {\color{green}{\mbf{\mu}}}$ and $\mbf{A} = {\color{green}{\mbf{\Sigma}^{-1}}}$

$$ \nabla_{\color{green}{\mbf{\mu}}} \sum_{i=1}^N \left[ -\frac{1}{2}(\bx_i - {\color{green}{\mbf{\mu}}})^{\mathrm{T}} {\color{green}{\mbf{\Sigma}^{-1}}} (\bx_i - {\color{green}{\mbf{\mu}}}) \right] $$

$$ \nabla_{\mbf{\hat{x}}} (\mbf{\hat{x}}^{\mathrm{T}} \mbf{A} \mbf{\hat{x}}) = \underbrace{(\mbf{A} + \mbf{A}^{\mathrm{T}})}_{\mbf{A} = \mbf{A}^{\mathrm{T}}} \mbf{\hat{x}} = 2 \mbf{A} \mbf{\hat{x}} $$

• Substituting back $\mbf{\hat{x}} = \bx_i - {\color{green}{\mbf{\mu}}}$ and $\mbf{A} = {\color{green}{\mbf{\Sigma}^{-1}}}$:

$$ \sum_{i=1}^N {\color{green}{\mbf{\Sigma}^{-1}}} (\bx_i - {\color{green}{\mbf{\mu}}}) = \mbf{0} $$

• So it is:

$$ \nabla_{\color{green}{\mbf{\mu}}} \sum_{i=1}^N \left[ -\frac{1}{2}(\bx_i - {\color{green}{\mbf{\mu}}})^{\mathrm{T}} {\color{green}{\mbf{\Sigma}^{-1}}} (\bx_i - {\color{green}{\mbf{\mu}}}) \right] = $$$$ = \sum_{i=1}^N -\frac{1}{2} \underbrace{\left[ 2 {\color{green}{\mbf{\Sigma}^{-1}}} (\bx_i - {\color{green}{\mbf{\mu}}}) \right]}_{\text{gradient}} \cdot \underbrace{(-1)}_{\text{chain rule}} $$$$ = \sum_{i=1}^N {\color{green}{\mbf{\Sigma}^{-1}}} (\bx_i - {\color{green}{\mbf{\mu}}}) $$

• Now we set the gradient equal to 0 :

• Factoring out and simplifying:

$$ {\color{green}{\mbf{\Sigma}^{-1}}} \left( \sum_{i=1}^N \bx_i \right) = {\color{green}{\mbf{\Sigma}^{-1}}} (N {\color{green}{\mbf{\mu}}}) $$

• Since the inverse covariance is non-singular, we simplify to:

$$ \sum_{i=1}^N \bx_i = N {\color{green}{\mbf{\mu}}} $$

• Final Maximum Likelihood Estimate:

$$ {\color{green}{\mbf{\mu}^*}} = \frac{1}{N} \sum_{i=1}^N \bx_i $$

Continue the proof by hand...

MLE for a single Gaussian

• Input: training data $X=\{\mbf{x}_1,\ldots,\mbf{x}_N \}$ in d-dimensional space
• Assumption: the underlying generative process is a single Gaussian
• MLE (estimate):
• $\bmf{\mu}=\frac{1}{N}\sum_i \mbf{x}_i$
• $\bmf{\Sigma}=\frac{1}{N}\sum_i (\mbf{x}_i-\bmf{\mu})(\mbf{x}_i-\bmf{\mu})^T$

Remember: MLE gives you an estimate, NOT the underlying distribution

Let's see a practical example (10 datapoints)

• Assumes we have 10 input data points training data $X=\{\mbf{x}_1,\ldots,\mbf{x}_N \}$
• True Gaussian generative process is $\mathcal{N}(0,1)$

from scipy.stats import norm
np.random.seed(1)
########################################
mu, sigma, N = 0, 1, 10  # mean and standard deviation of the generative process
# when we work in practice we have only points we do not know the generative pocess.
########################################
points = np.random.normal(mu, sigma, N)
X_plot = np.linspace(-5, 5, 1000)

true_dens = norm(mu, sigma).pdf(X_plot)

fig, ax = plt.subplots(figsize=(5, 5))
_ = ax.fill(X_plot, true_dens, fc="black",
            alpha=0.2, label="Input distribution (unknown)")
_ = ax.plot(points, -0.005 - 0.01 * np.random.random(points.shape[0]), "+k")
ax.legend()

<matplotlib.legend.Legend at 0x7f7881b48b20>

def estimate_gaussian_mle(points, X_plot, plot=True):
    # Now, we estimate with MLE in close form
    mu_mle = points.mean()
    std_mle = np.std(points, ddof=0)
    if plot:
        MLE_dens = norm(mu_mle, std_mle).pdf(X_plot)
        _ = ax.fill(X_plot, MLE_dens, fc="red",
                    alpha=0.2, label="estimated")
    return mu_mle, std_mle


mu_mle, std_mle = estimate_gaussian_mle(points, X_plot)
print(f'Estimated ({mu_mle}, {std_mle}) vs Ground-truth ({mu}, {sigma})')

Estimated (-0.09714089080609985, 1.190898552063902) vs Ground-truth (0, 1)

#help(np.std) #toggle to get help on the standard deviation function in numpy

fig, ax = plt.subplots(figsize=(5, 5))
_ = ax.fill(X_plot, true_dens, fc="black",
            alpha=0.2, label="input distribution")
_ = ax.plot(points, -0.005 - 0.01 * np.random.random(points.shape[0]), "+k")
mu_mle, std_mle = estimate_gaussian_mle(points, X_plot)
print(f'Estimated ({mu_mle}, {std_mle}) vs Ground-truth ({mu}, {sigma})')

Estimated (-0.09714089080609985, 1.190898552063902) vs Ground-truth (0, 1)

Let's see a practical example (100 datapoints)

• Assumes we have 100 input data points training data $X=\{\mbf{x}_1,\ldots,\mbf{x}_N \}$
• True Gaussian generative process is $\mathcal{N}(0,1)$

from scipy.stats import norm
np.random.seed(1)
########################################
mu, sigma, N = 0, 1, 100  # mean and standard deviation of the generative process
# when we work in practice, we have only points and we do not know the generative process.
########################################
# when we work in practice, we have only points and we do not know the generative process.
points = np.random.normal(mu, sigma, N)
X_plot = np.linspace(-5, 5, 1000)

true_dens = norm(mu, sigma).pdf(X_plot)

fig, ax = plt.subplots(figsize=(5, 5))
_ = ax.fill(X_plot, true_dens, fc="black",
            alpha=0.2, label="input distribution")
_ = ax.plot(points, -0.005 - 0.01 * np.random.random(points.shape[0]), "+k")

fig, ax = plt.subplots(figsize=(5, 5))
_ = ax.fill(X_plot, true_dens, fc="black",
            alpha=0.2, label="input distribution")
_ = ax.plot(points, -0.005 - 0.01 * np.random.random(points.shape[0]), "+k")
mu_mle, std_mle = estimate_gaussian_mle(points, X_plot)
print(f'Estimated ({mu_mle}, {std_mle}) vs Ground-truth ({mu}, {sigma})')

Estimated (0.060582852075698704, 0.885156213831585) vs Ground-truth (0, 1)

Better than before but the std. deviation is still wrong

Let's see a practical example (1000 data points)

• Assumes we have 1000 input data points training data $X=\{\mbf{x}_1,\ldots,\mbf{x}_N \}$
• True Gaussian generative process is $\mathcal{N}(0,1)$

from scipy.stats import norm
np.random.seed(1)
########################################
mu, sigma, N = 0, 1, 1000  # mean and standard deviation of the generative process
# when we work in practice, we have only points and we do not know the generative process.
########################################
# when we work in practice, we have only points and we do not know the generative process.
points = np.random.normal(mu, sigma, N)
X_plot = np.linspace(-5, 5, 1000)

true_dens = norm(mu, sigma).pdf(X_plot)

fig, ax = plt.subplots(figsize=(5, 5))
_ = ax.fill(X_plot, true_dens, fc="black",
            alpha=0.2, label="input distribution")
_ = ax.plot(points, -0.005 - 0.01 * np.random.random(points.shape[0]), "+k")

fig, ax = plt.subplots(figsize=(5, 5))
_ = ax.fill(X_plot, true_dens, fc="black",
            alpha=0.2, label="input distribution")
_ = ax.plot(points, -0.005 - 0.01 * np.random.random(points.shape[0]), "+k")
mu_mle, std_mle = estimate_gaussian_mle(points, X_plot)
print(f'Estimated ({mu_mle}, {std_mle}) vs Ground-truth ({mu}, {sigma})')

Estimated (0.03881247615960185, 0.9810041339322116) vs Ground-truth (0, 1)

With 1000 training points seems to be OK

But what about different shapes?

MLE tends to underestimate the deviation of the Gaussian

Machine Learning

---

5. Maximum Likelihood Estimation, Gaussian Mixture Model, Clustering

Prof. Iacopo Masi - Computer Science Department, Sapienza University of Rome

Today's lecture

Unsupervised Learning

Density Estimator

Single Gaussian

Maximum Likelihood Estimation (MLE)

Intro to GMM (Gaussian Mixture Model)

This lecture material is taken from

• Cimi Book - Chapter 15
• Bishop - Chapter 9.2
• MLE Principle - D2L.ai
• Stanford Kmeans
• Stanford Kmeans
• Illustrations
• Code

K-means as a loss minimization problem

Like PCA, we can define a loss or cost function that can better formalize the K-means algorithm.

$$ \mathcal{L}(\mbf{\mu},y;\mbf{D}) = \sum_{i=1}^N \left| \left| \mbf{x}_i - \mbf{\mu}_{y_i} \right|\right|_2^2 $$

• Criterion for separating samples in $K$ groups of equal variance, minimizing the "inertia" or "within-cluster sum-of-squares" (within cluster variance).
• Sum of squared distances from any data point to its assigned center

Homework 2a

Given the scalar function:

$$ \mathcal{L}(\mbf{\mu};\mbf{D}) = \sum_{i=1}^N \left( {x}_i - {\mu} \right)^2 $$

Write a proof sketch that ${\mu}$ as the average of $\mbf{D}\doteq\{{x}_i\}_{i=1}^N$ is a critical point of the function. (It is also a global minimum)

Homework 2b

Given the function vector to scalar function, where the assignments $\bar{y}$ are given somehow (they cannot "move"):

$$ \mathcal{L}(\mbf{\mu},\bar{y};\mbf{D}) = \sum_{i=1}^N \left| \left| \mbf{x}_i - \boldsymbol{\mu}_{\bar{y}_i} \right|\right|_2^2 $$

Write a proof sketch that $\boldsymbol{\mu}$ as the average of the vectors $\mbf{D}\doteq\{\mbf{x}_i\}_{i=1}^N$ is a critical point of the function. (It is also a global minimum)

Hint: set the derivative/gradient of mu to zero and solve; for 2b) pay attention to how the L2 norm is defined

Unsupervised Learning

Objective and Motivation: The goal of unsupervised learning is to find hidden patterns in unlabeled data.

\begin{equation} \underbrace{\{\mathbf{x}_i\}_{i=1}^N}_{\text{known}} \sim \underbrace{\mathcal{D}}_{\text{unknown}} \end{equation}

• Unlike in supervised learning, any data point is not paired with a label.

• As you can see the unsupervised learning problem is ill-posed (which hidden patterns?) and in principle more difficult than supervised learning.

• Unsupervised learning can be thought of as "finding structure" in the data.

Minimizing the Negative Log-Likelihood

We use log because:

• It is a strictly monotonic function
• It improves numerical stability for huge values
• There is a connection with the information theory loss function
• Math will simplify because with log and exp and $\prod$ will change to $\sum$ because of log properties

$$ \boldsymbol{\mu}^*, \boldsymbol{\Sigma}^* = \arg\max L(\boldsymbol{\mu}, \boldsymbol{\Sigma};\mbf{x}) = \arg\max \ln \big( L(\boldsymbol{\mu}, \boldsymbol{\Sigma};\mbf{x}) \big)$$$$ \boldsymbol{\mu}^*, \boldsymbol{\Sigma}^* = \arg\min \bbox[yellow, 5px]{-\ln \big( L(\boldsymbol{\mu}, \boldsymbol{\Sigma};\mbf{x}) \big)}$$

Minimizing the Negative Log-Likelihood

• Replacing with log works because we need argmax, not max

$$ \boldsymbol{\mu}^*, \boldsymbol{\Sigma}^* = \arg\min \bbox[yellow, 5px]{-\ln \big( L(\boldsymbol{\mu}, \boldsymbol{\Sigma};\mbf{x}) \big)} =-\ln \big( \prod_{i=1}^N L(\boldsymbol{\mu}, \boldsymbol{\Sigma};\mbf{x}_i) \big)=-\sum_{i=1}^N\ln L(\boldsymbol{\mu}, \boldsymbol{\Sigma};\mbf{x}_i)$$

General recipe for optimizing with MLE

Universal Template

$$ \boldsymbol{\theta}^* = \arg\min_{\boldsymbol{\theta}} -\sum_{i=1}^N\ln \big( L(\boldsymbol{\mu}, \boldsymbol{\Sigma};\mbf{x}_i)\big)$$

or

$$ \boldsymbol{\theta}^* = \arg\max_{\boldsymbol{\theta}} \sum_{i=1}^N\ln \big( L(\boldsymbol{\mu}, \boldsymbol{\Sigma};\mbf{x}_i)\big) $$

Today's lecture

Gaussian Mixture Model (GMM)

Density Estimator with GMM

Generate data with GMM

This lecture material is taken from

• Cimi Book - Chapter 15
• Bishop - Chapter 9.2
• Stanford Kmeans
• Stanford Kmeans
• Illustrations
• Code

Unsupervised Learning

Objective and Motivation: The goal of unsupervised learning is to find hidden patterns in unlabeled data.

\begin{equation} \underbrace{\{\mathbf{x}_i\}_{i=1}^N}_{\text{known}} \sim \underbrace{\mathcal{D}}_{\text{unknown}} \end{equation}

• Unlike in supervised learning, any data points is not paired with a label.

• As you can see the unsupervised learning problem is ill-posed (which hidden patterns?) and in principle more difficult than supervised learning.

• Unsupervised learning can be thought of as "finding structure" in the data.

Problem: What if the data is generated by a multi-modal distribution?

• For example: you are modeling the background for a video-surveillance application, and you have assumed that the background more or less behaves in a single mode during the whole day plus some minor variations.
  • For example, the background is always blackish with some minor variations of "black" given by less or more light.
• You then decide to estimate the parameters for each pixel of a single Gaussian using MLE, given a training set of background videos.
• At test-time, you can estimate how probable is a new pixel to be part of the background or foreground by simply estimating the probability of this new pixel to be part of your MLE estimate for that pixel.

Quick Note on Visualization in 1D/2D with Histograms

Gaussian

• We decide using a Gaussian distribution defined by $\theta=\bmu,\bsigma$
• Even if the training set increases, we are still going to store $\theta=\bmu,\bsigma$

Histogram

• We could decide to model it with a normalized histogram

• In practice, the histogram technique can be useful for obtaining a quick visualization of data in one or two dimensions but is unsuited to most density estimation applications

A single Gaussian for each pixel

Why model the background?

Problem: what if the data is multi-modal?

• By inspecting better the data, you realize that the single mode assumption is NOT correct.
• at 3 pm the background is black-ish but then suddenly in the data you have a tree that continuously moves the leaves so that those will be also part of the background
• Your single-mode program will detect each time the leaves as a false alarm and the video-surveillance app won't work well.

What can we do now?

A single Gaussian will not work properly

Mixture of Gaussians (MoG)

• A Gaussian mixture model represents a distribution as $$ p(\mathbf{x})=\sum_{k=1}^{K} \pi_{k} \mathcal{N}\left(\mathbf{x} \mid \mu_{k}, \Sigma_{k}\right) $$ with $\pi_{k}$ the mixing coefficients, where: $$ \sum_{k=1}^{K} \pi_{k}=1 \quad \text { and } \quad \pi_{k} \geq 0 \quad \forall k $$
• GMM is a density estimator

Mixture of Gaussians (MoG)

• A Gaussian mixture model represents a distribution as $$ p(\mathbf{x})=\sum_{k=1}^{K} \bbox[yellow, 5px]{\pi_{k}} \mathcal{N}\left(\mathbf{x} \mid \bbox[yellow, 5px]{\mu_{k}, \Sigma_{k}}\right) $$ with $\pi_{k}$ the mixing coefficients, where: $$ \sum_{k=1}^{K} \pi_{k}=1 \quad \text { and } \quad \pi_{k} \geq 0 \quad \forall k $$
• GMM is a density estimator like the single Gaussian case but handles multi-modal distributions
• You can interpret it as a linear combination of Gaussians defined by the $\pi_{k}$

Mixture of Gaussians (MoG)

• We will use a Mixture of Gaussians (MoG) to indicate the distribution used.
• We will use GMM (Gaussian Mixture Model) to indicate the machine learning method behind it.

• It is very powerful and expressive: GMMs are universal approximators of densities (if you have enough Gaussians).
• Even diagonal GMMs are universal approximators

[Slide credit: K. Kutulakos]

[Slide credit: K. Kutulakos]

Linear combination of 3 Gaussians in 2D space

• You can think of the first picture from left as $p(\mbf{x},z)$ the joint density.
• You can think of the second picture from left as $p(\mbf{x})$ density of all data (marginal).

Interpretation as a Generative Model

From the sum and product rules, we can express the data density $p(\mbf{x})$ as the marginal density where the marginalization is over the mode/Gaussians we have: $$ p(\mathbf{x})=\sum_{k=1}^{K} \underbrace{p(k)}_{\text{select mode}} \underbrace{p(\mathbf{x} \mid k)}_{\text{gen. point}} $$ which is equivalent to MoG equation, in which we can view $\pi_{k}=p(k)$ as the prior probability of picking the $k^{\text {th }}$ component, and the density $\mathcal{N}\left(\mathbf{x} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)=p(\mathbf{x} \mid k)$ as the probability of $\mathbf{x}$ conditioned on $k$.

Interpretation as a Generative Model

How can we generate points assuming we have the parameters of the model?

Interpretation as a Generative Model

How can we generate points assuming we have a model with $3$ Gaussian where mixing coefficients are $\pi=[0.5, 0.3, 0.2]$;

1. Gaussian k=1 is $\bmf{\mu}=[0,0]; \bmf{\Sigma}=[0.5, 0; 0, 0.15]$
2. Gaussian k=2 is $\bmf{\mu}=[3,2]; \bmf{\Sigma}=[1.25, 0; 0, 0.75]$
3. Gaussian k=3 is $\bmf{\mu}=[-2,-2]; \bmf{\Sigma}=[1, -0.74; -0.74, 1]$

Generative Model

A possible algorithm is:

1. sample which Gaussian $\hat{k}$ by inverse transform sampling on the mixing coefficients $\pi=[0.5, 0.3, 0.2]$
2. Once you know $\hat{k}$, then select the params of the appropriate gaussian and sample $\mbf{x}\sim\mathcal{N}_{\hat{k}}(\bmf{\mu},\bmf{\Sigma})$

Probabilistic Interpretation

A possible algorithm is:

1. sample $\mbf{z} \sim p(\mbf{z})$ from $\pi=[0.5, 0.3, 0.2]$
2. then sample $\mbf{x}$ from $p(\mbf{x}|\mbf{z})$ that is $\mbf{x}\sim\mathcal{N}_{\mbf{z}}(\bmf{\mu},\bmf{\Sigma})$

import matplotlib.pyplot as plt
np.random.seed(0)

def inverse_sampling(pmf):
    # import random
    # random.choices(range(1,k), pmf, k=1000)
    return np.argmin((np.random.rand(1)[:, None] > pmf.cumsum()), axis=1)


################### model params #######################
N_samples = 1000  # we sample this amount of points
z_to_gaussian = {0: ([0, 0], [[0.5, 0], [0, 0.15]], 'red'),
                 1: ([3, 2], [[1.25, 0], [0, 0.75]], 'green'),
                 2: ([-2, -2], [[1, -0.74], [-0.74, 1]], 'blue'),
                 }
mixing = np.array([0.05, 0.05, 0.9])  # mixing coefficients
########################################################
for _ in range(0, N_samples):  # for each sample
    z = inverse_sampling(mixing)  # sample k using mixing coefficients
    *normal, color = z_to_gaussian[z[0]]  # now sample the data from  N_k
    x, y = np.random.multivariate_normal(
        *normal, 1).T  # sample 1 point at a time for clarity
    plt.plot(x, y, '.', color=color)
plt.title('p(x,z)')
plt.axis('equal')
plt.show()

import matplotlib.pyplot as plt
np.random.seed(0)

def inverse_sampling(pmf):
    return np.argmin((np.random.rand(1)[:, None] > pmf.cumsum()), axis=1)


################### model params #######################
N_samples = 1000  # we sample this amount of points
z_to_gaussian = {0: ([0, 0], [[0.5, 0], [0, 0.15]], 'red'),
                 1: ([3, 2], [[1.25, 0], [0, 0.75]], 'green'),
                 2: ([-2, -2], [[1, -0.74], [-0.74, 1]], 'blue'),
                 }
mixing = np.array([0.2, 0.5, 0.3])  # mixing coefficients
########################################################
for _ in range(0, N_samples):  # for each sample
    z = inverse_sampling(mixing)  # sample k using mixing coefficients
    *normal, color = z_to_gaussian[z[0]]  # now sample the data from  N_k
    x, y = np.random.multivariate_normal(
        *normal, 1).T  # sample 1 point at time for clarity
    plt.plot(x, y, '.', color='black')
plt.axis('equal')
plt.title('p(x)')
plt.show()

How to fit a Mixture of Gaussians (MoG)

• Remember we always start with a data set of data points $X=\{\mbf{x}_1,\ldots,\mbf{x}_N \}$
• if we knew the "labels" of the Gaussian and their association with the data points, it would be EASY to recover the MoG
• Answer the question from which Gaussian is this point sampled? Alas, we do not have the answer to this question.

Supervised Setting (uncommon)

Let's assume $k=3$ Gaussian and $X=\{\mbf{x}_1,\ldots,\mbf{x}_N \}$ and $Y=\{{y}_1,\ldots,{y}_N \}$ where $Y$ denotes an index from 0 to k-1 that selects the Gaussian for that point.

If we had, we could do:

• Estimate $K$ parameters $\bmf{\mu}_k,\bmf{\Sigma}_k$ by:
  • Performing MLE on all points but selecting those of class $k=0$.
    • Obtain $\bmf{\mu}_0,\bmf{\Sigma}_0$. The red Gaussian before.
  • Performing MLE on all points but selecting those of class $k=1$.
    • Obtain $\bmf{\mu}_1,\bmf{\Sigma}_1$. The green Gaussian before.
  • Performing MLE on all points but selecting those of class $k=2$.
    • Obtain $\bmf{\mu}_2,\bmf{\Sigma}_2$. The blue Gaussian before.

$\bbox[yellow, 5px]{Problem:}$ we cannot assume to know the association of the points to Gaussians. (i..e we do not know $Y$).

How to fit a Mixture of Gaussians (MoG)

We can do MLE on a Mixture of Gaussians:

• Maximum likelihood maximizes $$ \ln p(\mathbf{X} \mid \pi, \mu, \Sigma)=\sum_{n=1}^{N} \ln \left(\sum_{k=1}^{K} \pi_{k} \mathcal{N}\left(\mathbf{x}^{(n)} \mid \mu_{k}, \Sigma_{k}\right)\right) $$ w.r.t $\Theta=\left\{\pi_{k}, \mu_{k}, \Sigma_{k}\right\}$

Problems:

• Singularities: Arbitrarily large likelihood when a Gaussian explains a single point
• Identifiability: the solution is up to permutations
• No close-form solution.

What is a singularity?

A (trivial) way of thinking about a data set of points to fit a parametric model or to estimate a density is to:

1. Fit a $\bmf{\delta}$ Delta Dirac function to each point.
2. Scale each delta for the inverse of the number of points we have (to adjust the total area)
3. The sum over the deltas will sum up to 1 and thus will be a probability distribution

This is related to the idea of overfitting the data (more on overfitting in supervised learning).

What is a singularity?

• Gaussian with infinitesimal std. deviation centered on the data point (i.e. Delta Dirac on the point)

$$ \mathcal{N}\left(\mathbf{x}_{n} \mid \mathbf{x}_{n}, \sigma_{j}^{2} \mathbf{I}\right)=\frac{1}{(2 \pi)^{1 / 2}} \frac{1}{\sigma_{j}} $$

Trivial Solution: Data Density is a Delta Dirac for each point

$$ p(\mbf{x}) = \frac{1}{N} \sum_{i=1}^N \delta(\mbf{x}-\mbf{x}_i) $$

How to fit a Mixture of Gaussians (MoG)

We can do MLE on a Mixture of Gaussians:

• Maximum likelihood maximizes $$ \ln p(\mathbf{X} \mid \pi, \mu, \Sigma)=\sum_{n=1}^{N} \ln \left(\sum_{k=1}^{K} \pi_{k} \mathcal{N}\left(\mathbf{x}^{(n)} \mid \mu_{k}, \Sigma_{k}\right)\right) $$ w.r.t $\Theta=\left\{\pi_{k}, \mu_{k}, \Sigma_{k}\right\}$

Assumptions:

• No Singularities
• No close-form solution. So we will implement an iterative solution.

How to fit a Mixture of Gaussians (MoG)

Trick 1: Use a latent variable $z$ to indicate the Gaussian selected

Trick 2: Express the data density as the marginalization of the latent variable

$$ p(\mathbf{x})=\sum_z \underbrace{p(\mbf{z})}_{\text{select mode}} \underbrace{p(\mathbf{x} \mid \mbf{z})}_{\text{gen. point}} $$

Latent Variable

• Some model variables may be unobserved, either at training, at test time, or both
• If occasionally unobserved they are missing, e.g., undefined inputs, missing class labels, erroneous targets

• Variables that are always unobserved are called latent variables, or sometimes hidden variables

• We may want to intentionally introduce latent variables to model complex dependencies between variables – this can simplify the model

• In a mixture model, the identity of the component that generated a given data point is a latent variable

Expectation Maximization (EM) for MoG

• EM is an elegant and powerful method for finding maximum likelihood solutions for models with latent variables in an interactive way.
• It can be applied to other models not only to MoG.
• You will find immediately that will resemble K-means

It is made of two steps:

1. Expectation step (E-step):
   • We assume we have a good initialization/version of the weights $\{\bmf{\mu}_k,\bmf{\Sigma}_k,\bmf{\pi_k}\}_k$ (same problem of initialization has K-means).
   • To adjust the parameters, we must first solve the inference problem. Which Gaussian generated each data point?
   • Compute the posterior probability that each Gaussian generates each data point (as this is unknown to us)
   • Given a data point, what is the probability $\mbf{x}$ comes from Gaussian $\mbf{z}==k$? $$ \text{Compute} \quad p(\mbf{z} | \mbf{x}; \bmf{\mu},\bmf{\Sigma},\bmf{\pi})$$

1. Maximization step (M-step):
   • Assuming that the data really was generated this way, change the parameters of each Gaussian to maximize the probability that it would generate the data it is currently responsible for
   • Each Gaussian gets a certain amount of posterior probability for each data point
   • At the optimum we shall satisfy $\frac{\partial \ln p(\mathbf{X} \mid \pi, \mu, \Sigma)}{\partial \Theta}=0$
   • We can derive closed-form updates for all parameters **in this step only**

E-step: Posterior Prob. and Responsibilities

Given a data point, what is the probability comes from Gaussian $\mbf{z}==k$?

• We apply Bayes Theorem in the middle

$$ p(\bz | \bx ) = \frac{p(\bx | \bz) p(\bz)}{p(\bx)} = \frac{p(\bx | \bz) p(\bz)}{\sum_z p(\bx | \bz) p(\bz)} $$

Given a data point, what is the probability comes from Gaussian $\mbf{z}==k$?

$$ \gamma_{k} = p(\bz=k | \bx) = \frac{p(\bx | \bz=k) p(\bz=k)}{p(\bx)} = \frac{p(\bx | \bz=k) p(\bz=k)}{\sum_j^K p(\bx | \bz=j) p(\bz=j)} $$

$$ \gamma_{k} = \frac{\bbox[yellow, 5px]{p(\bx | \bz=k) p(\bz=k)}}{\sum_j^K p(\bx | \bz=j) p(\bz=j)} $$

$$\underbrace{p(\bx | \bz=k)}_{\text{generation}}\underbrace{p(\bz=k)}_{\text{selection}}$$

$$ \underbrace{p(\bx | \bz=k)}_{\text{generation}}\underbrace{p(\bz=k)}_{\text{selection}} = \mathcal{N}(\bx;\bmu_k,\bsigma_k)\cdot\pi_k $$

Responsibilities

$$\gamma_{k} =\frac{ \mathcal{N}(\bx;\bmu_k,\bsigma_k)\cdot\pi_k} {\sum_j\mathcal{N}(\bx;\bmu_j,\bsigma_j)\cdot\pi_j} $$

• The responsibility that component $k$ takes for ‘explaining’ the observation x
• For a data point, $\bmf{\gamma}$ is $1\times K$ i.e. [0.8, 0.1, 0.1] for 3 Gaussians
  • It means that 80% of the data is explained by 1st Gaussian.
  • It means that 10% of the data is explained by 2nd Gaussian.
  • It means that 10% of the data is explained by 3rd Gaussian.

Responsibilities are on the right-most figure

Responsibilities are Soft-Assignments

• Soft-Assignment means an association of points to clusters
• For $N$ data point is a matrix $N\times K$
• Note the strong similarities in computing distances in K-means (but here we do soft-assignments)!

Maximization Step

• Log-likelihood: $$ \ln p(\mathbf{X} \mid \pi, \mu, \Sigma)=\sum_{n=1}^{N} \ln \left(\sum_{k=1}^{K} \pi_{k} \mathcal{N}\left(\mathbf{x}^{(n)} \mid \mu_{k}, \Sigma_{k}\right)\right) $$
• Set derivatives to zero: $$ \frac{\partial \ln p(\mathbf{X} \mid \pi, \mu, \Sigma)}{\partial \mu_{k}}=0$$
• for $\bmu_k,\bsigma_k,\pi_k$

Maximization Step

$$ -\sum_{n=1}^{N} \underbrace{\frac{\pi_{k} \mathcal{N}\left(\mathbf{x}^{(n)} \mid \mu_{k}, \Sigma_{k}\right)}{\sum_{j=1}^{K} \pi_{j} \mathcal{N}\left(\mathbf{x} \mid \mu_{j}, \Sigma_{j}\right)}}_{\gamma_{nk}}\Sigma_{k}^{-1}\left(\mathbf{x}^{(n)}-\mu_{k}\right) = 0 $$

• We used: $$ \mathcal{N}(\mathbf{x} \mid \mu, \Sigma)=\frac{1}{\sqrt{(2 \pi)^{d}|\Sigma|}} \exp \left(-\frac{1}{2}(\mathbf{x}-\mu)^{T} \Sigma^{-1}(\mathbf{x}-\mu)\right) $$ and: $$ \frac{\partial\left(\mathbf{x}^{T} A \mathbf{x}\right)}{\partial \mathbf{x}}=\mathbf{x}^{T}\left(A+A^{T}\right) $$

Maximization Step $\bmu_k$

• This gives $$ \mu_{k}=\frac{1}{N_{k}} \sum_{n=1}^{N} \gamma_{k}^{(n)} \mathbf{x}^{(n)} $$ with $N_{k}$ the effective number of points in cluster $k$ $$ N_{k}=\sum_{n=1}^{N} \gamma_{k}^{(n)} $$
• We just take the center-of-gravity of the data that the Gaussian is responsible for
• Just like in K-means, except the data is weighted by the posterior probability of the Gaussian.
• Guaranteed to lie in the convex hull of the data (Could be a big initial jump)

Maximization Step $\bsigma_k$

• We can similarly get the expression for the variance $$ \Sigma_{k}=\frac{1}{N_{k}} \sum_{n=1}^{N} \gamma_{k}^{(n)}\left(\mathbf{x}^{(n)}-\mu_{k}\right)\left(\mathbf{x}^{(n)}-\mu_{k}\right)^{T} $$
• We can also minimize w.r.t the mixing coefficients $$ \pi_{k}=\frac{N_{k}}{N}, \quad \text { with } \quad N_{k}=\sum_{n=1}^{N} \gamma_{k}^{(n)} $$
• The optimal mixing proportion to use (given these posterior probabilities) is just the fraction of the data that the Gaussian gets responsibility for.
• Note that this is not a closed-form solution of the parameters, as they depend on the responsibilities $\gamma_{k}^{(n)}$, which are complex functions of the parameters
• But we have a simple iterative scheme to optimize

EM Algorithm for GMM

• Initialize the means $\mu_{k}$, covariances $\Sigma_{k}$ and mixing coefficients $\pi_{k}$
• Iterate until convergence:
• E-step: Evaluate the responsibilities given current parameters $$ \gamma_{k}^{(n)}=p\left(z^{(n)} \mid \mathbf{x}\right)=\frac{\pi_{k} \mathcal{N}\left(\mathbf{x}^{(n)} \mid \mu_{k}, \Sigma_{k}\right)}{\sum_{j=1}^{K} \pi_{j} \mathcal{N}\left(\mathbf{x}^{(n)} \mid \mu_{j}, \Sigma_{j}\right)} $$
• M-step: Re-estimate the parameters given current responsibilities $$ \begin{aligned} &\mu_{k}=\frac{1}{N_{k}} \sum_{n=1}^{N} \gamma_{k}^{(n)} \mathbf{x}^{(n)} \\ &\Sigma_{k}=\frac{1}{N_{k}} \sum_{n=1}^{N} \gamma_{k}^{(n)}\left(\mathbf{x}^{(n)}-\mu_{k}\right)\left(\mathbf{x}^{(n)}-\mu_{k}\right)^{T} \\ &\pi_{k}=\frac{N_{k}}{N} \text { with } N_{k}=\sum_{n=1}^{N} \gamma_{k}^{(n)} \end{aligned} $$
• Evaluate log likelihood and check for convergence $$ \ln p(\mathbf{X} \mid \pi, \mu, \Sigma)=\sum_{n=1}^{N} \ln \left(\sum_{k=1}^{K} \pi_{k} \mathcal{N}\left(\mathbf{x}^{(n)} \mid \mu_{k}, \Sigma_{k}\right)\right) $$

Let's try this out on 2D sample dataset

We assume we know the generative process

...but the input is this!

import matplotlib.pyplot as plt
np.random.seed(0)

def inverse_sampling(pmf):
    return np.argmin((np.random.rand(1)[:, None] > pmf.cumsum()), axis=1)


################### model params #######################
N_samples = 1000  # we sample this amount of points
z_to_gaussian = {0: ([0, 0], [[0.5, 0], [0, 0.15]], 'red'),
                 1: ([3, 2], [[1.25, 0], [0, 0.75]], 'green'),
                 2: ([-2, -2], [[1, -0.74], [-0.74, 1]], 'blue'),
                 }
mixing = np.array([0.2, 0.5, 0.3])  # mixing coefficients
########################################################
for _ in range(0, N_samples):  # for each sample
    z = inverse_sampling(mixing)  # sample k using mixing coefficients
    *normal, color = z_to_gaussian[z[0]]  # now sample the data from  N_k
    x, y = np.random.multivariate_normal(
        *normal, 1).T  # sample 1 point at a time for clarity
    plt.plot(x, y, '.', color='black')
plt.axis('equal')
plt.title('p(x)')
plt.show()

Learning!

from sklearn.mixture import GaussianMixture
# X is 2xN
gmm = GaussianMixture(n_components=3, random_state=0).fit(X.T)
# >>> gm.means_ #Dxk
assignments = gmm.predict(X.T)

Now we can do multiple things: 1) clustering

z_to_col = {0: 'red', 1: 'green', 2: 'blue'}
_ = plt.scatter(X[0, ...], X[1, ...], color=[z_to_col[y] for y in assignments])
_ = plt.scatter(*gmm.means_.T,
                  s=250,
                  marker='*',
                  c='black',
                  label='centroids')

Note that we have recovered the labels automatically!

Now we can do multiple things: 1) density estimators 2) anomaly detection

You cannot do this with k-means

from matplotlib.colors import LogNorm
fig = plt.figure(figsize=(10, 10))
# display predicted scores by the model as a contour plot
x = np.linspace(-4.0, 6.0)
y = np.linspace(-4.0, 6.0)
X1, Y1 = np.meshgrid(x, y)
XX = np.array([X1.ravel(), Y1.ravel()]).T
Z = -gmm.score_samples(XX)  # negative log likelihood
Z = Z.reshape(X1.shape)

CS = plt.contourf(
    X1, Y1, Z, cmap='hot'
)
plt.grid(False)
CB = plt.colorbar(CS, shrink=0.8, extend="both")
plt.scatter(X[0, ...],  X[1, ...])

plt.title("Negative log-likelihood predicted by a GMM")
plt.axis("tight")
plt.show()

Now we can do multiple things: 3) we can also generate (sample) data as we showed before

We just have to sample from the parameters that we have now estimated!

Motivating GMM: Weaknesses of K-Means

This part is taken from this tutorial

%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
import numpy as np

# Generate some data
from sklearn.datasets import make_blobs
X, y_true = make_blobs(n_samples=400, centers=4,
                       cluster_std=0.60, random_state=0)
X = X[:, ::-1] # flip axes for better plotting
# Plot the data with K Means Labels
from sklearn.cluster import KMeans
kmeans = KMeans(4, random_state=0)
labels = kmeans.fit(X).predict(X)
plt.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis');

from sklearn.cluster import KMeans
from scipy.spatial.distance import cdist

def plot_kmeans(kmeans, X, n_clusters=4, rseed=0, ax=None):
    labels = kmeans.fit_predict(X)

    # plot the input data
    ax = ax or plt.gca()
    ax.axis('equal')
    ax.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis', zorder=2)

    # plot the representation of the KMeans model
    centers = kmeans.cluster_centers_
    radii = [cdist(X[labels == i], [center]).max()
             for i, center in enumerate(centers)]
    for c, r in zip(centers, radii):
        ax.add_patch(plt.Circle(c, r, fc='#CCCCCC', lw=3, alpha=0.5, zorder=1))
kmeans = KMeans(n_clusters=4, random_state=0)
plot_kmeans(kmeans, X)

rng = np.random.RandomState(13)
X_stretched = np.dot(X, rng.randn(2, 2))

kmeans = KMeans(n_clusters=4, random_state=0)
plot_kmeans(kmeans, X_stretched)

from matplotlib.patches import Ellipse

def draw_ellipse(position, covariance, ax=None, **kwargs):
    """Draw an ellipse with a given position and covariance"""
    ax = ax or plt.gca()
    
    # Convert covariance to principal axes
    if covariance.shape == (2, 2):
        U, s, Vt = np.linalg.svd(covariance)
        angle = np.degrees(np.arctan2(U[1, 0], U[0, 0]))
        width, height = 2 * np.sqrt(s)
    else:
        angle = 0
        width, height = 2 * np.sqrt(covariance)
    
    # Draw the Ellipse
    for nsig in range(1, 4):
        ax.add_patch(Ellipse(position, nsig * width, nsig * height,
                             angle, **kwargs))
        
def plot_gmm(gmm, X, label=True, ax=None):
    ax = ax or plt.gca()
    labels = gmm.fit(X).predict(X)
    if label:
        ax.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis', zorder=2)
    else:
        ax.scatter(X[:, 0], X[:, 1], s=40, zorder=2)
    ax.axis('equal')
    
    w_factor = 0.2 / gmm.weights_.max()
    for pos, covar, w in zip(gmm.means_, gmm.covariances_, gmm.weights_):
        draw_ellipse(pos, covar, alpha=w * w_factor)

from sklearn.mixture import GaussianMixture
gmm = GaussianMixture(n_components=4).fit(X)
labels = gmm.predict(X)
plt.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis');

gmm = GaussianMixture(n_components=4, random_state=42)
plot_gmm(gmm, X)

gmm = GaussianMixture(n_components=4, covariance_type='full', random_state=42)
plot_gmm(gmm, X_stretched)

Covariance Type

GMM can have assumptions even on the shape of the Gaussian based on the covariance matrix.

Remember:

• $\bmu$ indicates where to place the Gaussian (acts as a translation in space);
• $\bsigma$ instead models the shape of each Gaussian

$\bsigma$ types:

• full: each component has its own general covariance matrix (more parameters)
• tied: all components share the same general covariance matrix.
• diag each component has its own diagonal covariance matrix.
• spherical: each component has its own single variance. (less params)

Note: It is common to decorrelate the data with PCA so that then you can assume a diagonal matrix (less params!)

Machine Learning

---

5. GMM as Density Estimator

Prof. Iacopo Masi - Computer Science Department, Sapienza University of Rome

GMM as Density Estimator

from sklearn.datasets import make_moons
Xmoon, ymoon = make_moons(200, noise=.05, random_state=0)
plt.scatter(Xmoon[:, 0], Xmoon[:, 1]);

gmm2 = GaussianMixture(n_components=2, covariance_type='full', random_state=0)
plot_gmm(gmm2, Xmoon)

With two components does not fit the data well

let's try with 8 and full covariance matrices!

gmm16 = GaussianMixture(n_components=8, covariance_type='full', random_state=0)
plot_gmm(gmm16, Xmoon, label=False)

GMM models the overall distribution of the data

• Here the mixture of 8 Gaussians serves not to find separated clusters of data, but rather to model the overall distribution of the input data.
• This is a generative model of the distribution, meaning that the GMM gives us the recipe to generate new random data distributed similarly to our input.
• For example, here are 400 new points drawn from this 16-component GMM fit to our original data:

We can sample from the estimated distribution

Note that sample means generate

Xnew, _ = gmm16.sample(4000)
plot_gmm(gmm16, Xmoon, label=True)
plt.scatter(Xnew[:, 0], Xnew[:, 1], color='red')

we can even sample more points and plot their log-likelihood

Xnew, _ = gmm16.sample(40000)
loglike = gmm16.score_samples(Xnew)
plt.scatter(Xnew[:, 0], Xnew[:, 1], c=loglike, cmap='jet')

How many components? (similar to K-means)

• The fact that GMM is a generative model gives us a natural means of determining the optimal number of components for a given dataset.
• A generative model is inherently a probability distribution for the dataset, and so we can simply evaluate the likelihood of the data under the model, using cross-validation to avoid over-fitting.
• Another means of correcting for over-fitting is to adjust the model likelihoods using some analytic criterion such as the Akaike information criterion (AIC) or the Bayesian information criterion (BIC).
• Scikit-Learn's GMM estimator includes built-in methods that compute both of these, so it is very easy to operate on this approach.

n_components = np.arange(1, 21)
models = [GaussianMixture(n, covariance_type='full', random_state=0).fit(Xmoon)
          for n in n_components]

plt.plot(n_components, [m.bic(Xmoon) for m in models], label='BIC')
plt.plot(n_components, [m.aic(Xmoon) for m in models], label='AIC')
plt.legend(loc='best')
plt.xlabel('n_components');

Appendix - A generative model on real data

from sklearn.datasets import load_digits
digits = load_digits()
digits.data.shape

def plot_digits(data):
    fig, ax = plt.subplots(10, 10, figsize=(8, 8),
                           subplot_kw=dict(xticks=[], yticks=[]))
    fig.subplots_adjust(hspace=0.05, wspace=0.05)
    for i, axi in enumerate(ax.flat):
        im = axi.imshow(data[i].reshape(8, 8), cmap='gray')
        im.set_clim(0, 16)
plot_digits(digits.data)

The input dimension is D-64. We can try to compress it for better speed

• We have nearly 1,800 digits in 64 dimensions, and we can build a GMM on top of these to generate more.
• GMMs can have difficulty converging in such a high dimensional space, so we will start with an invertible dimensionality reduction algorithm on the data.
• Here we will use a straightforward PCA, asking it to preserve 99% of the variance in the projected data

from sklearn.decomposition import PCA
pca = PCA(0.99, whiten=False)
data = pca.fit_transform(digits.data)
data.shape

64 dimensions to 41

• The result is 41 dimensions, a reduction of nearly 1/3 with almost no variance loss.
• Given this projected data, let's use the BIC to get a gauge for the number of GMM components we should use.

Lets' now use BIC regularization to decide the number of components

• Note by doing this we have to run the algorithm for how many components we need to try.

n_components = np.arange(1, 40, 1)
models = [GaussianMixture(n, covariance_type='diag', random_state=0)
          for n in n_components]
aics = [model.fit(data).bic(data) for model in models]
plt.plot(n_components, aics);

We will use 10 components (because we know we have 10 digits!)

We could hope that each Gaussian model the digit class!

idx = (n_components==10).nonzero()[0][0]
model = models[idx]

Now let's generate (sample) 100 new data points

data_new, y_new = model.sample(100)
data_new.shape

Now from subspace to the Input space

digits_new = pca.inverse_transform(data_new)
plot_digits(digits_new)

These digits [should] not exist in the initial dataset!

GMM Model Debug

plt.bar(*np.unique(digits.target, return_counts=True));

plt.bar(range(10),model.weights_);

plt.bar(range(model.means_[0].shape[0]),model.means_[0]);

For example, you can check if the mean is very close to a training point

plt.bar(range(model.means_[0].shape[0]),model.covariances_[0]);

You could debug covariance to spot potential singularities (though sklearn may avoid that)

Homework 3

• Fit a GMM "manually" using the labels provided by digits.target with K=10
• Fit now the same GMM without labels and debug/confront the parameters of the models
• Try to draw some conclusions regarding the difference