11-DensityEstimationWithGaussianMixtureModels.Rmd

# Density Estimation with Gaussian Mixture Models

In addition to regression and dimensionality reduction, **density estimation** is a core pillar of machine learning.  Its goal is to represent data compactly by estimating the **underlying probability density function** that generated it.



Instead of storing all data points, density estimation models the data using a **parametric family of distributions**, such as a Gaussian.  

For example, we can represent a dataset using its mean and variance, found via maximum likelihood (MLE) or maximum a posteriori (MAP) estimation.  However, a single Gaussian distribution often provides a poor fit for complex or multimodal data (data with multiple clusters).  To handle this, we use **mixture models**.



<div class="definition">
A **mixture model** represents a probability density as a **convex combination** of simpler component distributions:
\[
p(\mathbf{x}) = \sum_{k=1}^{K} \pi_k p_k(\mathbf{x}),
\]
subject to:
\[
0 \leq \pi_k \leq 1, \quad \sum_{k=1}^{K} \pi_k = 1,
\]
where

- \( p_k(\mathbf{x}) \) are individual component distributions (e.g., Gaussian, Bernoulli),  and    
- \( \pi_k \): mixture weights representing the contribution of each component  
</div>

Mixture models can capture multimodal structures in data and are thus more expressive than single distributions.


---

## Gaussian Mixture Models (GMMs)




<div class="definition">
A **Gaussian Mixture Model (GMM)** combines multiple Gaussian distributions:
\[
p(\mathbf{x} | \boldsymbol{\boldsymbol{\theta}}) = \sum_{k=1}^{K} \pi_k \, \mathcal{N}(\mathbf{x} | \boldsymbol{\boldsymbol{\mu}}_k, \boldsymbol{\boldsymbol{\Sigma}}_k),
\]
where:

- \( \boldsymbol{\boldsymbol{\theta}} = \{\pi_k, \boldsymbol{\boldsymbol{\mu}}_k, \boldsymbol{\boldsymbol{\Sigma}}_k : k = 1, \dots, K\} \),  
- each \( \mathcal{N}(\mathbf{x} | \boldsymbol{\boldsymbol{\mu}}_k, \boldsymbol{\boldsymbol{\Sigma}}_k) \) is a Gaussian component, and    
- the mixture weights \( \pi_k \) satisfy \( \sum_{k=1}^{K} \pi_k = 1 \).
</div>


This convex combination of Gaussians provides far greater flexibility for modeling complex or clustered data. 



<div class="example">
A GMM represents a probability distribution as a weighted sum of Gaussian (normal) distributions. 


Let \( x \in \mathbb{R} \) be a continuous random variable.  A **2-component GMM** is defined as:
\[
p(x) = \pi_1 \, \mathcal{N}(x \mid \mu_1, \sigma_1^2)
     + \pi_2 \, \mathcal{N}(x \mid \mu_2, \sigma_2^2),
\]
where:

- \( \pi_1, \pi_2 \ge 0 \) are **mixing coefficients**  
- \( \pi_1 + \pi_2 = 1 \)  
- \( \mathcal{N}(x \mid \mu_k, \sigma_k^2) \) is a Gaussian pdf  


Let:
\[
\pi_1 = 0.4, \quad \pi_2 = 0.6
\]
\[
\mu_1 = -2, \quad \sigma_1^2 = 1
\]
\[
\mu_2 = 3, \quad \sigma_2^2 = 2
\]
Then the model becomes:
\[
p(x) = 0.4 \, \mathcal{N}(x \mid -2, 1)
     + 0.6 \, \mathcal{N}(x \mid 3, 2)
\]

So,

- With probability 0.4, a data point is generated from a Gaussian centered at −2
- With probability 0.6, a data point is generated from a Gaussian centered at 3
- The overall distribution is **bimodal**, with two peaks.


Given an observed value \( x \), the probability that it came from component \( k \) is:

\[
\gamma_k(x) = p(z = k \mid x)
= \frac{\pi_k \mathcal{N}(x \mid \mu_k, \sigma_k^2)}
{\sum_{j=1}^{2} \pi_j \mathcal{N}(x \mid \mu_j, \sigma_j^2)}
\]
These are called **responsibilities** and are used in the **Expectation–Maximization (EM)** algorithm.

**Example Calculation**

Suppose we observe \( x = 0 \).

- Component 1 likelihood:
\[
\mathcal{N}(0 \mid -2, 1) \approx \frac{0.1353}{2.5066}
\approx 0.054.
\]

- Component 2 likelihood:
\[
\mathcal{N}(0 \mid 3, 2)
\approx \frac{0.1054}{3.5449}
\approx 0.0297.
\]
After weighting by \( \pi_1 \) and \( \pi_2 \), we compute \( \gamma_1(0) \) and \( \gamma_2(0) \) to determine which Gaussian most likely generated the data point.
\[\gamma_1(0) = \dfrac{0.4(0.054)}{0.4(0.054) + 0.6(0.0297)} = 0.548 \quad \gamma_2(0) = \dfrac{0.6(0.0297)}{0.4(0.054) + 0.6(0.0297)} = 0.452.
\]
Therefore, it is more likely that the point came from Gaussian 1.
</div>



Unlike linear regression or PCA, GMMs do not have a closed-form MLE solution.  Instead, parameter estimation is achieved iteratively, most commonly using the Expectation-Maximization (EM) algorithm.


---


### Exercises {.unnumbered .unlisted}











## Parameter Learning via Maximum Likelihood

We are given a dataset \( \mathcal{X} = \{\textbf{x}_1, \dots, \textbf{x}_N\} \), where each data point \( \textbf{x}_n \) is drawn i.i.d. from an unknown distribution \( p(\textbf{x}) \).  Our goal is to approximate this distribution using a Gaussian mixture model with \( K \) components, each defined by:
\[
\boldsymbol{\theta} = \{\pi_k, \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k : k = 1, \dots, K\}
\]
where:

- \( \pi_k \) is the mixture weights for the $k^{th}$ distribution,  
- \( \boldsymbol{\mu}_k \) is the mean for the $k^{th}$ distribution, and  
- \( \boldsymbol{\Sigma}_k \) is the covariance associated with the $k^{th}$ distribution.

---

### Likelihood and Log-Likelihood

The likelihood of the dataset is:
\[
p(\mathcal{X} | \boldsymbol{\theta}) = \prod_{n=1}^{N} p(\textbf{x}_n | \boldsymbol{\theta}), \quad 
p(\textbf{x}_n | \boldsymbol{\theta}) = \sum_{k=1}^{K} \pi_k \mathcal{N}(\textbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k).
\]

The log-likelihood becomes:
\[
L(\boldsymbol{\theta}) = \sum_{n=1}^{N} \log \sum_{k=1}^{K} \pi_k \mathcal{N}(\textbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k).
\]

Maximizing this log-likelihood yields the maximum likelihood estimate (MLE) \( \boldsymbol{\theta}_{ML} \).  However, because of the log-sum term, there is no closed-form solution, so we rely on iterative optimization, leading to the **Expectation-Maximization (EM) algorithm**.

---

### Responsibilities



<div class="definition">
The **responsibility** \( r_{nk} \) as the probability that mixture component \( k \) generated data point \( \textbf{x}_n \):
\[
r_{nk} = \frac{\pi_k \mathcal{N}(\textbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)}{\sum_{j=1}^{K} \pi_j \mathcal{N}(\textbf{x}_n | \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j)}.
\]
</div>

Each row \( r_n = [r_{n1}, \dots, r_{nK}] \) forms a probability vector where \( \sum_k r_{nk} = 1 \).

---

### Updating Parameters

