# Expectation Maximization Algorithm

In this article, exctation maximization(EM) algorithm will be introduced. EM algorithm is introduced as early as 1950 by Ceppellini etal. in the context of gene frequency estimation, the expectation maximization algorithm was nanalyzed more generally by Hartley and by Baum et al. In this article, I will begin my story with an example of flipping coins and maximum likelihood estimator(MLE). This example is from the paper [what is the expectation maximization algorithm?](https://www.nature.com/articles/nbt1406)

## 1.1 Example of flipping coins

Suppose we have two coins A and B with unknown biases, $\theta_{A}$ and $\theta_{B}$ respectively. We do the experiemt that we fip the coins repeatedly. During the experiment, we keep track the record of the $x={x_{1},x_{2},...,x_{n}}$ and $z={z_{1},z_{2},...,z_{n}}$, where $x$ indicates the number of heads we observe when we fip the coin and $z$ imply the identity of coin we use to fip. Given $x$ and $z$, we want to estimated biases $\theta$ of these two coins. Obviously, we can use MLE to estimate the biase $\theta$ based on the bays theorem:

$$\begin{equation*}
p(\theta \vert x, z)=\frac{p(\theta)p(x,z\vert \theta)}{p(x,z)}
\end{equation*}$$

where we assume the normalization term and prior term are constant, thus to maximize the posterior is equivalent to maximize the likelihood $p(x,z\vert \theta)$. That's a **commom parameter estimation problem.** This kind of estimated problem is called parameter estimated with complete data.

Now, we consider a more challenging variant problem that we still we want to estimate baises of two coins A and B, but in this time, we only has the record $x$ and do not know about the $z$, where $z$ is called hidden parameter or latent factor. Parameter estimated problem with complete data is converted into parameter estimated problem with imcomplete data, so the MLE is not effective for the problem of imcomplete data. This is a big limitation of MLE and thus we need another method to sovle for imcomplete data problem.

However, if we have some way to know the value of $z$, we can use MLE to inference $\theta$ which is the same as parameter estimation with complete data. Thus, we can reduce parameter estimation for this problem with incomplete data to maximum likelihood estimation with complete data.

One iterative scheme for obtaining completions caould works as follows: starting from some initial parameters, $\widetilde{\theta}^{(t)}=(\widetilde{\theta}_{A}^{(t)},\widetilde{\theta}_{B}^{(t)})$, determine for each of the five sets whether coin $A$ or coin $B$ was more liekly to have generated the observed flips(using the current parameter estimates). Then, assume these completions (that is , guessed coin assignments) to be correct, and apply the regular maximum likelihood estimation procedure to get $\widetilde{\theta}^{(t+1)}$. Finally, repeat these two steps until convergence. 

The expectation maximization algorithm is a refinement on this basic idea that the expectation maximization algorithm computes probabilities for each possible completion of the missing data, using the current parameters $\widetilde{\theta}^{(t)}$. These probabilities are used to create a weighted training set consisting of all possible completions of the data. Then, a modified version of maximum likelihood estimation that deals with weighted training examples provides new parameter estimates, $\widetilde{\theta}^{(t+1)}$. 

The figure below show the procedure of maximum likelihood estimator and expectation maximization algorithm. The figure is from the paper:[ what is the expectation maximization algorithm?](https://www.nature.com/articles/nbt1406)



### Summary

**Intuitively, the expectation maximization algorithm alternates between the steps of guessing a probability distribution over completions of missing data given the current model( known as the E-step). and then re-estimating the model parameters using these completions(known as the M-step).**

## 1.2 Mathematical foundations(Convergence at intuitive level)

In the complete data case, the objective function $\log(P(x,z\vert \theta))$ has a single global optimum, which can often be found in closed form. In contrast, in the incomplete data case the function $\log(P(x\vert \theta))$ has multiple local maximum and no closed form solution.

To deal with this, the expectation maximization algorithm reduces the difficult task of optimization $\log(P(x\vert \theta))$ into a sequence of simpler optimization subproblems, whose objective functions have unique global maxima that can often be computed in closed form. These subproblems are chosen in a way that gurantees their corresponding solutions $\widetilde{\theta}^{(1)},\widetilde{\theta}^{(2)},...$ and will converge to a local optimmum of $\log(P(x\vert \theta))$.

More specifically, the expectation maximization algorithm alternates between two phases. During the E-step, expectation maximization chooses a function $g_{t}$, that lower bounds $\log(P(x\vert \theta))$ everywhere, and for which $g_{t}(\widetilde{\theta}^{(t)}=\log(P(x\vert \widetilde{\theta}^{(t)}))$. During the M-step, the expectation maximization algorithm moves to a new parameter set $\widetilde{\theta}^{(t+1)}$ that maximizes $g_{t}$. As the value of the lower-bound $g_{t}$ matches the objective function at $\widetilde{\theta}^{(t)}$, it follows that $\log(P(x\vert \widetilde{\theta}^{(t)}))=g_{t}(\widetilde{\theta}^{(t)})\leq g_{t}(\widetilde{\theta}^{(t+1)})=\log(P(x\vert \widetilde{\theta}^{(t+1)}))$. So the objective function monotonically increases during each iteration of expectation maximization.

## 1.3 Jenson's inequality

Jenson's inequality states that for convex function $f$, if $\lambda \in [0,1]$, then

$$
\begin{equation*}
\lambda f(x_{1})+(1-\lambda)f(x_{2}) \geq f(\lambda x_{1}+(1-\lambda)x_{2})
\end{equation*}
$$

Gerenrally, the formula above can be rewritten as:

$$
\begin{equation*}
\mathbb{E}[f(x)]\geq f(\mathbb{E}[x])
\end{equation*}
$$

The physical meaning is shown as below:



### EM Algorithm and Jensen's inequality

We re-write the complete log likelihood function by multiplication it by $\frac{q(z)}{q(z)}$, where $q(z)$ represents an arbitrary distribution for the random variable $Z$.

$$
\begin{align*}
\log p(x\vert \theta) &= \log \sum_{Z}p(z\vert \theta)p(x\vert z,\theta)\frac{q(z)}{q(z)} \\
&= \log \mathbb{E}_{q}\bigg[\frac{p(z\vert \theta)p(x\vert z, \theta)}{q(z)}\bigg] \\
&\geq \mathbb{E}_{q}\bigg[\log \frac{p(z\vert \theta)p(x\vert z, \theta)}{q(z)}\bigg] \\ 
\end{align*}
$$

**Note that this function is concave function, so the symbol of inequality changes the direction**, 

let' s $\mathcal{L}(\theta \vert q)= \mathbb{E}_{q}\bigg[\log \frac{p(z\vert \theta)p(x\vert z, \theta)}{q(z)}\bigg]$, then

$$
\begin{equation}
\mathcal{L}(\theta \vert q) = \mathbb{E}_{q}[\log p(z\vert \theta)]+\mathbb{E}_{q}[\log p(x\vert z,\theta)] -\mathbb{E}_{q}[\log q(z)] 
\end{equation}
$$

This is the objective function of EM algorithm. Note that the third term in the equation is known as entropy.

## 1.4 The EM algorithm

The EM algorithm proceeds by coordinate ascent. At each iteration $t$, we have update two parameters: $q^{(t)}$ and $\theta^{(t)}$.

At the E-step: we update the posterior value $q$ of the random variable given the observations while holding $\theta^{(t)}$ fixed.

$$
\begin{align*}
q^{(t+1)} &= 
\arg\max_{q} \mathcal{L}(q,\theta^{(t)}) \\ 
&= p(z\vert x, \theta^{(t)})
\end{align*}
$$

At the M-step, we update the model parameters to maximize the expected complete log likelihood function.

$$
\begin{equation}
\theta^{(t+1)} = \arg\max_{\theta} \mathcal{L}(q^{(t+1)},\theta)
\end{equation}
$$

In detail, let's see how it works:

$$
\begin{align*}
\mathcal{L}(p(z\vert x, \theta),\theta) &= \sum_{Z}p(z\vert x, \theta) \log \frac{p(x,z\vert \theta)}{p(z\vert x, \theta)} \\
&= \sum_{Z}p(z\vert x, \theta)\log \frac{p(x,z\vert \theta)p(x\vert \theta)}{p(x,z\vert \theta)} \\
&= \sum_{Z}p(z\vert x, \theta)\log p(x\vert \theta) \\
&= \log p(x\vert \theta)\sum_{Z}p(x\vert x, \theta) \\
&= \log p(x\vert \theta)
\end{align*}
$$

Therefore, as we maximize the objective function, we are also maximizing the log likelihood function.



## There some codes about coin toss example(these examples are from [sta663](https://people.duke.edu/~ccc14/sta-663/EMAlgorithm.html))


```python
import os
import sys
import glob
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
%matplotlib inline
plt.style.use('ggplot')
np.random.seed(1234)

np.set_printoptions(formatter={'all':lambda x: '%.3f' % x})

from IPython.display import Image
from numpy.core.umath_tests import matrix_multiply as mm

from scipy.optimize import minimize
from scipy.stats import bernoulli, binom
```

## Solving for complete likelihood


```python
def neg_log_likelihood(theta, n, xs, zs):
    return -np.sum([binom(n, theta[z]).logpmf(x) for (x, z) in zip(xs, zs)])
```


```python
# each time we flip the coin repeatedly 10 times
num = 10

theta_A = 0.8 # the baise of coin A
theta_B = 0.4 # the baise of coin B

coin_A = bernoulli(theta_A)
coin_B = bernoulli(theta_B)

xs = np.array(list(map(sum, [coin_A.rvs(num), coin_B.rvs(num), coin_A.rvs(num), coin_A.rvs(num), coin_B.rvs(num)])))
zs = [0,1,0,0,1]
```

### Sovle by hand


```python
ml_A = np.sum(xs[[0,2,3]])/(3.0*num)
ml_B = np.sum(xs[[1,4]])/(2.0*num)

print(ml_A, ml_B)
```

    0.833333333333 0.45


### Solve by iterative algorithm


```python
bounds = [(0,1),(0,1)]
minimize(neg_log_likelihood, (0.5, 0.5), args=(num,xs,zs), method='tnc',
        bounds=bounds, options={'maxiter': 100})
```




         fun: 8.7695482351748506
         jac: array([-0.000, -0.000])
     message: 'Converged (|f_n-f_(n-1)| ~= 0)'
        nfev: 24
         nit: 7
      status: 1
     success: True
           x: array([0.833, 0.450])



## Solve for incomplete likelihood


```python
# we only know the data xs, and we need to guess the thetas firstly.
xs = np.array([(5,5), (9,1), (8,2), (4,6), (7,3)])
thetas = np.array([[0.6, 0.4], [0.5, 0.5]])

tol = 0.01 # stop criterion
num_maxiter = 100
ll_old = 0 #LogLikeliHood
for i in range(num_maxiter):
    ws_A = []
    ws_B = []
    
    vs_A = []
    vs_B = []
    
    ll_new=0
    for x in xs:
        # E-step: compute the probability distribution over all possible completions
        ll_A = np.sum([x*np.log(thetas[0])])
        ll_B = np.sum([x*np.log(thetas[1])])
        
        demo = np.exp(ll_A)+np.exp(ll_B)
        w_A = np.exp(ll_A)/demo
        w_B = np.exp(ll_B)/demo
        
        ws_A.append(w_A)
        ws_B.append(w_B)
        
        v_A = np.dot(w_A, x)
        v_B = np.dot(w_B, x)
        
        vs_A.append(v_A)
        vs_B.append(v_B)
        ll_new += w_A*ll_A +w_B*ll_B
    
    # M-step: update the parameter of theta based on the current model
    thetas[0] = np.sum(vs_A, 0)/np.sum(vs_A)
    thetas[1] = np.sum(vs_B, 0)/np.sum(vs_B)
    
    print('the iteration is %d'% (i+1))
    print('theta_A=%.2f, theta_B=%.2f, ll=%.2f' % (thetas[0,0], thetas[1,0], ll_new))
    
    if np.abs(ll_new-ll_old)<tol:
        break
    ll_old=ll_new
```

    the iteration is 1
    theta_A=0.71, theta_B=0.58, ll=-32.69
    the iteration is 2
    theta_A=0.75, theta_B=0.57, ll=-31.26
    the iteration is 3
    theta_A=0.77, theta_B=0.55, ll=-30.76
    the iteration is 4
    theta_A=0.78, theta_B=0.53, ll=-30.33
    the iteration is 5
    theta_A=0.79, theta_B=0.53, ll=-30.07
    the iteration is 6
    theta_A=0.79, theta_B=0.52, ll=-29.95
    the iteration is 7
    theta_A=0.80, theta_B=0.52, ll=-29.90
    the iteration is 8
    theta_A=0.80, theta_B=0.52, ll=-29.88
    the iteration is 9
    theta_A=0.80, theta_B=0.52, ll=-29.87



```python
# we only know the data xs, and we need to guess the thetas firstly.
xs = np.array([(5,5), (9,1), (8,2), (4,6), (7,3)])
thetas = np.array([[0.6, 0.4], [0.5, 0.5]])

tol = 0.01 # stop criterion
num_maxiter = 100
ll_old = -np.infty #LogLikeliHood
for i in range(num_maxiter):
    ws_A = []
    ws_B = []
    
    vs_A = []
    vs_B = []
    
    ll_new=0
    for x in xs:
        # E-step: compute the probability distribution over all possible completions
        ll_A = np.sum([x*np.log(thetas[0])])
        ll_B = np.sum([x*np.log(thetas[1])])
        
        demo = np.exp(ll_A)+np.exp(ll_B)
        w_A = np.exp(ll_A)/demo
        w_B = np.exp(ll_B)/demo
        
        ws_A.append(w_A)
        ws_B.append(w_B)
        
        v_A = np.dot(w_A, x)
        v_B = np.dot(w_B, x)
        
        vs_A.append(v_A)
        vs_B.append(v_B)
        ll_new += w_A*ll_A +w_B*ll_B
    
    # M-step: update the parameter of theta based on the current model
    thetas[0] = np.sum(vs_A, 0)/np.sum(vs_A)
    thetas[1] = np.sum(vs_B, 0)/np.sum(vs_B)
    
    print('the iteration is %d'% (i+1))
    print('theta_A=%.2f, theta_B=%.2f, ll=%.2f' % (thetas[0,0], thetas[1,0], ll_new))
    
    if np.abs(ll_new-ll_old)<tol:
        break
    ll_old=ll_new
```

    the iteration is 1
    theta_A=0.71, theta_B=0.58, ll=-32.69
    the iteration is 2
    theta_A=0.75, theta_B=0.57, ll=-31.26
    the iteration is 3
    theta_A=0.77, theta_B=0.55, ll=-30.76
    the iteration is 4
    theta_A=0.78, theta_B=0.53, ll=-30.33
    the iteration is 5
    theta_A=0.79, theta_B=0.53, ll=-30.07
    the iteration is 6
    theta_A=0.79, theta_B=0.52, ll=-29.95
    the iteration is 7
    theta_A=0.80, theta_B=0.52, ll=-29.90
    the iteration is 8
    theta_A=0.80, theta_B=0.52, ll=-29.88
    the iteration is 9
    theta_A=0.80, theta_B=0.52, ll=-29.87



```python
# we only know the data xs, and we need to guess the thetas firstly.
xs = np.array([(5,5), (9,1), (8,2), (4,6), (7,3)])
thetas = np.array([[0.6, 0.4], [0.5, 0.5]])

tol = 0.01 # stop criterion
num_maxiter = 100
ll_old = -np.infty #LogLikeliHood
for i in range(num_maxiter):
    # E-step: compute the probability distribution over all possible completions
    ll_A = np.sum(xs*np.log(thetas[0]), axis=1)
    ll_B = np.sum(xs*np.log(thetas[1]), axis=1)
        
    demo = np.exp(ll_A)+np.exp(ll_B)
    w_A = np.exp(ll_A)/demo
    w_B = np.exp(ll_B)/demo
        
    vs_A = w_A[:, None] * xs
    vs_B = w_B[:, None] * xs
    # M-step: update the parameter of theta based on the current model
    thetas[0] = np.sum(vs_A, 0)/np.sum(vs_A)
    thetas[1] = np.sum(vs_B, 0)/np.sum(vs_B)
    
    ll_new = w_A.dot(ll_A) + w_B.dot(ll_B)-w_A.dot(np.log(w_A))-w_B.dot(np.log(w_B))
    
    print('the iteration is %d'% (i+1))
    print('theta_A=%.3f, theta_B=%.3f, ll=%.3f' % (thetas[0,0], thetas[1,0], ll_new))
    
    if np.abs(ll_new-ll_old)<tol:
        break
    ll_old=ll_new
```

    the iteration is 1
    theta_A=0.713, theta_B=0.581, ll=-29.628
    the iteration is 2
    theta_A=0.745, theta_B=0.569, ll=-28.394
    the iteration is 3
    theta_A=0.768, theta_B=0.550, ll=-28.257
    the iteration is 4
    theta_A=0.783, theta_B=0.535, ll=-28.162
    the iteration is 5
    theta_A=0.791, theta_B=0.526, ll=-28.119
    the iteration is 6
    theta_A=0.795, theta_B=0.522, ll=-28.107
    the iteration is 7
    theta_A=0.796, theta_B=0.521, ll=-28.105



```python
def em(xs, thetas, maxiter=100, tol=1e-6):
    '''Expectation-maximization for coin sample problem'''
    
    ll_old = -np.infty
    for i in range(maxiter):
        ll = np.array([np.sum(xs * np.log(theta), axis=1) for theta in thetas])
        lik = np.exp(ll)
        ws = lik/lik.sum(0)
        vs = np.array([w[:,None] * xs for w in ws])
        thetas = np.array([v.sum(0)/v.sum() for v in vs])
        ll_new = np.sum([w*l for w, l in zip(ws, ll)])
        if np.abs(ll_new - ll_old) < tol:
            break
        ll_old = ll_new
    return i, thetas, ll_new
```


```python
# we only know the data xs, and we need to guess the thetas firstly.
xs = np.array([(5,5), (9,1), (8,2), (4,6), (7,3)])
thetas = np.array([[0.6, 0.4], [0.5, 0.5]])

i, thetas, ll = em(xs, thetas)
print(i)
for theta in thetas:
    print(theta)
print(ll)
```

    18
    [0.797 0.203]
    [0.520 0.480]
    -29.868676155



```python
np.random.seed(1234)

n = 100
p0 = 0.8
p1 = 0.35

xs = np.concatenate([np.random.binomial(n, p0, int(n/2)), np.random.binomial(n, p1, int(n/2))])
xs = np.column_stack([xs, n-xs])
np.random.shuffle(xs)
```


```python
results = [em(xs, np.random.random((2,2))) for i in range(10)]
i, thetas, ll = sorted(results, key=lambda x: x[-1])[-1]
print(i)
for theta in thetas:
    print(theta)
print(ll)
```

    4
    [0.352 0.648]
    [0.798 0.202]
    -5756.59565198



```python
import scipy.stats as st

def f(x,y):
    z = np.column_stack([x.ravel(), y.ravel()])
    return (0.1*st.multivariate_normal([0,0], 1*np.eye(2)).pdf(z) +
           0.4*st.multivariate_normal([3,3], 2*np.eye(2)).pdf(z) +
           0.5*st.multivariate_normal([0,5], 3*np.eye(2)).pdf(z))
```


```python
f(np.arange(3), np.arange(3))
```




    array([0.017, 0.012, 0.023])




```python
s = 200
x = np.linspace(-3, 6, s)
y = np.linspace(-3, 8, s)
X, Y = np.meshgrid(x,y)
Z = np.reshape(f(X,Y),(s,s))

from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, cmap='jet')
plt.title('Gaussian Mixure, Model')
```


```python
from scipy.stats import multivariate_normal as mvn
```


```python
def em_gmm_orig(xs, pis, mus, sigmas, tol=0.01, max_iter=100):
    n, p = xs.shape
    k = len(pis)
    
    ll_old = 0
    for i in range(max_iter):
        exp_A = []
        exp_B = []
        ll_new = 0
        
        # E-step
        ws = np.zeros((k, n))
        for j in range(len(mus)):
            for i in range(n):
                ws[j, i] = pis[j] * mvn(mus[j], sigmas[j]).pdf(xs[i])
        ws /= ws.sum(0)
        
        # M-step
        pis = np.zeros(k)
        for j in range(len(mus)):
            for i in range(n):
                pis[j] += ws[j, i]
        pis /= n
        
        mus = np.zeros((k,p))
        for j in range(k):
            for i in range(n):
                mus[j] += ws[j, i] * xs[i]
            mus[j] /= ws[j, :].sum()
        
        sigmas = np.zeros((k,p,p))
        for j in range(k):
            for i in range(n):
                ys = np.reshape(xs[i]-mus[j],(2,1))
                sigmas[j] += ws[j,i] * np.dot(ys, ys.T)
            sigmas[j] /= ws[j,:].sum()
        
        #update complete log likelihood
        ll_new = 0.0
        for i in range(n):
            s = 0
            for j in range(k):
                s += pis[j] * mvn(mus[j], sigmas[j]).pdf(xs[i])
            ll_new += np.log(s)
        
        if np.abs(ll_new - ll_old) < tol:
            break
        ll_old = ll_new
        
    return ll_new, pis, mus, sigmas
```

### EM algorithm (Vector version)


```python
def em_gmm_vect(xs, pis, mus, sigmas, tol=0.01, max_iter=100):
    
    n, p = xs.shape
    k = len(pis)
    
    ll_old = 0
    for i in range(max_iter):
        exp_A = []
        exp_B = []
        ll_new = 0
        
        # E-step
        ws = np.zeros((k,n))
        for j in range(k):
            ws[j,:] = pis[j] * mvn(mus[j], sigmas[j]).pdf(xs)
        ws /= ws.sum(0)
        
        # M-step
        pis = ws.sum(axis=1)
        pis /= n
        
        mus = np.dot(ws, xs)
        mus /= ws.sum(1)[:, None]
        
        sigmas = np.zeros((k, p, p))
        for j in range(k):
            ys = xs - mus[j,:]
            sigmas[j] = (ws[j, :, None, None] * mm(ys[:,:,None], ys[:,None, :])).sum(axis=0)
        sigmas /= ws.sum(axis=1)[:,None,None]
        
        # update complete log likelihood
        ll_new = 0
        for pi, mu, sigma in zip(pis, mus, sigmas):
            ll_new += pi*mvn(mu, sigma).pdf(xs)
        ll_new = np.log(ll_new).sum()
        
        if np.abs(ll_new - ll_old) < tol:
            break
        ll_old = ll_new
        
    return ll_new, pis, mus, sigmas
```

### EM algorithm (Einstein summation)


```python
def em_gmm_eins(xs, pis, mus, sigmas, tol=0.01, max_iter=100):

    n, p = xs.shape
    k = len(pis)

    ll_old = 0
    for i in range(max_iter):
        exp_A = []
        exp_B = []
        ll_new = 0

        # E-step
        ws = np.zeros((k, n))
        for j, (pi, mu, sigma) in enumerate(zip(pis, mus, sigmas)):
            ws[j, :] = pi * mvn(mu, sigma).pdf(xs)
        ws /= ws.sum(0)

        # M-step
        pis = np.einsum('kn->k', ws)/n
        mus = np.einsum('kn,np -> kp', ws, xs)/ws.sum(1)[:, None]
        sigmas = np.einsum('kn,knp,knq -> kpq', ws,
            xs-mus[:,None,:], xs-mus[:,None,:])/ws.sum(axis=1)[:,None,None]

        # update complete log likelihoood
        ll_new = 0
        for pi, mu, sigma in zip(pis, mus, sigmas):
            ll_new += pi*mvn(mu, sigma).pdf(xs)
        ll_new = np.log(ll_new).sum()

        if np.abs(ll_new - ll_old) < tol:
            break
        ll_old = ll_new

    return ll_new, pis, mus, sigmas
```


```python

```


```python
np.random.seed(123)

# create data set
n = 1000
_mus = np.array([[0,4], [-2,0]])
_sigmas = np.array([[[3,0], [0,0.5]],[[1,0],[0,2]]])
_pis = np.array([0.6, 0.4])
xs = np.concatenate([np.random.multivariate_normal(mu, sigma, int(pi*n))
                     for pi, mu, sigma in zip(_pis, _mus, _sigmas)])

# initial guesses for paramters
pis = np.random.random(2)
pis /= pis.sum()
mus = np.random.random((2,2))
sigmas = np.array([np.eye(2)] * 2)
```


```python
%%time
ll1, pis1, mus1, sigmas1 = em_gmm_orig(xs, pis, mus, sigmas)
```

    CPU times: user 9.27 s, sys: 107 ms, total: 9.38 s
    Wall time: 9.4 s



```python
intervals = 101
ys = np.linspace(-8, 8, intervals)
X, Y = np.meshgrid(ys, ys)
_ys = np.vstack([X.ravel(), Y.ravel()]).T

z = np.zeros(len(_ys))
for pi, mu, sigma in zip(pis1, mus1, sigmas1):
    z += pi*mvn(mu, sigma).pdf(_ys)
z = z.reshape((intervals, intervals))

plt.figure(figsize=(10, 8))
ax = plt.subplot(111)
plt.scatter(xs[:,0], xs[:,1], alpha=0.2)
plt.contour(X, Y, z, N=10)
plt.axis([-8, 6, -6,8])
ax.axes.set_aspect('equal')
plt.tight_layout()
```


```python
%%time
ll2, pis2, mus2, sigmas2 = em_gmm_vect(xs, pis, mus, sigmas)
```

    CPU times: user 19.7 ms, sys: 3.68 ms, total: 23.4 ms
    Wall time: 21.1 ms



```python
intervals = 101
ys = np.linspace(-8, 8, intervals)
X, Y = np.meshgrid(ys, ys)
_ys = np.vstack([X.ravel(), Y.ravel()]).T

z = np.zeros(len(_ys))
for pi, mu, sigma in zip(pis2, mus2, sigmas2):
    z += pi*mvn(mu, sigma).pdf(_ys)
z = z.reshape((intervals, intervals))

plt.figure(figsize=(10, 8))
ax = plt.subplot(111)
plt.scatter(xs[:,0], xs[:,1], alpha=0.2)
plt.contour(X, Y, z, N=10)
plt.axis([-8, 6, -6,8])
ax.axes.set_aspect('equal')
plt.tight_layout()
```


```python
%%time
ll3, pis3, mus3, sigmas3 = em_gmm_eins(xs, pis, mus, sigmas)
```

    CPU times: user 18 ms, sys: 1.48 ms, total: 19.4 ms
    Wall time: 20 ms



```python
%timeit em_gmm_orig(xs, pis, mus, sigmas)
%timeit em_gmm_vect(xs, pis, mus, sigmas)
%timeit em_gmm_eins(xs, pis, mus, sigmas)
```

    9.26 s ± 378 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
    17.1 ms ± 6.38 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)
    17.7 ms ± 1.11 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)



```python
intervals = 101
ys = np.linspace(-8,8,intervals)
X, Y = np.meshgrid(ys, ys)
_ys = np.vstack([X.ravel(), Y.ravel()]).T

z = np.zeros(len(_ys))
for pi, mu, sigma in zip(pis3, mus3, sigmas3):
    z += pi*mvn(mu, sigma).pdf(_ys)
z = z.reshape((intervals, intervals))

plt.figure(figsize=(10,8))
ax = plt.subplot(111)
plt.scatter(xs[:,0], xs[:,1], alpha=0.2)
plt.contour(X, Y, z, N=10)
plt.axis([-8,6,-6,8])
ax.axes.set_aspect('equal')
plt.tight_layout()
```


```python

```

b76fe3c660cd7d0e569a346c2ba5ca51ea0e307c
ipynb
Jupyter Notebook
EM_and_VI/EM_Algrithm.ipynb
Gwan-Siu/BlogCode
446871a72b6c8bbfbba20827157910d63490b0eb
MIT
EM_and_VI/EM_Algrithm.ipynb
Gwan-Siu/BlogCode
446871a72b6c8bbfbba20827157910d63490b0eb
MIT
EM_and_VI/EM_Algrithm.ipynb
Gwan-Siu/BlogCode
446871a72b6c8bbfbba20827157910d63490b0eb
MIT
Qwen/Qwen-72B
1. YES
2. YES