import numpy as np
import scipy.stats as sps
import matplotlib.pyplot as plt
from scipy.misc import logsumexp
#%matplotlib qt5
%matplotlib inline
plt.rcParams.update({'figure.figsize': (10.0, 6.0), 'font.size': 18})

Gaussian Mixture Model (GMM)¶

probability density finction $p (x | {μ_{c}}, {σ_{c}^{2}}, {π_{c}}) = \sum_{c = 1}^{C} N (x | μ_{c}, σ_{c}^{2}) π_{c}$ $\large \DeclareMathOperator{\Norm}{\mathcal{N}} \DeclareMathOperator{\Gam}{Gam} \DeclareMathOperator{\e}{exp} p(x|\{\mu_c\}, \{\sigma_c^2\}, \{\pi_c\}) = \sum_{c=1}^{C} \Norm(x|\mu_c, \sigma^2_c) \pi_c$

where

$\{\mu_c\}$ is the set of $C$ means
$\{\sigma_c^2\}$ is the set of $C$ variances
$\{\pi_c\}$ is set of $C$ weights such that $\sum_{c=1}^C \pi_c = 1$

and single variate Gasussian distribution

p (x ∣ μ, σ^{2}) = N (x ∣ μ, σ^{2}) = \frac{1}{\sqrt{2 π σ^{2}}} \exp {\frac{- (x - μ)^{2}}{2 σ^{2}}}

