Expectation Maximization is a popular generative algorithm, which elegantly handles Multi-Variate Distributions and tries to find close approximations to underlying data-generating distribution.

In following I will consider a particular family of exponential distributions,
the Normal Distribution. (=$N(x | \mu,\Sigma)=\frac{1}{\sqrt{det(2.\pi.\Sigma)}}.e^{\frac{-(x-\mu)^T.\Sigma^{-1}.(x-\mu)}{2}})$
where, $\mu$: is the mean and $\Sigma$ is the covariance matrix.

To fully understand the underlying assumptions behind an algorithm, it usually helps to see how it was derived. Following calculations are a bit tedious with some Aha! moments along the way (and I promise to highlight interesting points!)

Let us consider that we are trying to explain some data distribution using a mix of K Gaussians.
So that is an important assumption right there, that the data is explained by such a mixture of Gaussians.

This derivation is done using Bayesian framework, so first comes the priors.
Say, the prior probability of $k_{th}$ gaussian or the probability that the data is generated by $k_{th}$ gaussian is
$p(z_{k}=1) = \pi_{k}$ … (1)

Now, in such a world of K Gaussians, what is the probability of seeing some data point x.
$p(x | \pi,\mu,\Sigma)=\sum_{k} p(x, z_{k}=1) = \sum_{k} p(x | z_{k}=1) * p(z_{k}=1)$
where,
$p(x | z_{k}=1)$ = probability of seeing data x, in a distribution defined by $k_{th}$ gaussian
= $N(x | \mu_{k},\Sigma_{k})$ defined above.

So, $p(x | \pi,\mu,\Sigma) = \sum_{k=1}^{K} \pi_k * N(x | \mu_{k},\Sigma_{k})$ … (2)

It will also be useful to have reverse conditional handy (basically comes from Bayes Rule)
$\gamma_{k}=p(z_{k}=1|x, \pi,\mu,\Sigma) =\\ \frac{p(x|z_{k}=1).p(z_{k}=1)}{p(x | \pi,\mu\Sigma)} = \frac{\pi_k.N(x | \mu_{k},\Sigma_{k})}{\sum_{k=1}^{K} \pi_k * N(x | \mu_{k},\Sigma_{k})}$ … (3)

Now lets define what we are actually trying to accomplish here
We are basically trying to maximize the likelihood of all data, $X={x_1, x_2, x_3...x_N}$, in a world characterized by some mixture of Gaussians.
That is, $maximize(p(X | \pi,\mu,\Sigma)) = \prod_{n=1}^{N} p(x_n | \pi,\mu,\Sigma) = \prod_{n=1}^{N} (\sum_{k} p(z_{k}=1) * N(x_n | \mu_{k},\Sigma_{k}))$
It is noteworthy, that there is another assumption, that the data points ${x_1, x_2, ... x_N}$ are independent, which is generally not a bad assumption, but something to keep in mind.

maximizing $p(X | \pi,\mu,\Sigma)$ is same as maximizing $ln(p(X | \pi,\mu,\Sigma))$ (as logarithm(x) is a monotonic function of x)
That is, the goal is
$maximize(ln(p(X | \pi,\mu,\Sigma))) = \sum_{n=1}^{N} ln(\sum_{k} p(z_{k}=1) * N(x_n | \mu_{k},\Sigma_{k}))$ … (4)

We have three set of variables,
$\pi={\pi_1, \pi_2...\pi_K}\\ \mu={\mu_1, \mu_2...\mu_K}\\ \Sigma={\Sigma_1, \Sigma_2...\Sigma_K}$
We need to maximize $ln(p(X | \pi,\mu,\Sigma))\ \forall k\in {1,2,...K}$

w.r.t. $\pi_k$
Note $\sum_{k=1}^{K} \pi_k = 1$, which can be incorporated as a constraint using Lagrangian as follows
$L(X | \pi,\mu,\Sigma)=ln(p(X | \pi,\mu,\Sigma)) + \lambda.(\sum_{k=1}^{K} \pi_k - 1) = \sum_{n=1}^{N} ln(\sum_{k=1}^{K} p(z_{k}=1) * N(x_n | \mu_{k},\Sigma_{k})) + \lambda.(\sum_{k=1}^{K} \pi_k - 1)\\ \frac{dL}{d\pi_{k}} = \sum_{n=1}^{N} \frac{N(x_{n} | \mu_{k},\Sigma_{k})}{\sum_{k=1}^{K} \pi_{k}.N(x_n | \mu_{k},\Sigma_{k})} + \lambda = 0\\ \sum_{n=1}^{N} \frac{\pi_{k}.N(x_{n} | \mu_{k},\Sigma_{k})}{\sum_{k=1}^{K} \pi_{k}.N(x_n | \mu_{k},\Sigma_{k})} + \lambda.\pi_{k} = 0\\ Using\ (3),\\ \sum_{n=1}^{N} \gamma_{nk} + \lambda.\pi_{k} = 0\\ \sum_{k=1}^{K} \sum_{n=1}^{N} \gamma_{nk} + \lambda.\sum_{k=1}^{K} \pi_{k} = 0\\ \sum_{n=1}^{N} \sum_{k=1}^{K} \gamma_{nk} + \lambda.\sum_{k=1}^{K} \pi_{k} = 0\\ ==> \lambda=-N\\ ==> \sum_{n=1}^{N} \gamma_{nk} - N.\pi_k=0 ==> \pi_k = \frac{N_k}{N}\ ...\ (5)$ where we have defined $N_k= \sum_{n=1}^{N} \gamma_{nk}$, which can be interpreted as Expected number of data points belonging to Gaussian k.


