#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Fri Dec 30 12:52:50 2022 @author: ojacques """ # For plotting import matplotlib.pyplot as plt import seaborn as sns sns.set_style("white") #%matplotlib inline #for matrix math import numpy as np #for normalization + probability density function computation from scipy import stats #for data preprocessing import pandas as pd from math import sqrt, log, exp, pi from random import uniform print("import done") #============================================================================= random_seed=36788765 np.random.seed(random_seed) Mean1 = 2.0 # Input parameter, mean of first normal probability distribution Standard_dev1 = 4.0 #@param {type:"number"} Mean2 = 9.0 # Input parameter, mean of second normal probability distribution Standard_dev2 = 2.0 #@param {type:"number"} # generate data y1 = np.random.normal(Mean1, Standard_dev1, 1000) y2 = np.random.normal(Mean2, Standard_dev2, 500) data=np.append(y1,y2) # For data visiualisation calculate left and right of the graph Min_graph = min(data) Max_graph = max(data) x = np.linspace(Min_graph, Max_graph, 2000) # to plot the data print(f'Input Gaussian 1: μ = {Mean1:.2}, σ = {Standard_dev1:.2}') print(f'Input Gaussian 1: μ = {Mean2:.2}, σ = {Standard_dev2:.2}') sns.distplot(data, bins=20, kde=False); #============================================================================== class Gaussian: "Model univariate Gaussian" def __init__(self, mu, sigma): #mean and standard deviation self.mu = mu self.sigma = sigma #probability density function def pdf(self, datum): "Probability of a data point given the current parameters" u = (datum - self.mu) / abs(self.sigma) y = (1 / (sqrt(2 * pi) * abs(self.sigma))) * exp(-u * u / 2) return y def __repr__(self): return f'Gaussian({self.mu:4.6}, {self.sigma:4.6}' print("done") #============================================================================== #gaussian of best fit best_single = Gaussian(np.mean(data), np.std(data)) print('Best single Gaussian: μ = {:.2}, σ = {:.2}'.format(best_single.mu, best_single.sigma)) #fit a single gaussian curve to the data g_single = stats.norm(best_single.mu, best_single.sigma).pdf(x) sns.distplot(data, bins=20, kde=False, norm_hist=True); plt.plot(x, g_single, label='single gaussian'); plt.legend(); (rigt) # compute new means self.one.mu = sum(w * d for (w, d) in zip(left, data)) / one_den self.two.mu = sum(w * d for (w, d) in zip(rigt, data)) / two_den # compute new sigmas self.one.sigma = sqrt(sum(w * ((d - self.one.mu) ** 2) for (w, d) in zip(left, data)) / one_den) self.two.sigma = sqrt(sum(w * ((d - self.two.mu) ** 2) for (w, d) in zip(rigt, data)) / two_den) # compute new mix self.mix = one_den / len(data) def iterate(self, N=1, verbose=False): "Perform N iterations, then compute log-likelihood" for i in range(1, N+1): self.Mstep(self.Estep()) #The heart of the algorith, perform E-stepand next M-step if verbose: print(f'{i:2} {self}') self.Estep() # to freshen up self.loglike def pdf(self, x): return (self.mix)*self.one.pdf(x) + (1-self.mix)*self.two.pdf(x) def __repr__(self): return f'GaussianMixture({self.one}, {self.two}, mix={self.mix:.03}' def __str__(self): return f'Mixture: {self.one}, {self.two}, mix={self.mix:.03}' print("done") #============================================================================== class GaussianMixture_self: "Model mixture of two univariate Gaussians and their EM estimation" def __init__(self, data, mu_min=min(data), mu_max=max(data), sigma_min=1, sigma_max=1, mix=.5): self.data = data #todo the Algorithm would be numerical enhanced by normalizing the data first, next do all the EM steps and do the de-normalising at the end #init with multiple gaussians self.one = Gaussian(uniform(mu_min, mu_max), uniform(sigma_min, sigma_max)) self.two = Gaussian(uniform(mu_min, mu_max), uniform(sigma_min, sigma_max)) #as well as how much to mix them self.mix = mix def Estep(self): "Perform an E(stimation)-step, assign each point to gaussian 1 or 2 with a percentage" # compute weights self.loglike = 0. # = log(p = 1) for datum in self.data: # unnormalized weights wp1 = self.one.pdf(datum) * self.mix wp2 = self.two.pdf(datum) * (1. - self.mix) # compute denominator den = wp1 + wp2 # normalize wp1 /= den wp2 /= den # wp1+wp2= 1, it either belongs to gaussian 1 or gaussion 2 # add into loglike self.loglike += log(den) #freshening up self.loglike in the process # yield weight tuple yield (wp1, wp2) def Mstep(self, weights): "Perform an M(aximization)-step" # compute denominators (left, rigt) = zip(*weights) one_den = sum(left) two_den = sum(rigt) # compute new means self.one.mu = sum(w * d for (w, d) in zip(left, data)) / one_den self.two.mu = sum(w * d for (w, d) in zip(rigt, data)) / two_den # compute new sigmas self.one.sigma = sqrt(sum(w * ((d - self.one.mu) ** 2) for (w, d) in zip(left, data)) / one_den) self.two.sigma = sqrt(sum(w * ((d - self.two.mu) ** 2) for (w, d) in zip(rigt, data)) / two_den) # compute new mix self.mix = one_den / len(data) def iterate(self, N=1, verbose=False): "Perform N iterations, then compute log-likelihood" for i in range(1, N+1): self.Mstep(self.Estep()) #The heart of the algorith, perform E-stepand next M-step if verbose: print(f'{i:2} {self}') self.Estep() # to freshen up self.loglike def pdf(self, x): return (self.mix)*self.one.pdf(x) + (1-self.mix)*self.two.pdf(x) def __repr__(self): return f'GaussianMixture({self.one}, {self.two}, mix={self.mix:.03}' def __str__(self): return f'Mixture: {self.one}, {self.two}, mix={self.mix:.03}' print("done") #============================================================================== # See the algorithm in action ============================================================================== n_iterations = 20 best_mix = None best_loglike = float('-inf') mix = GaussianMixture_self(data) for _ in range(n_iterations): try: #train! mix.iterate(verbose=True) if mix.loglike > best_loglike: best_loglike = mix.loglike best_mix = mix except (ZeroDivisionError, ValueError, RuntimeWarning): # Catch division errors from bad starts, and just throw them out... print("one less") pass