Training data

In real scenario, we observed a set of values, often called features, and we try to learn some model that can "explain" the observed features. Generally, we do not know the true distribution of the features. For this work we will simulate a set of single dimensional features sampled from a GMM. We will call this GMM the "true model" and in this notebook we will see different strategies to approximate the true model from the simulated training data.

Gaussian density

We can first try to model the observed data with a Gaussian density function. The Gaussian density function is defined as: $$ \DeclareMathOperator{\Norm}{\mathcal{N}} \DeclareMathOperator{\Gam}{Gam} \DeclareMathOperator{\e}{exp} p(x \mid \mu, \sigma^2) = \Norm(x \mid \mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \e\{ \frac{-(x - \mu)^2}{2\sigma^2} \} $$

• $\mu$ is the mean of the Gaussian density
• $\sigma^2$ is the variance of the Gaussian density

TODO

Try to play with the mean and the variance of the Gaussian. How each parameter change the density ?

Alternatively, the Gaussian density function can be parameterized with a precision $\lambda$ which is the inverse of the variance: $$ \DeclareMathOperator{\e}{exp} p(x | \mu, \lambda) = \Norm(x \mid \mu, \lambda^{-1}) = \frac{\sqrt{\lambda}}{\sqrt{2\pi}} \e\{ \frac{-\lambda(x - \mu)^2}{2} \} $$

• $\mu$ is the mean of the Gaussian density
• $\lambda = \frac{1}{\sigma^2}$ is the precision of the Gaussian density

This alternative parameterization will prove to be useful when dealing with Bayesian inference for the Gaussian/GMM density.

Maximum Likelihood estimation

We will first try to model our data with a simple Gaussian. We need to find the parameters (mean and variance) that fit the best the data.

$$ \begin{align} \mu &= \frac{1}{N} \sum_i x_i \\ \sigma^2 &= \frac{1}{N} \sum_i x_i^2 \end{align} $$

TODO

From the above equations, implement the maximum-likelihood solution for the Gaussian. Does the Gaussian seems to be a good model for the true density ?

Bayesian Inference

We can work within the Bayesian framework with the Gaussian density by putting a prior over the mean and precision. A common an convenient choice of prior for the Gaussian is the Normal-Gamma prior:

$$ p(\mu, \lambda \mid m_0, \kappa_0, a_0, b_0) = \Norm(\mu \mid m_0, (\kappa_0 \lambda)^{-1}) \Gam(\lambda \mid a_0, b_0) $$

where:

$$ \Gam(\lambda \mid a_0, b_0) = \frac{1}{\Gamma(a_0)} b_0^{a_0} \lambda^{a_0 - 1} \e \{ -b_0 \lambda\} $$

$m_0$, $\kappa_0$, $a_0$ and $b_0$ are called hyper-parameters. They are the parameters of the prior distribution.

TODO

Play with the parameters of the Normal-Gamma. How each parameter affects the density ?

Because the Normal-Gamma is the conjugate prior of the Normal density, the posterior distribution $p(\mu, \lambda \mid \mathbf{x})$ has a closed form solution:

$$ p(\mu, \lambda \mid \mathbf{x}) = \Norm(\mu \mid m_n, (\kappa_n \lambda)^{-1}) \Gam(\lambda \mid a_n, b_n) $$

where:

$$ \begin{align} m_n &= \frac{\kappa_0 m_0 + N \bar{x}} {\kappa_0 + N} \\ \kappa_n &= \kappa_0 + N \\ a_n &= a_0 + \frac{N}{2} \\ b_n &= b_0 + \frac{N}{2} ( s + \frac{\kappa_0 (\bar{x} - m_0)^2}{\kappa_0 + N} ) \\ \bar{x} &= \frac{1}{N} \sum_i x_i \\ s &= \frac{1}{N} \sum_i (x_i - \bar{x})^2 \end{align} $$

$N$ is the total number of point in the the training data and $m_n$, $\kappa_n$, $a_n$ and $b_n$ are the parameters of the posterior. Note that they are different from the hyper-parameters !!

TODO

Compute the posterior distribution with 1, 5, 10, 20 and 50 data points from the training data. What do you observe ?

Predictive probability

Now that we have our posterior distibution we can predict the probability of a new data point given the training data:

$$ p(x' \mid \mathbf{x}) = \int p(x' \mid \theta) p(\theta \mid \mathbf{x}) d \theta $$

For the Gaussian with Normal-Gamma prior the marginal predictive distribution is the Student's t distribution:

$$ \newcommand{\diff}{\mathop{}\!d} \DeclareMathOperator{\St}{\mathcal{St}} p(x' \mid \mathbf{x}) = \int_{-\infty}^{\infty} p(x|\mu, \lambda) p(\mu, \lambda \mid \mathbf{x}) \diff \mu \diff \lambda = \St(x' \mid \mu_n, \nu, \gamma) $$

where:

$$ \begin{align} \nu &= 2a_n \\ \gamma &= \frac{a_n \kappa_n}{b_n(\kappa_n + 1)} \\ \St(x' \mid \mu_n, \nu, \gamma) &= \frac{\Gamma(\frac{\nu}{2} + \frac{1}{2})}{\Gamma(\frac{\nu}{2})} \Big( \frac{\gamma}{\pi \nu} \Big)^{\frac{1}{2}} \Big[ 1 + \frac{\gamma (x - \mu_n)^2}{\nu} \Big]^{-\frac{\nu}{2} - \frac{1}{2}} \end{align} $$

TODO

Compute the log-likelihood of the predictive distribution derived from a posterior trained with 1, 5, 10, 20 and 50 data points. How does it compare to the log-likelihood of the Gaussian trained with maximum likelihood ?

Gaussian Mixture Model (GMM)

The Gaussian density is a very simple function, however, in most cases of interest the density we try to model has a complex shape that cannot be expressed with a simple formula. A solution is to assume that our complex density is made of $K$ Gaussian densities. This is called a Gaussian Mixture Model (GMM) and it is defined as: $$ p(x|\boldsymbol{\mu}, \boldsymbol{\lambda}, \boldsymbol{\pi}) = \sum_{k=1}^{K} \pi_k \Norm(x|\mu_k, (\lambda_k)^{-1}) $$

• $\boldsymbol{\mu}$ is the vector of $K$ means
• $\boldsymbol{\lambda}$ is the vector of $K$ precisions
• $\boldsymbol{\pi}$ is the vector of $K$ weights such that $\sum_{k=1}^K \pi_k = 1$

TODO

Observe the influence of each parameters on the density and add/remove some components. Does this model suffer of the same drawbacks of the Gaussian density ?

Maximum Likelihood estimation

The GMM parameters can be estimated with the Expectation-Maximization (EM) algorithm. The EM algorithm is an iterative algorithm that converges toward a (local) maximum of the log-likelihood of the data given the model. The EM training is as follows:

• initialize the parameters of the GMM
• iterate until convergence:
  • Expectation (E-step): compute the probability of the latent variable for each data point
  • Maximization (M-step): update the parameters from the statistics of the E-step.

TODO

• We have implemented a simple version of the EM algorithm that iterates 10 times. Change the following code to iterates until convergence of the log-likelihood. We will consider that the log-liklihood has converged if $ \log p(\mathbf{x} | \theta^{new}) - \log p(\mathbf{x} | \theta^{old}) \le 0.01$.
• Run the EM algorithm with different initialization: is the final log-likelihood always the same ? What can you conclude ? Is the log-likelihood greater than the one from the Gaussian estimated with Maximum-Likelihood ?

Bayesian GMM

We can also work the Bayesian inference for GMM by putting a prior over the GMM parameters. Let $\Theta$ be the set of parameters of the GMM: $$ \Theta = \{ \boldsymbol{\mu}, \boldsymbol{\lambda}, \boldsymbol{\pi} \} $$

The prior over the weights $\boldsymbol{\pi}$ will be a Dirichlet distribution:

$$ \DeclareMathOperator{\Dir}{Dir} p(\boldsymbol{\pi}) = \Dir(\boldsymbol{\pi} \mid \boldsymbol{\alpha}_0) $$

and the prior of the mean and precision of the component $k$ of the mixture will be a Normal-Gamma distribution:

$$ p(\mu_k, \lambda_k) = \Norm(\mu_k \mid m_0, (\kappa_0 \lambda_k)^{-1}) \Gam(\lambda_k \mid a_0, b_0) $$

The joint distribution of the data, the latent variables and the parameters can be written as:

$$ \begin{align} p(\mathbf{x}, \mathbf{z}, \Theta) &= p(\mathbf{x}, \mathbf{z} \mid \Theta)p(\Theta) \\ &= \Bigg[ \prod_{i=0}^{N} p(x_i \mid \boldsymbol{\mu}, \boldsymbol{\lambda}, \boldsymbol{\pi}) \Bigg] \Bigg[ \Dir(\boldsymbol{\pi} \mid \boldsymbol{\alpha}_0) \prod_{k=0}^{K} \Norm(\mu_k \mid m_0, (\kappa_0 \lambda_k)^{-1}) \Gam(\lambda_k \mid a_0, b_0) \Bigg] \end{align} $$

Gibbs Sampling:

The EM algorithm cannot be applied directly to train the Bayesian GMM. Instead, we will use Gibbs Sampling (GS) to learn the posterior distribution of the parameters. Gibbs Sampling is a simple way to sample values from a complex distribution. Let's say that we want to sample values for the random variable $A$ and $B$ conditioned on $C$. Unfortunately, $p(A, B \mid C)$ is too complex to be sampled from directly. Alternatively we can sample in turn $a$ from $p(A | B = b, C)$ and $b$ from $p(B | A=a, C)$. It can be proven that if we keep sampling long enough, the set of $a$ and $b$ will be distributed according to $p(A, B | C)$. This is the Gibbs Sampling algorithm.

NOTE: Gibbs Sampling is not always applicable as we cannot always sample from the conditional distributions (it may be intractable).

In the following we just show the 3 conditional distributions we need for the Gibbs Sampling of the Bayesian GMM.

latent variable $z_i$

$$ \begin{align} p(z_i \mid \mathbf{x}, \Theta) &= p(z_i \mid x_i, \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\lambda}) \\ p(z_i = k \mid x_i, \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\lambda}) &= \frac{\pi_k \Norm(x_i \mid \mu_k, \lambda_k)} {\sum_{j=0}^{K} \pi_j \Norm(x_i \mid \mu_j, \lambda_j)} \end{align} $$

