# Gaussian Mixture Models Implementation


```python
import matplotlib.pyplot as plt
import numpy as np
import random
from sklearn.cluster import KMeans
import scipy.stats as stats
import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning) 
```

Since we do not have any real data, we generate it randomly from normal distributions using the below fuction.


```python
def generate_random_sample(num_gaussians, num_observations):
    
    samples = []
    samples_tup = []
    parameters = []
    size = int(num_observations/num_gaussians)
    for x in range(num_gaussians):
        mean, sd = random.randint(-30, 30), random.randint(1, 3)
        samples.append(np.random.normal(loc=mean, scale=sd, size=size))
        parameters.append((mean, sd))
        samples_tup.append((mean, samples[x]))
        
    print('Parameters of the random distributions: ' +str(parameters))
    samples_array = np.array([x for sample_x in samples for x in sample_x])
    
    return {'parameters': parameters, 'X': samples_array.reshape((len(samples_array), 1)), 'samples': samples, 'samp_tup': samples_tup}
```


```python
def gaussian(X, mu, cov):
    n = X.shape[1]
    diff = (X - mu).T
    return np.diagonal(1 / ((2 * np.pi) ** (n / 2) * np.linalg.det(cov) ** 0.5) * np.exp(-0.5 * np.dot(np.dot(diff.T, np.linalg.inv(cov)), diff))).reshape(-1, 1)
```

Let us now implement the Gaussian density function. Ahem... I know you can use numpy functions for that, but I believe it is actually interesting to see how things work internally. Our aim is to create a function that implements this:

\begin{equation}
\large
p(\mathbf x | \mathbf\mu, \mathbf\Sigma) = \frac 1 {({2\pi})^{n/2}|\Sigma|^{1/2}}\exp\left(-\frac 1 2 (\mathbf x -\mathbf\mu)^T\mathbf\Sigma^{-1}(\mathbf x -\mathbf\mu)\right)
\end{equation}


__Step 1__

This is the initialization step of the GMM. At this point, we must initialise our parameters $\pi_k$, $\mu_k$, and $\Sigma_k$. In this case, we are going to use the results of KMeans as an initial value for $\mu_k$, set $\pi_k$ to one over the number of clusters and $\Sigma_k$ to the identity matrix. We could also use random numbers for everything, but using a sensible initialisation procedure will help the algorithm achieve better results.


```python
def initialize_clusters(X, n_clusters):
    clusters = []
    idx = np.arange(X.shape[0])
    
    # We use the KMeans centroids to initialise the GMM
    
    kmeans = KMeans(n_clusters, n_init = 100).fit(X)
    
    cluster_id = kmeans.labels_
    
    mu_k = kmeans.cluster_centers_
        
    for i in range(n_clusters):
        clusters.append({
            'pi_k': 1.0 / n_clusters,
            'mu_k': mu_k[i],
            'cov_k': np.identity(X.shape[1], dtype=np.float64)
        })
        
    return clusters, cluster_id
```

**Step 2 (Expectation step)**

We should now calculate $\gamma(z_{nk})$. We can achieve this by means of the following expression:

\begin{equation}
\large
\gamma{(z_{nk})}=\frac {\pi_k\mathcal N(\mathbf x_n| \mathbf\mu_k, \mathbf\Sigma_k)}{\sum_{j=1}^K\pi_j\mathcal N(\mathbf x_n| \mathbf\mu_j, \mathbf\Sigma_j)}
\end{equation}

For convenience, we just calculate the denominator as a sum over all terms in the numerator, and then assign it to a variable named __totals__


```python
def expectation_step(X, clusters):
    totals = np.zeros((X.shape[0], 1), dtype=np.float64)
    
    for cluster in clusters:
        pi_k = cluster['pi_k']
        mu_k = cluster['mu_k']
        cov_k = cluster['cov_k']
        
        gamma_nk = (pi_k * gaussian(X, mu_k, cov_k)).astype(np.float64)
        
        for i in range(X.shape[0]):
            totals[i] += gamma_nk[i]
        
        cluster['gamma_nk'] = gamma_nk
        cluster['totals'] = totals
        
    
    for cluster in clusters:
        cluster['gamma_nk'] /= cluster['totals']
```

