title: Expectation Maximization Algorithm (EM) for Gaussian Mixture Models (GMMs)
body: <p>I'm trying to apply the Expectation Maximization Algorithm (EM) to a Gaussian Mixture Model (GMM) using Python and NumPy. The PDF document I am basing my implementation on can be found <a href="https://www.ics.uci.edu/%7Esmyth/courses/cs274/notes/EMnotes.pdf" rel="nofollow noreferrer">here</a>.
  Below are the equations:</p>
  <p><a href="https://i.stack.imgur.com/YKGPM.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/YKGPM.png" alt="enter image description here" /></a></p>
  <p>When applying the algorithm I get the mean of the first and second cluster equal to:</p>
  <pre><code>array([[2.50832195],
         [2.51546208]])
  </code></pre>
  <p>When the actual vector means for the first and second cluster are, respectively:</p>
  <pre><code>array([[0],
         [0]])
  </code></pre>
  <p>and:</p>
  <pre><code>array([[5],
         [5]])
  </code></pre>
  <p>The same thing happens when getting the values of the covariance matrices I get:</p>
  <pre><code>array([[7.05168736, 6.17098629],
         [6.17098629, 7.23009494]])
  </code></pre>
  <p>When it should be:</p>
  <pre><code>array([[1, 0],
         [0, 1]])
  </code></pre>
  <p>for both clusters.
  Here is the code:</p>
  <pre><code>np.random.seed(1)
  
  # first cluster
  X_11 = np.random.normal(0, 1, 1000)
  X_21 = np.random.normal(0, 1, 1000)
  
  # second cluster
  X_12 = np.random.normal(5, 1, 1000)
  X_22 = np.random.normal(5, 1, 1000)
  
  X_1 = np.concatenate((X_11,X_12), axis=None)
  X_2 = np.concatenate((X_21,X_22), axis=None)
  
  # data matrix of k x n dimensions (2 x 2000 dimensions)
  X = np.concatenate((np.array([X_1]),np.array([X_2])), axis=0)
  
  # multivariate normal distribution function gives n x 1 vector (2000 x 1 vector)
  def normal_distribution(x, mu, sigma):
    mvnd = []
    for i in range(np.shape(x)[1]):
      gd = (2*np.pi)**(-2/2) * np.linalg.det(sigma)**(-1/2) * np.exp((-1/2) * np.dot(np.dot((x[:,i:i+1]-mu).T, np.linalg.inv(sigma)), (x[:,i:i+1]-mu)))
      mvnd.append(gd)
    return np.reshape(np.array(mvnd), (np.shape(x)[1], 1))
  
  # Initialized parameters
  sigma_1 = np.array([[10, 0],
                      [0, 10]])
  sigma_2 = np.array([[10, 0],
                      [0, 10]])
  mu_1 = np.array([[10], 
                   [10]])
  mu_2 = np.array([[10], 
                   [10]])
  pi_1 = 0.5
  pi_2 = 0.5
  
  Sigma_1 = np.empty([2000, 2, 2])
  Sigma_2 = np.empty([2000, 2, 2])
  
  for i in range(10):
    # E-step:
    w_i1 = (pi_1*normal_distribution(X, mu_1, sigma_1))/(pi_1*normal_distribution(X, mu_1, sigma_1) + pi_2*normal_distribution(X, mu_2, sigma_2))
    w_i2 = (pi_2*normal_distribution(X, mu_2, sigma_2))/(pi_1*normal_distribution(X, mu_1, sigma_1) + pi_2*normal_distribution(X, mu_2, sigma_2))
    # M-step:
    pi_1 = np.sum(w_i1)/2000
    pi_2 = np.sum(w_i2)/2000
    mu_1 = np.array([(1/(np.sum(w_i1)))*np.sum(w_i1.T*X, axis=1)]).T
    mu_2 = np.array([(1/(np.sum(w_i2)))*np.sum(w_i2.T*X, axis=1)]).T
    for i in range(2000):
      Sigma_1[i:i+1, :, :] = w_i1[i:i+1,:]*np.dot((X[:,i:i+1]-mu_1), (X[:,i:i+1]-mu_1).T)
      Sigma_2[i:i+1, :, :] = w_i2[i:i+1,:]*np.dot((X[:,i:i+1]-mu_2), (X[:,i:i+1]-mu_2).T)
      sigma_1 = (1/(np.sum(w_i1)))*np.sum(Sigma_1, axis=0)
      sigma_2 = (1/(np.sum(w_i2)))*np.sum(Sigma_2, axis=0)
  </code></pre>
  <p>Would really appreciate if someone could point out the mistake in my code or in my misunderstanding of the algorithm..</p>