mean and variance

We define the set of all data point $x_i$ that are assigned to the component $k$ of the mixture as follows: $$ \mathbf{x}_{(k)} = \{ x_i : z_i = k, \forall i \in \{1,... , N \} \} $$ and similarly for the latent variables $\mathbf{z}$: $$ \mathbf{z}_{(k)} = \{ z_i : z_i = k, \forall i \in \{1,... , N \} \} $$

$$ \begin{align} p(\mu_k, \lambda_k \mid \mathbf{x}, \mathbf{z}, \Theta_{\smallsetminus \{ \mu_k, \lambda_k \} } ) &= p(\mu_k, \lambda_k \mid \mathbf{x}_{(k)}, \mathbf{z}_{(k)}, \Theta_{\smallsetminus \{ \mu_k, \lambda_k \} } ) \\ &= \Norm(\mu_k \mid m_{n,k}, (\kappa_{n,k} \lambda_k)^{-1}) \Gam(\lambda_k \mid a_{n,k}, b_{n,k}) \end{align} $$

where:

$$ \begin{align} m_{n,k} &= \frac{\kappa_0 m_0 + N_k \bar{x}_k} {\kappa_0 + N_k} \\ \kappa_{n,k} &= \kappa_0 + N_k \\ a_{n,k} &= a_0 + \frac{N_k}{2} \\ b_{n,k} &= b_0 + \frac{N_k}{2} ( s + \frac{\kappa_0 (\bar{x}_k - m_0)^2}{\kappa_0 + N_k} ) \\ N_k &= \left\vert \mathbf{x}_{(k)} \right\vert \\ \bar{x}_k &= \frac{1}{N_k} \sum_{\forall x \in \mathbf{x}_{(k)}} x \\ s_n &= \frac{1}{N} \sum_{\forall x \in \mathbf{x}_{(k)}} (x_i - \bar{x})^2 \end{align} $$

NOTE: these equations are very similar to the Bayesian Gaussian estimate. However, it remains some difference

weights

$$ \begin{align} p( \boldsymbol{\pi} \mid \mathbf{x}, \mathbf{z}, \Theta_{\smallsetminus \{ \boldsymbol{\pi} \} } ) &= p( \boldsymbol{\pi} \mid \mathbf{z}) \\ &= \Dir(\boldsymbol{\pi} \mid \boldsymbol{\alpha}) \end{align} $$

where: $$ \alpha_{n,k} = \alpha_{0,k} + N_k \; ; \; \forall \, k = 1\dots K $$

TODO

• Compare the following implementation of the Gibbs Sampling to the EM algorithm. What are their differences/similarities ?
• Run the Gibbs Sampling algorithm for 1, 2, 5, 10, 20, ... iterations. How many iterations do we need to get a reasonable estimate of the true density ?
• Change the hyper-parameters and run the Gibbs Sampling algorithm for 3 iterations. Are the initial sampled GMM very similar ? Run more iterations, does the Gibbs Sampling converges to the same solution regardless of the initialization ?