**Step 3 (Maximization step):**

Let us now implement the maximization step. Since $\gamma(z_{nk})$ is common to the expressions for $\pi_k$, $\mu_k$ and $\Sigma_k$, we can simply define:

\begin{equation}
\large
N_k=\sum_{n=1}^N\gamma({z_{nk}})
\end{equation}

And then we can calculate the revised parameters by using:

\begin{equation}
\large
\pi_k^*=\frac {N_k} N
\end{equation}


\begin{equation}
\large
\mu_k^*=\frac 1 {N_k} \sum_{n=1}^N\gamma({z_{nk}})\mathbf x_n
\end{equation}


\begin{equation}
\large
\Sigma_k^*=\frac 1 {N_k} \sum_{n=1}^N\gamma({z_{nk}})(\mathbf x_n-\mathbf\mu_k)(\mathbf x_n-\mathbf\mu_k)^T
\end{equation}

Note: To calculate the covariance, we define an auxiliary variable __diff__ that contains $(x_n-\mu_k)^T$.


```python
def maximization_step(X, clusters):
    N = float(X.shape[0])
    Q = []
    for cluster in clusters:
        gamma_nk = cluster['gamma_nk']
        cov_k = np.zeros((X.shape[1], X.shape[1]))
        
        N_k = np.sum(gamma_nk, axis=0) 
        pi_k = N_k / N
        mu_k = np.sum(gamma_nk * X, axis=0) / N_k
        
        for j in range(X.shape[0]):
            diff = (X[j] - mu_k).reshape(-1, 1)
            cov_k += gamma_nk[j] * np.dot(diff, diff.T)
            
        cov_k /= N_k
        
        cluster['pi_k'] = pi_k
        cluster['mu_k'] = mu_k
        cluster['cov_k'] = cov_k
        
        Q_K = np.sum(np.multiply(gamma_nk, np.log(pi_k)))

```

Let us now determine the log-likelihood of the model. It is given by:

\begin{equation}
\large
\ln p(\mathbf X)=\sum_{n=1}^N\ln\sum_{k=1}^K\pi_k\mathcal N(\mathbf x_n|\mu_k,\Sigma_k)
\end{equation}

However, the second summation has already been calculated in the __expectation_step__ function and is available in the __totals__ variable. So we just make use of it.


```python
def get_likelihood(X, clusters):
    sample_likelihoods = np.log(np.array([cluster['totals'] for cluster in clusters]))
    return np.sum(sample_likelihoods), sample_likelihoods
```

Finally, let's put everything together! First, we are going to initialise the parameters by using the __initialise_clusters__ function, and then perform several expectation-maximization steps. In this case, we set the number of iterations of the training procedure to a fixed __n_epochs__ number. I have done this on purpose to generate graphs of the log-likelihood later.


```python
def train_gmm(X, n_clusters, n_epochs):
    
    clusters, cluster_id = initialize_clusters(X, n_clusters)
    likelihoods = np.zeros((n_epochs, ))
    scores = np.zeros((X.shape[0], n_clusters))
    history = []

    for i in range(n_epochs):
        clusters_snapshot = []
        
        # This is just for our later use in the graphs
        for cluster in clusters:
            clusters_snapshot.append({
                'mu_k': cluster['mu_k'].copy(),
                'cov_k': cluster['cov_k'].copy()
            })
            
        history.append(clusters_snapshot)
      
        expectation_step(X, clusters)
        maximization_step(X, clusters)

        likelihood, sample_likelihoods = get_likelihood(X, clusters)
        likelihoods[i] = likelihood
        
    for i, cluster in enumerate(clusters):
        scores[:, i] = np.log(cluster['gamma_nk']).reshape(-1)
    return clusters, likelihoods, scores, sample_likelihoods, history, cluster_id
```

# Find the best number of Gaussians

The below function estimates the fit error of each gaussian using the Kolmogorov–Smirnov test (stats.ks_2samp). It returns the Kolmogorov–Smirnov test values and the parameters of each Gaussian (mean and standard deviation)