The EM algorithm alternates between computing responsibilities (E-step) and updating parameters (M-step):

**1. Update Means**
\[
\boldsymbol{\mu}_k^{new} = \frac{\sum_{n=1}^{N} r_{nk} \textbf{x}_n}{\sum_{n=1}^{N} r_{nk}} = \frac{1}{N_k} \sum_{n=1}^{N} r_{nk} \textbf{x}_n,
\]
where \( N_k = \sum_{n=1}^{N} r_{nk} \) is the effective number of points assigned to component \( k \).

**2. Update Covariances**
\[
\boldsymbol{\Sigma}_k^{new} = \frac{1}{N_k} \sum_{n=1}^{N} r_{nk} (\textbf{x}_n - \boldsymbol{\mu}_k)(\textbf{x}_n - \boldsymbol{\mu}_k)^{\top}.
\]
This represents a **responsibility-weighted covariance**.

**3. Update Mixture Weights**
\[
\pi_k^{new} = \frac{N_k}{N}.
\]
This ensures that the mixture weights sum to 1, reflecting the proportion of data explained by each component.

---

### Interpretation

At each update step, three things happen in the algorithm:

- It re-estimates model parameters based on the soft assignment of data points (responsibilities).  
- It moves the Gaussian components toward regions of high data density.  
- It increases the log-likelihood, converging to a **local optimum**.

Together, these iterative updates form the **Expectation-Maximization (EM)** algorithm for GMMs, which alternates between inferring hidden assignments (E-step) and maximizing parameters (M-step) until convergence.






---


### Exercises {.unnumbered .unlisted}





<div class="exercise">
 If we were to consider a single Gaussian as the desired density, the sum over $k$ in \[\sum_{n=1}^N \log \sum_{k=1}^K \pi_kN(\mathbf{x}_n|\mathbf{\mu}_k,\mathbf{\Sigma}_k)\] vanishes.  So, $\log p(X|\mathbf{\theta})$ equals what?
 
<div style="text-align: right;">
[Solution]( )
</div> 
</div>




<div class="exercise">
Given $p(X|\mathbf{\theta}) = \prod_{n=1}^N p(\mathbf{x}_n|\mathbf{\theta})$ and $p(\mathbf{x}_n|\mathbf{\theta}) = \sum_{k=1}^K \pi_k N(\mathbf{x}_n|\mathbf{\mu}_k, \mathbf{\Sigma}_k)$, derive the log likelihood and the necessary conditions for a local optimum.

<div style="text-align: right;">
[Solution]( )
</div> 
</div>




<div class="exercise">
 Prove the update of the mean parameters $\mathbf{\mu}_k$ of the GMM is given by \[\mathbf{\mu}_k^{new} = \dfrac{\sum_{n=1}^N r_{nk}\mathbf{x}_n}{\sum_{n=1}^N r_{nk}},\] where $r_{nk}$ is the responsibility.
 
<div style="text-align: right;">
[Solution]( )
</div> 
</div>




<div class="exercise">
Prove the update of the covariance parameters $\pi_k$ of the GMM is given by \[\pi_k^{new} = \dfrac{N_k}{N},\] where $N$ is the number of data points and $N_k$ the column sum of the responsibility matrix.

<div style="text-align: right;">
[Solution]( )
</div> 
</div>




<div class="exercise">
Prove the update of the mixture weights $\mathbf{\Sigma}_k$ of the GMM is given by \[\mathbf{\Sigma}_k^{new} = \dfrac{1}{N_k}\sum_{n=1}^N r_{nk}(\mathbf{x}_n-\mathbf{\mu}_k)(\mathbf{x}_n - \mathbf{\mu}_k)^T,\] where $r_{nk}$ is the responsibility and $N_k$ the column sum of the responsibility matrix.

<div style="text-align: right;">
[Solution]( )
</div> 
</div>




## Expectation-Maximization (EM) Algorithm

