Gaussian_mixture_model.py

from sklearn.datasets import fetch_openml
from sklearn.model_selection import train_test_split
from sklearn.decomposition import PCA
import numpy as np
import pandas as pd
import pickle
import os
import matplotlib.pyplot as plt
CLASS_SIZE=40
MIXTURE_COUNT=10
class multivariate_gaussian:
    def __init__(self,mean,cov):
        self.mean = mean
        self.cov=cov
        self.n=len(mean)
    def logpdf(self,x):
        if isinstance(x, np.ndarray):
            jitter = 1e-6 * np.eye(self.n)
            cov_reg = self.cov + jitter
            try:
                cov_inv = np.linalg.inv(cov_reg)
            except np.linalg.LinAlgError:
                return -np.inf
            sign, log_det_cov = np.linalg.slogdet(cov_reg)
            if sign <= 0:
                return -np.inf
            diff = x - self.mean
            numenartor = diff.T @ cov_inv @ diff
            log_prob = -0.5 * numenartor - 0.5 * log_det_cov - (self.n / 2) * np.log(2 * np.pi)
            return log_prob
        else: 
            return 0
def expectation(X,gm,alpha):
    m=X.shape[0]
    gamma=np.zeros((m,MIXTURE_COUNT))
    for i in range(m): 
        for k in range(MIXTURE_COUNT):
            # p(xi,ki|params)
            LP=np.log(alpha[k])+gm[k].logpdf(X[i,:])
            gamma[i][k]=LP
    max_gamma = np.max(gamma)
    for i in range(m):
        gamma[i] = np.exp(gamma[i] - np.max(gamma[i]))
        gamma[i] /= np.sum(gamma[i])
    gamma[gamma<1e-300]=1e-300
    total_responsibilities_k=np.sum(gamma,axis=1)
    for i in range(m):
                gamma[i,:]/=total_responsibilities_k[i]
    P=np.sum(np.log(total_responsibilities_k))+m*max_gamma
    return gamma,P
def re_estimate_parameters(X,gm,gamma):
    m,n=X.shape
    row_sum_gamma=np.sum(gamma,axis=0)
    alpha=row_sum_gamma/np.sum(row_sum_gamma)
    for k in range(MIXTURE_COUNT):
        mean=np.zeros(n)
        for i in range(m):
            mean+=gamma[i,k]*X[i,:]
        mean/=row_sum_gamma[k]
        cov=np.zeros((n,n))
        for i in range(m):
            X1=np.reshape(X[i,:]-mean,(n,1))
            cov+=gamma[i,k]*(X1@X1.T)
        cov/=row_sum_gamma[k]
        gm[k]=multivariate_gaussian(mean,cov)
    return gm,alpha
def train_guassian_mixture(X):
    alpha=np.random.uniform(0,1,MIXTURE_COUNT)
    alpha=alpha/np.sum(alpha)
    mean=np.mean(X,0)
    cov=np.cov(X.T)
    gm = np.array([multivariate_gaussian(mean, cov) for _ in range(MIXTURE_COUNT)])
    prev_P=-np.inf
    for i in range(7):
        gamma,P = expectation(X,gm,alpha)
        if P - prev_P < -1e4:
            break
        gm,alpha = re_estimate_parameters(X,gm,gamma)
        prev_P=P
        print(f"i: {i}")
    return gm,alpha
def train_guassian_mixtures_model(X_train,y_train):
    class_labels=np.unique(y_train)
    n=class_labels.size
    alphas=[]
    gms=[]
    for i in range(10):
        gm,alpha=train_guassian_mixture(X_train[np.where(y_train==f'{i}')])
        print(f"{i} class done")
        alphas.append(alpha)
        gms.append(gm)
    return class_labels,gms,alphas
def test_guassian_mixture_model(X_test, class_labels, gms, alphas):
    m = X_test.shape[0]
    num_classes = len(class_labels)
    log_likelihoods = np.zeros((m, num_classes))
    print("\n--- Testing Model ---")
    for i in range(m):
        if i > 0 and i % 2000 == 0:
            print(f"Testing sample {i}/{m}")
        x_i = X_test[i, :]
        for c_idx, label in enumerate(class_labels):
            gm_c = gms[c_idx]
            alpha_c = alphas[c_idx]
            if gm_c is None:
                log_likelihoods[i, c_idx] = -np.inf
                continue
            log_probs_per_component = np.zeros(MIXTURE_COUNT)
            for k in range(MIXTURE_COUNT):
                log_probs_per_component[k] = np.log(alpha_c[k]) + gm_c[k].logpdf(x_i)
            max_log_prob = np.max(log_probs_per_component)
            sum_exp = np.sum(np.exp(log_probs_per_component - max_log_prob))
            log_sum = np.log(sum_exp)
            log_likelihoods[i, c_idx] = max_log_prob + log_sum
    predicted_indices = np.argmax(log_likelihoods, axis=1)
    y_pred = class_labels[predicted_indices]
    return y_pred,log_likelihoods

def roc_curve(log_likelihoods, y, class_labels):
    m = log_likelihoods.shape[0]
    plt.figure(figsize=(15, 10))

    for i, label in enumerate(class_labels):
        y_binary = (y == label).astype(int)
        scores = log_likelihoods[:, i]
        sorted_indices = np.argsort(scores)[::-1] 
        y_sorted = y_binary[sorted_indices]
        P = np.sum(y_sorted == 1)
        N = np.sum(y_sorted == 0)
        tp = np.cumsum(y_sorted == 1)
        fp = np.cumsum(y_sorted == 0)
        pd = np.concatenate([[0], tp / P]) if P > 0 else np.zeros(len(tp) + 1)
        pfa = np.concatenate([[0], fp / N]) if N > 0 else np.zeros(len(fp) + 1)

        plt.subplot(2, 5, i + 1)
        plt.plot(pfa, pd, label=f"Class {label}")
        plt.plot([0, 1], [0, 1], 'k--', lw=0.8)
        plt.title(f"Class {label}")
        plt.xlabel("Pfa")
        plt.ylabel("Pd")
        plt.grid(True)
    
    plt.suptitle("ROC Curves", fontsize=16)
    plt.tight_layout(rect=[0, 0, 1, 0.95])
    plt.show()


mnist = fetch_openml('mnist_784',version=1,as_frame=False)
(X,y) = (mnist.data,mnist.target)
X_train,X_test,y_train,y_test=train_test_split(X,y,test_size=0.3,random_state=10)
pca = PCA(n_components=CLASS_SIZE)
pca.fit(X_train)
X_train = pca.transform(X_train)
X_test = pca.transform(X_test)
MODEL_PATH = "gmm_model.pkl"

if os.path.exists(MODEL_PATH):
    print("Loading saved GMM model...")
    with open(MODEL_PATH, "rb") as f:
        class_labels, gms, alphas = pickle.load(f)
else:
    print("Training GMM model...")
    class_labels, gms, alphas = train_guassian_mixtures_model(X_train, y_train)
    with open(MODEL_PATH, "wb") as f:
        pickle.dump((class_labels, gms, alphas), f)
    print("Model saved successfully.")
y_pred,log_likelihoods = test_guassian_mixture_model(X_test, class_labels, gms, alphas)
roc_curve(log_likelihoods,y_test,class_labels)
accuracy = np.mean(y_pred == y_test)
error_rate = 1.0 - accuracy

print("\n--- Results ---")
print(f"Predicted labels (first 20): {y_pred[:20]}")
print(f"Actual labels (first 20):   {y_test[:20]}")
print(f"\nAccuracy: {accuracy * 100:.2f}%")
print(f"Error Rate: {error_rate * 100:.2f}%")