```python
def compare(sol, clusters, n_clusters, num_observations, cluster_Xs):
    
    bins = 50
    fitness = []
    parameters = []
    
    for i in range(n_clusters):
        
        X = cluster_Xs[cluster_Xs[:,1] == i][:,0]    
    size = int(num_observations/n_clusters)
    for n in range(n_clusters):
        mean = float(clusters[n]['mu_k'])
        standard_deviation = np.sqrt(float(clusters[n]['cov_k']))
           
        X = cluster_Xs[cluster_Xs[:,1] == n][:,0]
        rvs1 = stats.norm.rvs(size=size, loc=mean, scale=standard_deviation)
        fitness.append(stats.ks_2samp(rvs1, X)[0])
        x_values = np.arange(-50, 50, 0.01)
        y_values = stats.norm(mean, standard_deviation)
                
        parameters.append((mean, standard_deviation))
        
    return np.mean(fitness), parameters

```


```python
def show_best(X,n_clusters):
    left, right = min(X), max(X)
    n_epochs = 2
    clusters, likelihoods, scores, sample_likelihoods, history, cluster_Xs = train_gmm(X, n_clusters, n_epochs)
    cluster_Xs = cluster_Xs.reshape((len(cluster_Xs), 1))
    cluster_Xs =  np.column_stack((X, cluster_Xs))
    
    bins = 20
 
    plt.figure(figsize=(15, 6), dpi=80)
    for i in range(n_clusters):
        
        X = cluster_Xs[cluster_Xs[:,1] == i][:,0] 
        plt.hist(X, density=True, bins = bins, alpha = 0.4)[0]
    
    
    size = int(num_observations/n_clusters)
    for n in range(n_clusters):
        mean = float(clusters[n]['mu_k'])
        standard_deviation = np.sqrt(float(clusters[n]['cov_k']))
           
        X = cluster_Xs[cluster_Xs[:,1] == n][:,0]
        rvs1 = stats.norm.rvs(size=size, loc=mean, scale=standard_deviation)
        x_values = np.arange(left*1.2, right*1.2, 0.01)
        y_values = stats.norm(mean, standard_deviation)
        
        plt.plot(x_values, y_values.pdf(x_values), linewidth=4)
        
    
    plt.show()
    
    
```


```python
def gaussian_mixture(X, max_gaussians = 10):
    num_observations = len(X)
    fitness = []
    parameters = []
    for i in range(max_gaussians):
        n_clusters = i+1
        n_epochs = 2
        clusters, likelihoods, scores, sample_likelihoods, history, cluster_Xs = train_gmm(X, n_clusters, n_epochs)
        cluster_Xs = cluster_Xs.reshape((len(cluster_Xs), 1))
        cluster_Xs =  np.column_stack((X, cluster_Xs))
        fit, param = compare(sol, clusters, n_clusters, num_observations, cluster_Xs)
        fitness.append(fit)
        parameters.append(param)
    best_n_clusters = fitness.index(min(fitness))+1
    best_parameters = parameters[fitness.index(min(fitness))]
    print('The best Nº of Guassian fits is ' + str(fitness.index(min(fitness))+1))
    print('The paramters are: ' + str(parameters[fitness.index(min(fitness))]))
    print('The fit looks like: ')
    show_best(X,best_n_clusters)
    
    return parameters[fitness.index(min(fitness))]
```


```python
num_gaussians = 5
num_observations = 1000
sol = generate_random_sample(num_gaussians,num_observations)
X = sol['X']

parameters = gaussian_mixture(X)
```

e03ad8d767e23c2609a7cdedfb4d37c6d42d292e
ipynb
Jupyter Notebook
hyperactive/HyperActive_Multiple_Gaussian.ipynb
nevinadalal/Hyperactive
3232ffeda70c5d4853b9e71aaf5d1e761c0db9c2
MIT
hyperactive/HyperActive_Multiple_Gaussian.ipynb
nevinadalal/Hyperactive
3232ffeda70c5d4853b9e71aaf5d1e761c0db9c2
MIT
hyperactive/HyperActive_Multiple_Gaussian.ipynb
nevinadalal/Hyperactive
3232ffeda70c5d4853b9e71aaf5d1e761c0db9c2
MIT
Qwen/Qwen-72B
1. YES
2. YES