The **Expectation-Maximization (EM) algorithm**, introduced by *Dempster et al. (1977)*, is an iterative method for estimating parameters in models with **latent variables**, such as Gaussian mixture models.  In GMMs, parameters \(\mu_k, \boldsymbol{\Sigma}_k, \pi_k\) cannot be solved in closed form because the responsibilities \(r_{nk}\) depend on them in a complex, nonlinear way.  Therefore, the EM algorithm alternates between estimating responsibilities and updating parameters until convergence.

---

### Algorithm Steps

1. **Initialize parameters**    
   Select initial values for each parameter we want to find:
   \[
   \{\pi_k, \mu_k, \boldsymbol{\Sigma}_k\}_{k=1}^{K}.
   \]

2. **E-step (Expectation)**  
   Compute the responsibility of component \(k\) for data point \(\mathbf{x}_n\):
   \[
   r_{nk} = 
   \frac{\pi_k \, \mathcal{N}(\mathbf{x}_n | \mu_k, \boldsymbol{\Sigma}_k)}
   {\sum_{j=1}^{K} \pi_j \, \mathcal{N}(\mathbf{x}_n | \mu_j, \boldsymbol{\Sigma}_j)}.
   \]

3. **M-step (Maximization)**  
   Update the parameters using the new responsibilities:
   \[
   N_k = \sum_{n=1}^{N} r_{nk}.
   \]

   **Update means:**
   \[
   \mu_k = \frac{1}{N_k} \sum_{n=1}^{N} r_{nk} \mathbf{x}_n.
   \]

   **Update covariances:**
   \[
   \boldsymbol{\Sigma}_k = \frac{1}{N_k} \sum_{n=1}^{N} r_{nk}(\mathbf{x}_n - \mu_k)(\mathbf{x}_n - \mu_k)^{\top}.
   \]

   **Update mixture weights:**
   \[
   \pi_k = \frac{N_k}{N}.
   \]

Each EM iteration increases the log-likelihood, and convergence can be monitored by checking either the change in log-likelihood or parameter stability.

---

### Example: GMM Fit

After several iterations (typically few), the EM algorithm converges to a stable mixture model.  For instance, a fitted GMM might be:
\[
p(\mathbf{x}) = 0.29\,\mathcal{N}(\mathbf{x}|-2.75, 0.06) 
      + 0.28\,\mathcal{N}(\mathbf{x}|-0.5, 0.25) 
      + 0.43\,\mathcal{N}(\mathbf{x}|3.64, 1.63).
\]
Plots of negative log-likelihood across iterations demonstrate steady improvement until convergence.







---


### Exercises {.unnumbered .unlisted}






<div class="exercise">
Use Excel to analyze the GMM for 2 Gaussians and data $(-3, -1,0, 1, 3, 4)$ 

<div style="text-align: right;">
[Solution]( )
</div> 
</div>




<div class="exercise">
Use Excel to analyze the GMM for 3 Gaussians and data $(-3, -1,0, 1, 3, 4)$ 
    
<div style="text-align: right;">
[Solution]( )
</div> 
</div>




<div class="exercise">
Use Excel to analyze the GMM for 2 Gaussians and data $(-4,-3, -1,0, 1, 3, 4)$

<div style="text-align: right;">
[Solution]( )
</div> 
</div>




<div class="exercise">
Use Excel to analyze the GMM for 3 Gaussians and data $(-4,-3, -1,0, 1, 3, 4)$ 
    
<div style="text-align: right;">
[Solution]( )
</div> 
</div>





<div class="exercise">
Use Excel to analyze the GMM for 3 Gaussians and data $(-4,-3, -2.5, 0, 1, 3, 4)$

<div style="text-align: right;">
[Solution]( )
</div> 
</div>




<div class="exercise">
Use Excel to analyze the GMM for 3 Gaussians and data $(-4,-3, -3, 0, 1, 3, 4)$
    