w.r.t. $\mu$
$\frac{dL}{d\mu_{k}} = \sum_{n=1}^{N} \frac{\pi_{k}.N(x_{n} | \mu_{k},\Sigma_{k})}{\sum_{k=1}^{K} \pi_{k}.N(x_n | \mu_{k},\Sigma_{k})}.(-2.(x_n-\mu_k).\Sigma_k^{-1})=0\\ ==> \sum_{n=1}^{N}\gamma_{nk}.(x_n-\mu_k)=0\\ ==> \mu_k = \frac{\sum_{n=1}^{N}\gamma_{nk}.x_n}{\sum_{n=1}^{N}\gamma_{nk}} = \frac{\sum_{n=1}^{N}\gamma_{nk}.x_n}{N_k}\\$
which basically is the weighted mean of Gaussian k.


w.r.t. $\Sigma$
Now is the hardest part, taking derivative w.r.t. $\Sigma$. Since there are some non-trivial Jacobains involved, I will make use of couple of shortcuts using The Matrix Cook Book
$\frac{dL}{d\Sigma_k} = \sum_{n=1}^{N} \frac{\pi_k.\frac{d(\frac{1}{det(2.\pi.\Sigma_k)}.e^{\frac{-(x-\mu)^T.\Sigma_k^{-1}.(x-\mu)}{2}})}{d\Sigma_k}}{\sum_{k=1}^{K} \pi_{k}.N(x_n | \mu_{k},\Sigma_k)} = 0\\ ==> \frac{\sum_{n=1}^{N} \pi_k.(\frac{-1}{2}.((det(2.\pi.\Sigma_k))^{-\frac{3}{2}}.e^{-\frac{1}{2}.(x_n-\mu_k)^T.\Sigma_k^{-1}.(x_n-\mu_k)}.\frac{d(det(2.\pi.\Sigma_k))}{d\Sigma_k}\ +\ ((det(2.\pi.\Sigma_k))^{-\frac{1}{2}}.e^{-\frac{1}{2}.(x_n-\mu_k)^T.\Sigma_k^{-1}.(x_n-\mu_k)}.\frac{-1}{2}.\frac{d((x-\mu)^T.\Sigma_k^{-1}.(x-\mu))}{d\Sigma_k})}{\sum_{k=1}^{K} \pi_k * N(x | \mu_{k},\Sigma_k)} = 0\\ Now,\ using\\$ $\frac{d(det(2.\pi.\Sigma_k))}{d\Sigma_k} = det(2.\pi.\Sigma_k).\Sigma_k^{-1}\ using\ equation\ 49\ from\ MatrixCookBook\\ and, \frac{d((x-\mu)^T.\Sigma_k^{-1}.(x-\mu))}{d\Sigma_k} = -\Sigma_k^{-1}.(x_n-\mu_k)(x_n-\mu_k)^T.\Sigma_k^{-1}\ using\ equation\ 61\ from\ MatrixCookBook\\ ==> \frac{\sum_{n=1}^{N} \pi_k.N(x | \mu_k,\Sigma_k).(\frac{-1}{2}\Sigma_k^{-1}+\frac{1}{2}.\Sigma_k^{-1}.(x_n-\mu_k)(x_n-\mu_k)^T.\Sigma_k^{-1})}{\sum_{k=1}^{K} \pi_k * N(x | \mu_{k},\Sigma_k)} = 0\\ ==> \sum_{n=1}^{N} \gamma_{nk}.(-I + (x_n-\mu_k)(x_n-\mu_k)^T.\Sigma_k^{-1}) = 0\\ ==> \Sigma_k.\sum_{n=1}^{N} \gamma_{nk} = \sum_{n=1}^{N} \gamma_{nk}.(x_n-\mu_k)(x_n-\mu_k)^T \\ ==> \Sigma_k = \frac{\sum_{n=1}^{N} \gamma_{nk}.(x_n-\mu_k)(x_n-\mu_k)^T}{\sum_{n=1}^{N} \gamma_{nk}} = \frac{\sum_{n=1}^{N} \gamma_{nk}.(x_n-\mu_k)(x_n-\mu_k)^T}{N_k}\\$