$\large p(x \mid \mu, \sigma^2) = \Norm(x \mid \mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \e\left\{ \frac{-(x - \mu)^2}{2\sigma^2} \right\}$

# Plot GMM pdf together with the individual components pdfs; return the GMM pdf line
def plot_GMM(t, mus, sigmas, pis):
  p_xz = sps.norm.pdf(t[:,np.newaxis], mus, sigmas) * pis # all GMM components are evaluated at once
  px = np.sum(p_xz, axis=1)
  plt.plot(t, p_xz, ':')
  plt.plot(t, px, 'k')
  return px


#Handcraft some GMM parameter 
mus = [-4.0, 0.0, 4.0, 5]
sigmas = [1.0, 1.4, 1.2, 1]
pis = [0.1, 0.4, 0.2, 0.3]

t = np.linspace(-10,10,1000)
true_GMM_pdf = plot_GMM(t, mus, sigmas, pis)

# Generate N datapoints from this GMM
N = 100
Nc = sps.multinomial.rvs(N, pis) # Draw observation counts for each component from multinomial distribution
x = sps.norm.rvs(np.repeat(mus, Nc), np.repeat(sigmas, Nc))
np.random.shuffle(x)
plt.plot(x, np.zeros_like(x), '+k');

GMM - EM algorithm¶

E-step

γ_{n c} = P (z_{n} = c | x_{n}, η^{o l d}) = \frac{p (x_{n} | z_{n} = c, η^{o l d}) P (z_{n} = c | η^{o l d})}{p (x_{n} | η^{o l d})} = \frac{N (x_{n} | μ_{c}^{o l d} {σ_{c}^{2}}^{o l d}) π_{c}^{o l d}}{\sum_{k} N (x_{n} | μ_{k}^{o l d}, {σ_{k}^{2}}^{o l d}) π_{k}^{o l d}}

$\large \DeclareMathOperator{\Norm}{\mathcal{N}} \DeclareMathOperator{\eeta}{\boldsymbol{\eta}} \gamma_{nc}=P(z_n=c|x_n,\eeta^{old})=\frac{p(x_n|z_n=c,\eeta^{old})P(z_n=c|\eeta^{old})}{p(x_n|\eeta^{old})} =\frac{\mathcal{N}(x_n|\mu_c^{old}{\sigma_c^2}^{old})\pi_c^{old}}{\sum_{k}\mathcal{N}(x_n|\mu_{k}^{old},{\sigma_{k}^2}^{old})\pi_{k}^{old}}$

M-step

\begin{aligned} μ_{c}^{n e w} & = \frac{\sum_{n} γ_{n c} x_{n}}{\sum_{n} γ_{n c}} \\ {σ_{c}^{2}}^{n e w} & = \frac{\sum_{n} γ_{n c} (x_{n} - μ_{c}^{n e w})^{2}}{\sum_{n} γ_{n c}} \\ π_{c} & = \frac{\sum_{n} γ_{n c}}{N} \end{aligned}

$\large \begin{align} \mu_c^{new} & =\frac{\sum_n \gamma_{nc}x_n}{\sum_n \gamma_{nc}}\\ {\sigma_c^2}^{new} & =\frac{\sum_n \gamma_{nc}(x_n-\mu_c^{new})^2}{\sum_n \gamma_{nc}}\\ \pi_c & =\frac{\sum_n \gamma_{nc}}{N} \end{align}$

#Choose some initial parameters
C = 3        # number of GMM components 
mus = x[:C]  # we choose few first observations as the initial means
sigmas = np.repeat(np.std(x), C) # sigma for all components is set to std of the the training data
pis = np.ones(C)/C

plt.clf()
plt.plot(t, true_GMM_pdf, 'gray')
plot_GMM(t, mus, sigmas, pis);

for _ in range(50):
  #E-step
  log_p_xz = sps.norm.logpdf(x[:,np.newaxis], mus, sigmas) + np.log(pis)
  log_p_x = logsumexp(log_p_xz, axis=1, keepdims=True)
  print "Training data log likelihood:", log_p_x.sum()

  gammas = np.exp(log_p_xz - log_p_x)
  #M-step
  Nc = gammas.sum(axis=0)
  mus =  x.dot(gammas) / Nc
  sigmas =  np.sqrt((x**2).dot(gammas) / Nc - mus**2) # we use std, not variance!
  pis = Nc / Nc.sum()
    
plot_GMM(t, mus, sigmas, pis)

plt.clf()
plt.plot(t, true_GMM_pdf, 'gray')
plot_GMM(t, mus, sigmas, pis);
plt.plot(x, np.zeros_like(x), '+k');

Training data log likelihood: -275.687446231
Training data log likelihood: -257.756537213
Training data log likelihood: -256.73326682
Training data log likelihood: -255.418278557
Training data log likelihood: -253.868553645
Training data log likelihood: -252.588253235
Training data log likelihood: -251.9249919
Training data log likelihood: -251.682616958
Training data log likelihood: -251.604724967
Training data log likelihood: -251.57672686
Training data log likelihood: -251.562926418
Training data log likelihood: -251.553611058
Training data log likelihood: -251.546153868
Training data log likelihood: -251.539755108
Training data log likelihood: -251.534101033
Training data log likelihood: -251.529025116
Training data log likelihood: -251.52441792
Training data log likelihood: -251.520199681
Training data log likelihood: -251.516309709
Training data log likelihood: -251.51270094
Training data log likelihood: -251.509336448
Training data log likelihood: -251.506186972
Training data log likelihood: -251.503229061
Training data log likelihood: -251.500443693
Training data log likelihood: -251.497815217
Training data log likelihood: -251.495330568
Training data log likelihood: -251.492978667
Training data log likelihood: -251.490749977
Training data log likelihood: -251.488636164
Training data log likelihood: -251.486629848
Training data log likelihood: -251.484724409
Training data log likelihood: -251.482913848
Training data log likelihood: -251.48119267
Training data log likelihood: -251.479555809
Training data log likelihood: -251.47799856
Training data log likelihood: -251.476516528
Training data log likelihood: -251.475105591
Training data log likelihood: -251.473761871
Training data log likelihood: -251.472481707
Training data log likelihood: -251.471261635
Training data log likelihood: -251.470098376
Training data log likelihood: -251.468988817
Training data log likelihood: -251.467930006
Training data log likelihood: -251.466919136
Training data log likelihood: -251.465953544
Training data log likelihood: -251.465030697
Training data log likelihood: -251.464148189
Training data log likelihood: -251.463303737
Training data log likelihood: -251.462495171
Training data log likelihood: -251.461720434