<div style="text-align: right;">
[Solution]( )
</div> 
</div>




---

##  Latent-Variable Perspective

A GMM can be interpreted as a latent-variable model where each data point \(\mathbf{x}_n\) is associated with a hidden variable \(z_n \in \{0,1\}\) indicating if the $n^{th}$ mixture component generated it.


### Generative Model

Each component \(k\) is selected according to the mixture probabilities \(\pi = [\pi_1, \dots, \pi_K]\), and data is sampled as:

1. Sample \(z \sim \text{Categorical}(\pi)\)  
2. Sample \(\mathbf{x} \sim \mathcal{N}(\mu_k, \boldsymbol{\Sigma}_k)\) given \(z_k = 1.\)  
Thus, the joint distribution is:
\[
p(\mathbf{x}, z_k = 1) = \pi_k \, \mathcal{N}(\mathbf{x} | \mu_k, \boldsymbol{\Sigma}_k).
\]

---

###   Likelihood

The marginal likelihood is obtained by summing over latent states:
\[
p(\mathbf{x} | \boldsymbol{\theta}) = \sum_{k=1}^{K} \pi_k \, \mathcal{N}(\mathbf{x} | \mu_k, \boldsymbol{\Sigma}_k)
\]

For a dataset \(\mathcal{X} = \{\mathbf{x}_1, \dots, \mathbf{x}_N\}\), the total likelihood is:
\[
p(\mathcal{X} | \boldsymbol{\theta}) = \prod_{n=1}^{N} \sum_{k=1}^{K} \pi_k \, \mathcal{N}(\mathbf{x}_n | \mu_k, \boldsymbol{\Sigma}_k).
\]
This formulation is identical to the GMM likelihood derived earlier.

---

###   Posterior Distribution

Using Bayes’ theorem, the posterior probability (responsibility) that component \(k\) generated \(\mathbf{x}_n\) is:
\[
p(z_k = 1 | \mathbf{x}_n) = 
\frac{\pi_k \, \mathcal{N}(\mathbf{x}_n | \mu_k, \boldsymbol{\Sigma}_k)}
{\sum_{j=1}^{K} \pi_j \, \mathcal{N}(\mathbf{x}_n | \mu_j, \boldsymbol{\Sigma}_j)} 
= r_{nk}.
\]
Thus, the responsibilities from the EM algorithm have a probabilistic interpretation as posterior probabilities.

---

###  Extension to Full Dataset

Each data point \(\mathbf{x}_n\) has its own latent variable \(z_n\), forming a set of hidden assignments:
\[
z_n = [z_{n1}, \dots, z_{nK}]^{\top}.
\]
The same prior \(\pi\) applies to all, and the joint conditional distribution factorizes as:
\[
p(\mathbf{x}_1, \dots, \mathbf{x}_N | z_1, \dots, z_N) = \prod_{n=1}^{N} p(\mathbf{x}_n | z_n).
\]
Responsibilities \(r_{nk}\) again represent \(p(z_{nk} = 1 | \mathbf{x}_n)\).

---

###  EM Algorithm Revisited

From the latent-variable perspective, EM can be derived as maximizing the expected complete-data log-likelihood:
\[
Q(\boldsymbol{\theta} | \boldsymbol{\theta}^{(t)}) = \mathbb{E}_{z | \mathbf{x}, \boldsymbol{\theta}^{(t)}}[\log p(\mathbf{x}, z | \boldsymbol{\theta})].
\]

- **E-step:** Compute the expected value of the log-likelihood under the posterior \(p(z | \mathbf{x}, \boldsymbol{\theta}^{(t)})\).  
- **M-step:** Maximize this expectation with respect to \(\boldsymbol{\theta}\).

Each iteration increases the log-likelihood but may converge to a local maximum, depending on initialization.

---




### Exercises {.unnumbered .unlisted}


Put some exercises here.