title: How can implement EM-GMM in python?
body: <p>I have implemented <a href="https://en.wikipedia.org/wiki/Expectation%E2%80%93maximization_algorithm" rel="nofollow noreferrer">EM algorithm</a> for <a href="https://scikit-learn.org/stable/modules/mixture.html" rel="nofollow noreferrer">GMM</a> using this post <a href="https://towardsdatascience.com/how-to-code-gaussian-mixture-models-from-scratch-in-python-9e7975df5252" rel="nofollow noreferrer">GMMs and Maximum Likelihood Optimization Using NumPy</a> unsuccessfully as follows:</p>
  <pre><code>import numpy as np
  
  def PDF(data, means, variances):
      return 1/(np.sqrt(2 * np.pi * variances) + eps) * np.exp(-1/2 * (np.square(data - means) / (variances + eps)))
  
  def EM_GMM(data, k, iterations):
      weights = np.ones((k, 1)) / k # shape=(k, 1)
      means = np.random.choice(data, k)[:, np.newaxis] # shape=(k, 1)
      variances = np.random.random_sample(size=k)[:, np.newaxis] # shape=(k, 1)
  
      data = np.repeat(data[np.newaxis, :], k, 0) # shape=(k, n)
  
      for step in range(iterations):
          # Expectation step
          likelihood = PDF(data, means, np.sqrt(variances)) # shape=(k, n)
  
          # Maximization step
          b = likelihood * weights # shape=(k, n)
          b /= np.sum(b, axis=1)[:, np.newaxis] + eps
  
          # updage means, variances, and weights
          means = np.sum(b * data, axis=1)[:, np.newaxis] / (np.sum(b, axis=1)[:, np.newaxis] + eps)
          variances = np.sum(b * np.square(data - means), axis=1)[:, np.newaxis] / (np.sum(b, axis=1)[:, np.newaxis] + eps)
          weights = np.mean(b, axis=1)[:, np.newaxis]
          
      return means, variances
  </code></pre>
  <p>when I run the algorithm on a 1-D time-series dataset, for k equal to 3, it returns an output like the following:</p>
  <pre><code>array([[0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
      3.05053810e-003, 2.36989898e-025, 2.36989898e-025,
      1.32797395e-136, 6.91134950e-031, 5.47347807e-001,
      1.44637007e+000, 1.44637007e+000, 1.44637007e+000,
      1.44637007e+000, 1.44637007e+000, 1.44637007e+000,
      1.44637007e+000, 1.44637007e+000, 1.44637007e+000,
      1.44637007e+000, 1.44637007e+000, 1.44637007e+000,
      1.44637007e+000, 2.25849208e-064, 0.00000000e+000,
      1.61228562e-303, 0.00000000e+000, 0.00000000e+000,
      0.00000000e+000, 0.00000000e+000, 3.94387272e-242,
      1.13078186e+000, 2.53108878e-001, 5.33548114e-001,
      9.14920432e-001, 2.07015697e-013, 4.45250680e-038,
      1.43000602e+000, 1.28781615e+000, 1.44821615e+000,
      1.18186109e+000, 3.21610659e-002, 3.21610659e-002,
      3.21610659e-002, 3.21610659e-002, 3.21610659e-002,
      2.47382844e-039, 0.00000000e+000, 2.09150855e-200,
      0.00000000e+000, 0.00000000e+000],
     [5.93203066e-002, 1.01647068e+000, 5.99299162e-001,
      0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
      0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
      0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
      0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
      0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
      0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
      0.00000000e+000, 0.00000000e+000, 2.14690238e-010,
      2.49337135e-191, 5.10499986e-001, 9.32658804e-001,
      1.21148135e+000, 1.13315278e+000, 2.50324069e-237,
      0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
      0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
      0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
      0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
      0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
      0.00000000e+000, 1.73966953e-125, 2.53559290e-275,
      1.42960975e-065, 7.57552338e-001],
     [0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
      3.05053810e-003, 2.36989898e-025, 2.36989898e-025,
      1.32797395e-136, 6.91134950e-031, 5.47347807e-001,
      1.44637007e+000, 1.44637007e+000, 1.44637007e+000,
      1.44637007e+000, 1.44637007e+000, 1.44637007e+000,
      1.44637007e+000, 1.44637007e+000, 1.44637007e+000,
      1.44637007e+000, 1.44637007e+000, 1.44637007e+000,
      1.44637007e+000, 2.25849208e-064, 0.00000000e+000,
      1.61228562e-303, 0.00000000e+000, 0.00000000e+000,
      0.00000000e+000, 0.00000000e+000, 3.94387272e-242,
      1.13078186e+000, 2.53108878e-001, 5.33548114e-001,
      9.14920432e-001, 2.07015697e-013, 4.45250680e-038,
      1.43000602e+000, 1.28781615e+000, 1.44821615e+000,
      1.18186109e+000, 3.21610659e-002, 3.21610659e-002,
      3.21610659e-002, 3.21610659e-002, 3.21610659e-002,
      2.47382844e-039, 0.00000000e+000, 2.09150855e-200,
      0.00000000e+000, 0.00000000e+000]])
  </code></pre>
  <p>which I believe is working wrong since the outputs are two vectors which one of them represents <code>means</code> values and the other one represents <code>variances</code> values. The vague point which made me doubtful about implementation is it returns back <code>0.00000000e+000</code> for most of the outputs as it can be seen and it doesn't need really to visualize these outputs. BTW the input data are <strong>time-series</strong> data. I have checked everything and traced multiple times but no bug shows up.</p>
  <p>Here are my input data:</p>
  <pre><code>[25.31      , 24.31      , 24.12      , 43.46      , 41.48666667,
     41.48666667, 37.54      , 41.175     , 44.81      , 44.44571429,
     44.44571429, 44.44571429, 44.44571429, 44.44571429, 44.44571429,
     44.44571429, 44.44571429, 44.44571429, 44.44571429, 44.44571429,
     44.44571429, 44.44571429, 39.71      , 26.69      , 34.15      ,
     24.94      , 24.75      , 24.56      , 24.38      , 35.25      ,
     44.62      , 44.94      , 44.815     , 44.69      , 42.31      ,
     40.81      , 44.38      , 44.56      , 44.44      , 44.25      ,
     43.66666667, 43.66666667, 43.66666667, 43.66666667, 43.66666667,
     40.75      , 32.31      , 36.08      , 30.135     , 24.19      ]
  </code></pre>
  <p>I was wondering if there is an elegant way to implement it via <code>numpy</code> or <code>SciKit-learn</code>. Any helps will be appreciated.</p>
  <p><strong>Update</strong>
  Following is current output and expected output:
  <img src="https://i.imgur.com/4Otlrsx.png" alt="img" /></p>