The Expectation-Maximization Algorithm: A Method for Modelling Mixtures of Distributions

The expectation-maximization (EM) algorithm provides a method for modelling mixtures of distributions. To explain the EM algorithm, I do so in the context of a coin-flipping example where, for each flip, one of two coins are used, and so a mixture of binomial distributions underlie the data. Although maximum likelihood estimation cannot estimate mixture models, the EM algorithm can because it optimizes the likelihood function indirectly. In the expectation (E) step, a lower-bounding function is used to obtain responsibilities. In the maximization step (M), the lower-bounding function is optimized, with the responsibilities being used to obtain new parameter estimates. By optimizing the lower-bounding function, likelihood function increases by at least as much as the lower-bounding function, thus necessitating another E step. The E and M steps repeat until the parameter values stop updating.

Published
28 April 2023

To write this post on the expectation-maximization algorithm, I used the following academic sources: Citation: Bishop, 2006 Bishop, C. (2006). Pattern recognition and machine learning. Springer New York. Retrieved from bit.ly/411YnEq and Citation: Do & Batzoglou, 2008 Do, C. & Batzoglou, S. (2008). What is the expectation maximization algorithm?. Nature Biotechnology, 26(8), 897–899. https://doi.org/10.1038/nbt1406 . I also used the following lectures/blogs as aids: EM blog and EM lecture.

Introduction to Mixture Models

Consider a situation where a researcher entrusts a colleague to flip coins and record the results of 10 flips. Flips that result in heads are recorded as 1 and flips that result in tails are recorded as 0. Importantly, before each flip, the colleague picks one of two coins (according to another probability mass function) but does not tell the researcher which coin was flipped. When the colleague records the results of the 10 coin flips, they provide the researcher with the following data: $𝐱 = [1, 1, 1, 1, 0, 0, 0, 0, 0, 0]$ . With these data, the researcher wants to estimate the following three parameters:

The probability of picking each of the two coins (i.e., a probability mass function), $𝜇 1, 𝜇 2$ . Given that both these probability values sum to one, only the probability of picking one coin must be estimated, $𝜇 1$ , with $𝜇 2 = 1 - 𝜇 1$ .
The first coin’s probability of heads, $𝑝 1$ .
The second coin’s probability of heads, $𝑝 2$ .

Mixture Models Cannot be Estimated With Maximum Likelihood Estimation

One way to estimate the parameters ( $𝜇 1$ , $𝜇 2$ , $𝑝 1$ , $𝑝 2$ ) is to use maximum likelihood estimation (for a review, see my post on maximum likelihood estimation). In maximum likelihood estimation, we solve for the parameter values $𝜽 = 𝜇 1, 𝑝 1, 𝑝 2$ that maximize the probability of observing the data, $𝑃 (𝐱 | 𝜽)$ . Because the values being maximized are not technically probabilities and are instead likelihoods, $𝐿 (𝜽 | 𝐱)$ , maximum likelihood estimation solves for the parameter values with the highest likelihood, as shown below in Equation $2.1$ :

$𝜽 𝑀 𝐿 𝐸 = a r g m a x 𝜽 𝐿 (𝜽 | 𝐱) . (2.1)$ In the current example, the likelihood values can be calculated by using Equation $2.2$ shown below:

$𝐿 (𝜽 | 𝐱) = 10 \prod 𝑛 = 1 2 \sum 𝑘 = 1 𝜇 𝑘 𝐵 (𝑥 𝑛 | 𝑝 𝑘), (2.2)$

where $𝐵 (𝑥 𝑛 | 𝑝 𝑘)$ is the binomial probability of the $𝑛 𝑡 ℎ$ data point given the the $𝑘 𝑡 ℎ$ coin, and this probability is weighted by the corresponding probability of selecting the $𝑘 𝑡 ℎ$ coin. Importantly, because the researcher does not know which of the two coins produces the result of any flip, then any flip could be the result of flipping the first or second coin. To model this uncertainty, the calculation of the likelihood for each coin flip result, $𝑥 𝑛$ , computes the sum of weighted binomial probabilities, $\sum 2 𝑘 = 1 𝜇 𝑘 𝐵 (𝑥 𝑛 | 𝑝 𝑘)$ . The lack of information surrounding the identity of the coin that produces each flip result also explains why Equation $2.2$ above is often called the incomplete-data likelihood. To prevent underflow (the rounding of small numbers to zero in computers), the log-likelihood is taken, resulting in the incomplete-data log-likelihood shown below in Equation $2.3$ :

$l o g 𝐿 (𝜽 | 𝐱) = 10 \sum 𝑛 = 1 l o g (2 \sum 𝑘 = 1 𝜇 𝑘 𝐵 (𝑥 𝑛 | 𝑝 𝑘)) . (2.3)$

To find maximum likelihood estimates for the parameters, partial derivatives can be computed with respect to each parameter ( $𝜇 1, 𝑝 1, 𝑝 2$ ) and then set to equal zero. In computing the partial derivatives of the incomplete-data log-likelihood (Equation $2.3$ with respect to the parameters, it is important to note that the existence of the summation symbol within the logarithm will oftentimes yield a complex and lengthy derivative because the chain rule has to be applied for each data point. Although computing the partial derivatives does not yield overly complex equations in the current example, solutions (often called closed-form solutions) for the parameters cannot be obtained. As an example, I will show how maximum likelihood estimation would be implemented for estimating the probability of selecting coin 1, $𝜇 1$ . To compute the likelihood, we can expand the binomial term in the incomplete-data log-likelihood function above (Equation $2.3$ ) to produce

$l o g 𝐿 (𝜽 | 𝐱) = 10 \sum 𝑛 = 1 l o g (2 \sum 𝑘 = 1 𝜇 𝑘 (𝑛 ℎ) 𝑝 ℎ 𝑘 (1 - 𝑝 𝑘) 𝑛 - ℎ) . (2.4)$ To allow the partial derivative to be computed, I will apply the binomial calculation on each flip. Thus, $𝑛 = 1$ and $ℎ = {0, 1}$ , which means that $(𝑛 ℎ) = 1$ . In expanding Equation $2.4$ over the summation sign within the logarithm, Equation $2.5$ is obtained

$l o g 𝐿 (𝜽 | 𝐱) = 10 \sum 𝑛 = 1 l o g (𝜇 1 𝑝 ℎ 1 (1 - 𝑝 𝑘) 𝑛 - ℎ + (1 - 𝜇 1) 𝑝 ℎ 2 (1 - 𝑝 2) 𝑛 - ℎ) . (2.5)$ Because $ℎ = {0, 1}$ for any given $𝑛$ flip, the term inside the logarithm will only ever take on one of the two following forms:

$l o g 𝐿 (𝜽 | 𝐱 𝑖) = ⎧ {⎨ {⎩ If 𝐱 𝑖 = ℎ, l o g (𝜇 1 𝑝 1 + (1 - 𝜇 1) 𝑝 2), If 𝐱 𝑖 = 𝑡, l o g (𝜇 1 (1 - 𝑝 1) + (1 - 𝜇 1) (1 - 𝑝 2)) (2.6),$ where $ℎ$ indicates a coin flip that results in ‘heads’ and t indicates a coin flip that results in ’tails’ (i.e., $ℎ = 0$ ). To expand Equation $2.5$ over the summation sign outside the logarithm, we can apply Equation $2.6$ and obtain a simplified expression of the incomplete-data log-likelihood shown below in Equation $2.7$ :

$l o g 𝐿 (𝜽 | 𝐱) = 𝐻 l o g (𝜇 1 𝑝 1 + (1 - 𝜇 1) 𝑝 2) + 𝑇 l o g (𝜇 1 (1 - 𝑝 1) + (1 - 𝜇 1) (1 - 𝑝 2)), (2.7)$ where $𝐻$ indicates the total number of heads and $𝑇$ indicates the total number of tails. In the current data set, $𝐱 = [1, 1, 1, 1, 0, 0, 0, 0, 0, 0]$ , four heads and six tails are obtained. Although we can compute the partial derivative of Equation $2.7$ with respect to $𝜇 1$ and obtain a closed-form solution, Equation $2.8$ shows that it is inadmissible because it always yields negative values for $𝜇 1$ .

$l o g 𝐿 (𝜽 | 𝐱) = 4 l o g (𝜇 1 𝑝 1 + (1 - 𝜇 1) 𝑝 2) + 6 l o g (𝜇 1 (1 - 𝑝 1) + (1 - 𝜇 1) (1 - 𝑝 2)) 𝜕 𝐿 (𝜽 | 𝐱) 𝜕 𝜇 1 = 𝜕 𝜕 𝜇 1 (4 l o g (𝜇 1 𝑝 1 + (1 - 𝜇 1) 𝑝 2) + 6 l o g (𝜇 1 (1 - 𝑝 1) + (1 - 𝜇 1) (1 - 𝑝 2))) = 4 (𝑝 1 - 𝑝 2) 𝜇 1 𝑝 1 + (1 - 𝜇 1) 𝑝 2 + 6 (𝑝 2 - 𝑝 1) 𝜇 1 (1 - 𝑝 1) + (1 - 𝜇 1) (1 - 𝑝 2) Set 𝜕 𝐿 (𝜽 | 𝐱) 𝜕 𝜇 1 = 0 0 = 4 (𝑝 1 - 𝑝 2) 𝜇 1 𝑝 1 + (1 - 𝜇 1) 𝑝 2 + 6 (𝑝 2 - 𝑝 1) 𝜇 1 (1 - 𝑝 1) + (1 - 𝜇 1) (1 - 𝑝 2) = 4 (𝑝 1 - 𝑝 2) 𝜇 1 𝑝 1 - 𝜇 1 𝑝 2 + 𝑝 2 + 6 (𝑝 2 - 𝑝 1) 𝜇 1 - 𝜇 1 𝑝 1 + 1 - 𝜇 1 + 𝜇 1 𝑝 2 - 𝑝 2 = 4 (𝑝 1 - 𝑝 2) 𝜇 1 (𝑝 1 - 𝑝 2) + 𝑝 2 + 6 (𝑝 2 - 𝑝 1) 𝜇 1 (𝑝 2 - 𝑝 1) + 1 - 𝑝 2 = 4 𝜇 1 + 𝑝 2 + 6 𝜇 1 + 1 - 𝑝 2 = 4 (𝜇 1 + 1 - 𝑝 2) + 6 (𝜇 1 + 𝑝 2) = 10 𝜇 1 + 2 𝑝 2 + 4 𝜇 1 = - 𝑝 2 - 2 5 (2.8)$ Therefore, although the summation symbol within the logarithm does not result in an overly complex partial derivative in the current example, maximum likelihood estimation results in inadmissible estimates for parameter values, and so is not a viable method for modelling mixture distributions.¹

Mixture Models Can be Estimated with the Expectation-Maximization (EM) Algorithm

Unlike maximum likelihood estimation, the expectation-maximization (EM) algorithm provides viable parameter estimates for modelling mixture distributions. The EM algorithm works because it indirectly maximizes the incomplete-data log-likelihood; that is, it does not directly operate on the incomplete-data log-likelihood. To act as an indirect estimation method, the EM algorithm begins by modelling the uncertainty of the coin’s identity on each flip as a latent variable: A variable that is assumed to exist but has not been directly measured, whether by choice or because direct measurement is impossible. To model the coin’s identity on each $𝑛$ flip as a latent variable, one-hot (or 1-of-K) encoding is used. In one-hot encoding, the levels of a categorical variable can be represented numerically with a binary vector that sums to one and has a length equal to the number of levels in the categorical variable. In the current example, the categorical variable is the coin’s identity, and the two levels (i.e., coin 1, coin 2) can be represented in each $𝑛$ flip by $𝐳 𝑛$ , as shown below in Equation $3.1$ below:

$𝐳 𝑛 = ⎧ {⎨ {⎩ If coin 1, [1, 0] If coin 2, [0, 1] (3.1) .$ By modelling the coin’s status as a latent variable with one-hot encoding, the incomplete-data likelihood (Equation $2.2$ ) can be modified to produce the complete-data likelihood shown below in Equation $3.2$ :

$𝐿 (𝜽 | 𝐱, 𝐳) = 10 \prod 𝑛 = 1 2 \prod 𝑘 = 1 𝜇 𝑘 𝐵 (𝑥 𝑛 | 𝑝 𝑘) 𝑧 𝑛 𝑘, (3.2)$

where $𝑘$ is used to index each value of $𝐳 𝑛$ . Note that, because the likelihood of each $𝑥 𝑛$ data point is raised to an exponent value of either 0 or 1, we now take the product over the K classes. Thus, data points that do not belong to a mixture do not contribute to the likelihood computation; these likelihoods are raised to the power of zero, and, therefore, contribute a redundant multiplication factor of one.

Although including a binary exponent on the likelihood computation seems inconsequential, it provides two desirable outcomes. First, it allows viable parameter estimates to be computed. In taking the logarithm of the complete-data likelihood to prevent underflow, the complete-data log-likelihood is obtained and is shown below in Equation $3.3$ :

$l o g 𝐿 (𝜽 | 𝐱, 𝐳) = 10 \sum 𝑛 = 1 2 \sum 𝑘 = 1 𝑧 𝑛 𝑘 (l o g (𝜇 𝑘) + l o g (𝐵 (𝑥 𝑛 | 𝑝 𝑘))) . (3.3)$ Because the summation over K is now outside the logarithm, this means that the partial derivatives will be much less complex and yield admissible solutions.

Second, using one-hot encoding creates an equivalence between the complete- and incomplete-data log-likelihoods. As shown below, the complete-data log-likelihood becomes the incomplete-data log-likelihood after summing over all $𝑧$ possible latent variable representations for each $𝑛$ data point (see Equation $3.4$ below). Therefore, we can optimize the incomplete-data log-likelihood by optimizing the complete-data log-likelihood.

$l o g 𝐿 (𝜽 | 𝐱) = l o g (\sum 𝐳 𝐿 (𝜽 | 𝐱, 𝐳) = l o g (\sum 𝐳 10 \prod 𝑛 = 1 2 \prod 𝑘 = 1 𝜇 𝑘 𝐵 (𝑥 𝑛 | 𝑝 𝑘) 𝑧 𝑛 𝑘) = l o g (𝜇 1 𝐵 (𝑥 1 | 𝑝 1) 1 \cdot 𝜇 2 𝐵 (𝑥 1 | 𝑝 2) 0 \cdot \dots \cdot 𝜇 1 𝐵 (𝑥 10 | 𝑝 1) 1 \cdot 𝜇 2 𝐵 (𝑥 10 | 𝑝 2) 0 + 𝜇 1 𝐵 (𝑥 1 | 𝑝 1) 0 \cdot 𝜇 2 𝐵 (𝑥 1 | 𝑝 2) 1 \cdot \dots \cdot 𝜇 1 𝐵 (𝑥 10 | 𝑝 1) 0 \cdot 𝜇 2 𝐵 (𝑥 10 | 𝑝 2) 1) = l o g (10 \prod 𝑛 = 1 2 \sum 𝑘 = 1 𝜇 𝑘 𝐵 (𝑥 𝑛 | 𝑝 𝑘)) = 10 \sum 𝑛 = 1 l o g (2 \sum 𝑘 = 1 𝜇 𝑘 𝐵 (𝑥 𝑛 | 𝑝 𝑘)) . (3.4)$ As an aside, in computing the sum over $𝑧$ , we are marginalizing over $𝑧$ , which explains why the incomplete-data log-likelihood is often called the marginal log-likelihood and the complete-data log-likelihood is often called the joint log-likelihood.

Despite the benefits of using one-hot encoding, one challenge remains. The complete-data log-likelihood cannot yet be computed because the latent variable memberships for each $𝑛$ data point, $𝐳 𝑛$ , are unknown. Recall that the researcher does not know which coin produced the result of any flip. Related, and as an aside, Equation $3.2$ is called the complete-data likelihood because it can only be computed if we know the mixture membership of each $𝑥 𝑛$ data point; that is, we must have the complete data.

Although the complete-data log-likelihood cannot be directly computed, the EM algorithm finds a clever way to circumvent this problem in the E step.

Expectation (E) Step: Using Expectations to Estimate the Distribution of the Latent Variable

The E step takes advantage of the equivalence between the incomplete- and complete-data log-likelihoods and applies two clever tricks to work with the complete-data log-likelihood. First, it applies what appears to be an inconsequential algebraic manipulation, whereby the complete-data log-likelihood is multiplied and divided by some distribution on the latent variable $𝑞 (𝐳)$ as shown below in Equation $3.5$ :

$l o g 𝐿 (𝜽 | 𝐱) = l o g \sum 𝐳 𝐿 (𝜽 | 𝐱, 𝐳) = l o g \sum 𝐳 𝑞 (𝐳) 𝐿 (𝜽 | 𝐱, 𝐳) 𝑞 (𝐳) . (3.5)$

Second, given that the logarithm is a concave function, Jensen’s inequality can be applied to convert Equation $3.5$ into an inequality. Briefly, Jensen’s inequality states that, for concave functions, the function of the expected value, $𝑓 (𝔼 [𝑥])$ ), is greater than or equal to the expected value of the function, $𝔼 [𝑓 (𝑥)]$ (for a proof, see Appendix A), and is shown below in Equation $3.6$ :

$𝑓 (𝔼 [𝑥]) \geq 𝔼 [𝑓 (𝑥)] . (3.6)$

In other words, Equation $3.6$ above shows that $𝔼 [𝑓 (𝑥)]$ is a lower bound on $𝑓 (𝔼 [𝑥])$ . Applying Jensen’s inequality to Equation $3.5$ , a lower bound for the incomplete-data log-likelihood is obtained, $L (𝑞, 𝜽)$ , resulting in the inequality shown below in Equation $3.7$ :

$l o g 𝐿 (𝜽 | 𝐱) \geq \sum 𝐳 𝑞 (𝐳) l o g (𝐿 (𝜽 | 𝐱, 𝐳) 𝑞 (𝐳)) \geq 𝔼 𝑞 (𝐳) l o g (𝐿 (𝜽 | 𝐱, 𝐳) 𝑞 (𝐳)) = L (𝑞, 𝜽), (3.7)$

where the summation over $𝑧$ of $𝑞 (𝐳)$ is the expected value of $𝑞 (𝐳)$ , and can be represented as $𝔼 𝑞 (𝐳)$ . As a note, the lower bound, $L (𝑞, 𝜽)$ , is often called the evidence lower bound because it is a lower bound on the marginal log-likelihood, which is often called the evidence in Bayesian inference.

Although we still do not have a way for computing $𝑞 (𝐳)$ , a closer inspection of the above inequality provides a way forward. In the E step, it is in our best interest to obtain the most accurate approximation of the distribution $𝑞 (𝐳)$ . By obtaining the best estimate of $𝑞 (𝐳)$ , the greatest improvements in the parameter can be realized in the M step. To obtain the best estimate of $𝑞 (𝐳)$ , and thus maximize the potential for improvement in the parameter estimates in the M step, we need to maximize the evidence lower bound with respect to $𝑞 (𝐳)$ . Thus, the inequality of Equation $3.7$ must be transformed into an equality. To do so, we can compute $𝑞 (𝐳)$ such that the logarithm, $l o g (𝐿 (𝜽 | 𝐱, 𝐲) 𝑞 (𝐳))$ , returns constant values, and this can be accomplished if the probability values computed for the latent variable $𝐳$ are proportional to the numerator, $𝑞 (𝐳) \propto 𝐿 (𝜽 | 𝐱, 𝐲)$ . Fortunately, Bayes’ theorem provides one way for us to compute $𝑞 (𝐳)$ to maximize the evidence lower bound such that

where I have used likelihood notation in Equation $3.9$ to highlight the equivalence with the calculation of probabilities and likelihoods. It is important to note that, because latent variable memberships exist for each $𝑛$ data point for each $𝑘$ mixture in the complete-data log-likelihood (Equation $3.3$ ), Equation $3.9$ above is computed for each $𝑛$ data point such that

$𝑃 (𝑧 𝑛 𝑘 | 𝑥 𝑛, 𝜽 = [𝜇 𝑘, 𝑝 𝑘]) = 𝛾 (𝑧 𝑛 𝑘) = 𝜇 𝑘 𝐵 (𝑥 𝑛 | 𝑝 𝑘) \sum 2 𝑘 𝜇 𝑘 𝐵 (𝑥 𝑛 | 𝑝 𝑘) . (3.10)$

Because these values represent the (posterior) probability of membership to each $𝑘$ mixture, they are often called responsibilities. Note that the responsibilities are often represented as $𝛾 (𝑧 𝑛 𝑘)$ , which is simply the scalar form of $𝔼 𝑃 (𝐳 | 𝐱, 𝜽)$ that I use throughout this post. Therefore, by setting $𝑞 (𝐳) = 𝑃 (𝐳 | 𝐱, 𝜽)$ , we can compute $𝑞 (𝐳)$ and also obtain a lower bound, $L (𝑃 (𝐳 | 𝐱, 𝜽), 𝜽)$ that is equal to the complete-data log-likelihood. Using Equation $3.8$ , we can rewrite the inequality of Equation $3.7$ as an equality in Equation $3.12$ below:

To show that the lower bound is equal to the complete-data log-likelihood when $𝑞 (𝐳) = 𝑃 (𝐳 | 𝐱, 𝜽)$ (see Equation $3.11$ ), I provide the Python code block below (see lines 1–109). In order to better understand the Python code, I provide the function for the incomplete-data log-likelihood ( Equation $2.3$ ) and the expansion for the evidence lower bound (Equation $2.3$ ). Importantly, and as I will discuss later on in this post, the first term in Equation $3.14$ is the expected complete-data log-likelihood, and the second term is the entropy of the responsibilities. Recall that the researcher’s data set is $𝐱 = [1, 1, 1, 1, 0, 0, 0, 0, 0, 0]$ .

$l o g 𝐿 (𝜽 | 𝐱) = 10 \sum 𝑛 = 1 l o g (2 \sum 𝑘 = 1 𝜇 𝑘 𝐵 (𝑥 𝑛 | 𝜇 𝑘)) L (𝑃 (𝐳 | 𝐱, 𝜽), 𝜽)) = 𝔼 𝑃 (𝐳 | 𝐱, 𝜽) l o g (𝐿 (𝜽 | 𝐱, 𝐳)) ⏟_⏟_⏟ Expected complete-data log-likelihood 𝑒 𝑥 - 𝔼 𝑃 (𝐳 | 𝐱, 𝜽) l o g (𝑃 (𝐳 | 𝐱, 𝜽)) ⏟_⏟_⏟ Entropy = 10 \sum 𝑛 = 1 2 \sum 𝑘 = 1 𝛾 (𝑧 𝑛 𝑘) (l o g (𝜇 𝑘) + 𝑥 𝑛 l o g (𝑝 𝑘) + (1 - 𝑥 𝑛) l o g (1 - 𝑝 𝑘)) - 𝛾 (𝑧 𝑛 𝑘) l o g (𝛾 (𝑧 𝑛 𝑘)) (2.3 revisited) (3.13) (3.14)$

 import numpy as np
import pandas as pd
from scipy.stats import binom

def e_step(data, mu, p, n = 1):
  """
  Compute expectations (i.e., responsibilities) for each data point's membership to each mixture
  Parameters:
      - data: data set 
      - mu: Probability of each component 
      - p: Probabilities of success for each binomial distribution
  Returns:
      - pandas dataframe
  """
    
  assert len(mu) == len(p), "Number of estimates in mu is equal to the number of sucsess probabilities"
  assert sum(mu) == 1, "Sum of mu should be equal to 1"
  
  #unnormalized responsibilities for each data point for each mixture
  unnormalized_responsibilities = [mu * binom.pmf(x, n=n, p= np.array(p)) for x in data]
  
  #normalized responsibilities (i.e., probabilities)
  normalized_responsibilities = [rp / np.sum(rp) for rp in unnormalized_responsibilities]
  
  column_names = ['resp_mixture_{}'.format(mix+1) for mix in range(len(normalized_responsibilities[0]))]

  df_responsibilities = pd.DataFrame(np.vstack(normalized_responsibilities), 
                                    columns = column_names)
  
  #insert data column as the first one
  df_responsibilities.insert(0, 'data', data)                

  return(df_responsibilities)


#incomplete/complete-data log-likelihood
def compute_incomplete_log_like(data, mu, p):
  """
  Compute incomplete-data log-likelihood 
  Parameters:
      - data: data set 
      - mu: Probability of each component 
      - p: Probability of success for each binomial distribution
  """
  
  #probability of each data point coming from each distribution
  mixture_sums = [np.sum(mu * binom.pmf(flip_result, n=1, p= np.array(p))) for flip_result in binom_mixture_data]
  
  #log of mixture_sums
  log_mixture_sums = np.log(mixture_sums)
  
  #sums of log of mixture_sums
  incomplete_like = np.sum(log_mixture_sums)

  return(incomplete_like)


#lower bound = expected complete-data log-likelihood + entropy 
def compute_lower_bound(responsibilities, mu, p):
  
  #expected complete-data log-likelihood 
  expected_complete_data_like = responsibilities.apply(compute_expected_complete_like, mu = mu, p = p, axis=1).sum()

  ##compute entropy
  entropy = compute_entropy(responsibilities = responsibilities)

  return expected_complete_data_like + entropy


#entropy: sum of rs*log(rs) for all rs (responsibilities)
def compute_entropy(responsibilities):
  
  ##extract responsibility columns and then compute entropy
  resp_colummns = responsibilities.filter(like = 'resp_mixture')
  
  ##take sum of x*log(x) for each responsibility
  entropy = -np.sum(resp_colummns.values * np.log(resp_colummns.values))
  
  return entropy
  

#expected complete-data log-likelihood
def compute_expected_complete_like(row, mu, p):
  resp_columns = [col for col in row.index if 'resp_mixture' in col]
  resp_values = [row[col] for col in resp_columns]
  
  return np.sum(
      [resp_values * (np.log(mu) + 
      row['data'] * np.log(p) + #non-zero if flip result is success (i.e., 'heads')
      (1 - row['data']) * np.log(1 - np.array(p)) #non-zero if flip result is failure (i.e., 'tails')
      )]
  )
    
    
#data given to researcher
binom_mixture_data = [1, 1, 1, 1, 0, 0, 0, 0, 0, 0]

#initial guesses for E step 
mu = [0.3, 0.7] #mixture probabilities 
p = [0.6, 0.8] #success probabilities

#E step
responsibilities = e_step(data = binom_mixture_data, mu = mu, p = p)

evidence_lower_bound = np.round(compute_incomplete_log_like(data = binom_mixture_data, mu = mu, p = p), 5)
incomplete_log_likelihood = np.round(compute_lower_bound(responsibilities = responsibilities,  mu = mu, p = p), 5)

print('Incomplete-data log-likelihood:', evidence_lower_bound)
print('Lower bound:', incomplete_log_likelihood)

 Incomplete-data log-likelihood: -9.28686
Lower bound: -9.28686

Therefore, by introducing a distribution of the latent variable, $𝑞 (𝐳)$ ), its expectation is taken so that the evidence lower bound becomes equal to the incomplete-data log-likelihood. More importantly, however, responsibilities are obtained for each $𝑛$ data point for each $𝑘$ mixture, $𝛾 (𝑧 𝑛 𝑘)$ , which allow new parameter estimates to be obtained in the M step.

Maximization (M) Step: Using Expectations to Obtain New Parameter Estimates

In this section, I first discuss the intuition of computing new parameter estimates with the expectations (see intuition of M step) and then show how these estimates are computed by computing the appropriate partial derivatives (see computation of new estimates).

Understanding the Intuition of the M Step

Beginning with the intuition of the M step, it uses responsibilities obtained in the E step to compute new parameter estimates for $𝜽$ . As in the E step, the M step also indirectly optimizes the incomplete-data log-likelihood by optimizing the evidence lower bound. In the M step, however, instead of optimizing the evidence lower bound with respect to $𝑞 (𝐳)$ , the lower bound is optimized with respect to the parameters, $𝜽$ , resulting in new parameter estimates. Because new parameter estimates are obtained in the M step, I will represent them with $𝜽 𝑛 𝑒 𝑤$ and the old estimates with $𝜽 𝑜 𝑙 𝑑$ . Thus, in optimizing the evidence lower bound with respect to the parameter values in the M step, we can say the lower bound is optimized with respect to $𝜽 𝑜 𝑙 𝑑$ . In optimizing the evidence lower bound with respect to $𝜽 𝑜 𝑙 𝑑$ , the incomplete-data log-likelihood is obtained with new parameter values, $l o g 𝐿 (𝜽 𝑛 𝑒 𝑤 | 𝐱)$ , and increases by at least as much as the lower bound increases when optimized with respect to $𝜽 𝑜 𝑙 𝑑$ , as shown below in Equation $3.15$ :

To understand the logic behind Equation $3.15$ above, a brief discussion of entropy, cross-entropy, and the Kullback-Liebler (KL) divergence is necessary and I provide an overview in the following section.

A Brief Review of Entropy and the Kullback-Liebler (KL) Divergence

Consider an example where an analyst for a wine merchant records wine preferences among customers over a two-year period. Specifically, the analyst asks customers for their favourite type of wine from red, white, and sparkling. Figure 1 shows the customers’ preferences in each year.

Figure 1

Favourite Wine Types by Customers in Each of Two Years

To quantify the difference between the probability distributions with a single value, the analyst uses the Kullback-Liebler (KL) divergence shown below in Equation $3.16$

$𝐾 𝐿 (𝑃 ‖ 𝑄) = \sum 𝐱 𝑃 (𝐱) l o g (𝑃 (𝐱) 𝑄 (𝐱)) . (3.16)$

To understand the KL divergence, it is helpful to understand each of its three computations (for an excellent explanation, see KL divergence) that are presented below:

$𝑃 (𝐱) 𝑄 (𝐱)$ : measures the change in each wine type relative to Year 1.
$l o g$ : gives equal weightings to reciprocals. As an example, consider the change in preferences across the two years for red and white wine. Across the two years, the preference for red wine across increases from 20% to 40%, whereas the preference for white wine decreases from 40% to 20%. Given that these changes are exactly the same, they should contribute the same amount to the total difference between the years. Using logarithm accomplishes this goal; whereas $24 \neq 42$ , $l o g (0.4 0.2) = l o g (0.2 0.4)$ .
$𝑃 (𝐱)$ : each value of $𝐱$ is weighed by its current probability (i.e., Year 2).

Thus, the KL divergence measures the difference between two probability distributions, $𝑃 (𝐱)$ and $𝑄 (𝐱)$ , by computing the sum of weighted logarithmic ratios. Intuitively, if the two distributions are the same, the KL divergence is zero, and if the distributions are different, the KL divergence is positive. Therefore, the KL divergence is always non-negative, $𝐾 𝐿 \geq 0$ (for a proof, see Appendix B).

To understand why the KL divergence is always non-negative, it is important to understand entropy and cross-entropy (for an excellent explanation, see entropy & cross-entropy). If we expand the KL divergence expression of Equation $3.16$ , we obtain

$𝐾 𝐿 (𝑃 ‖ 𝑄) = \sum 𝐱 𝑃 (𝐱) l o g (𝑃 (𝐱)) ⏟_⏟_⏟ (Negative) Entropy - \sum 𝐱 𝑃 (𝐱) l o g (𝑄 (𝐱)) ⏟_⏟_⏟ Cross-entropy . (3.17)$

The first term of Equation $3.17$ represents entropy,² which can be conceptualized as the amount of information or surprise obtained for a given $𝑥$ wine type from the distribution in Year 2, $𝑃 (𝐱)$ , when encoding it by itself.

The second term of Equation $3.17$ represents cross-entropy, which can be conceptualized as the amount of information or surprise obtained for a given $𝑥$ wine type from the distribution in Year 1, $𝑄 (𝐱)$ , when encoded by the distribution in Year 2, $𝑃 (𝐱)$ . Because the distributions in each year are different, it is intuitive to think that cross-entropy is greater than the entropy; that is, it should be more surprising to encode values of one distribution, $𝑃 (𝐱)$ , with values of another distribution, $𝑄 (𝐱)$ , than with values of the same distribution. The conceptualization that cross-entropy is greater than entropy is formally represented by Gibbs’ inequality (for a proof, see Appendix C) below:

$- \sum 𝐱 𝑃 (𝐱) l o g (𝑄 (𝐱)) \geq - \sum 𝐱 𝑃 (𝐱) l o g (𝑃 (𝐱)) . (3.18)$ Using Gibbs’ inequality, it then becomes clear looking at the expression in Equation $3.17$ that the KL divergence will always be non-negative because the larger value of the cross-entropy is added to negative entropy, which has a smaller value.

Computing New Parameter Estimates Increases The Incomplete Log-Likelihood More Than the Evidence Lower Bound

Returning to the M step, the incomplete-data log-likelihood increases by at least as much as the evidence lower bound increases. The inequality between the increase in the incomplete-data log-likelihood and the evidence lower bound after the optimization of the lower bound with respect to $𝜽 𝑜 𝑙 𝑑$ occurs because the optimization results in a larger-value cross-entropy term, $𝔼 𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑) l o g (𝑃 (𝐳 | 𝐱, 𝜽 𝑛 𝑒 𝑤))$ . Importantly, the cross-entropy term is absorbed by the incomplete-data log-likelihood and not the evidence lower bound. To understand how, after the M step, the increase in the incomplete-data log-likelihood is at least the increase in the evidence lower bound, I will prove the following two points:

Point 1: After the M step, the evidence lower bound only increases as much as the expected complete-data log-likelihood.
Point 2: After the M step, the incomplete-data log-likelihood increases by as much as the evidence lower bound and the cross-entropy of the new responsibilities, $𝑃 (𝐳 | 𝐱, 𝜽 𝑛 𝑒 𝑤)$ , with respect to the old responsibilities, $𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑)$ .

Point 1: The Increase in the Evidence Lower Bound is Equal to the Increase in the Expected Complete-Data Log-Likelihood

In computing new parameter estimates in the M step, $𝜽 𝑛 𝑒 𝑤$ , the evidence lower bound increases by the amount that the expected complete-data log-likelihood increases. To show this, I repeat the function for the evidence lower bound below (Equation $3.13$ ) and set $𝜽 = 𝜽 𝑜 𝑙 𝑑$ to keep track of the iteration index.

Below, I show that, when determining the parameters values, $𝜽 𝑜 𝑙 𝑑 = [𝜇 𝑜 𝑙 𝑑 𝑘, 𝑝 𝑜 𝑙 𝑑 𝑘]$ , that maximize the evidence lower bound, the entropy term does not contribute to the derivative (Equation $3.19$ ), and so the maximization is equivalent to maximizing the expected data complete-data log-likelihood (Equation $3.21$ ).

$a r g m a x 𝜽 𝑜 𝑙 𝑑 = [𝜇 𝑜 𝑙 𝑑 𝑘, 𝑝 𝑜 𝑙 𝑑 𝑘] L (𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑), 𝜽 𝑜 𝑙 𝑑)) = a r g m a x 𝜽 𝑜 𝑙 𝑑 = [𝜇 𝑜 𝑙 𝑑 𝑘, 𝑝 𝑜 𝑙 𝑑 𝑘] (10 \sum 𝑛 = 1 2 \sum 𝑘 = 1 𝛾 (𝑧 𝑛 𝑘) (l o g (𝜇 𝑜 𝑙 𝑑 𝑘) + 𝑥 𝑛 l o g (𝑝 𝑜 𝑙 𝑑 𝑘) + (1 - 𝑥 𝑛) l o g (1 - 𝑝 𝑜 𝑙 𝑑 𝑘)) 𝑒 𝑥 - 𝛾 (𝑧 𝑛 𝑘) l o g (𝛾 (𝑧 𝑛 𝑘) ⏟_⏟_⏟ =0 (i.e., constant)) = 10 \sum 𝑛 = 1 2 \sum 𝑘 = 1 𝛾 (𝑧 𝑛 𝑘) (l o g (𝜇 𝑛 𝑒 𝑤 𝑘) + 𝑥 𝑛 l o g (𝑝 𝑛 𝑒 𝑤 𝑘) + (1 - 𝑥 𝑛) l o g (1 - 𝑝 𝑛 𝑒 𝑤 𝑘) a r g m a x 𝜽 𝑜 𝑙 𝑑 = [𝜇 𝑜 𝑙 𝑑 𝑘, 𝑝 𝑜 𝑙 𝑑 𝑘] 𝔼 𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑) l o g (𝐿 (𝜽 𝑜 𝑙 𝑑 | 𝐱, 𝐳)) = 𝔼 𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑) l o g (𝐿 (𝜽 𝑛 𝑒 𝑤 | 𝐱, 𝐳)) (3.19) (3.20) (3.21)$

Therefore, optimizing the evidence lower bound with respect to $𝜽 𝑜 𝑙 𝑑$ is equivalent to maximizing the expected complete-data log-likelihood with respect to $𝜽 𝑜 𝑙 𝑑$ . In other words, to obtain new parameter estimates, we only need to compute partial derivatives with respect to each parameter of the expected complete-data log-likelihood (see computation of estimates for these partial derivatives). Importantly, to compute the value of the lower bound after it has been optimized with respect to $𝜽 𝑜 𝑙 𝑑$ , $L (𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑), 𝜽 𝑛 𝑒 𝑤)$ , the entropy of $𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑)$ is included, as shown below in Equation $3.22$ .

To more concisely represent the difference between the new and old evidence lower bounds, I use auxiliary function notation in Equation $3.23$ , where $𝑄 (𝜽 𝑛 𝑒 𝑤 | 𝜽 𝑜 𝑙 𝑑)$ is the new lower bound and $𝑄 (𝜽 𝑜 𝑙 𝑑 | 𝜽 𝑜 𝑙 𝑑)$ is the old lower bound.

$L (𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑), 𝜽 𝑛 𝑒 𝑤) - L (𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑), 𝜽 𝑜 𝑙 𝑑) = 𝔼 𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑) l o g 𝐿 (𝜽 𝑛 𝑒 𝑤 | 𝐱, 𝐳)) - 𝔼 𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑) l o g 𝐿 (𝜽 𝑜 𝑙 𝑑 | 𝐱, 𝐳)) = 𝑄 (𝜽 𝑛 𝑒 𝑤 | 𝜽 𝑜 𝑙 𝑑) - 𝑄 (𝜽 𝑜 𝑙 𝑑 | 𝜽 𝑜 𝑙 𝑑) (3.23)$ The Python code block below (see lines 112–158) shows that, after the M step, the evidence lower bound indeed only increases by as much as the expected complete-data log-likelihood. Note that, although I have not yet shown how to derive new parameter estimates, I do so in the section on computing new parameter estimates. To better understand how new values are computed for $𝜇 𝑘$ (probability of selecting each $𝑘$ coin) and $𝑝 𝑘$ (probability of heads for each $𝑘$ coin), I provide previews of the solutions in Equations $3.36$ and $3.35$ .

$𝜇 𝑜 𝑙 𝑑 𝑘 = \sum 10 𝑛 = 1 𝛾 (𝑧 𝑛 𝑘) 𝑁 = 𝑁 𝑘 𝑁 = 𝜇 𝑛 𝑒 𝑤 𝑘 𝑝 𝑜 𝑙 𝑑 𝑘 = \sum 10 𝑛 = 1 𝑥 𝑛 𝛾 (𝑧 𝑛 𝑘) \sum 10 𝑛 = 1 𝛾 (𝑧 𝑛 𝑘) = \sum 10 𝑛 = 1 𝑥 𝑛 𝛾 (𝑧 𝑛 𝑘) 𝑁 𝑘 = 𝑝 𝑛 𝑒 𝑤 𝑘 (3.35 𝑠 previewed) (3.36 𝑠 previewed)$

Also note that I fix the values of $𝜇 𝑘$ so that convergence does not occur in one trial.

 def m_step(responsibilities, n = 1):
  
  #isolate columns that contain responsibilities
  resp_cols = responsibilities.filter(like = 'resp_mixture')

  #New estimate for the probability of heads
  """Weigh each data point by the corresponding responsibility and divide by the sum
  of responsibilities for each coin (N_k)."""
  #specify axis=1 so that operations are conducted along rows 
  p_new = np.sum(responsibilities.filter(regex='^resp_mixture').mul(responsibilities['data'], 
  axis=0))/(np.sum(resp_cols)*n)

  #new mixture probabilities 
  mu_new = resp_cols.sum()/resp_cols.sum().sum()

  return pd.DataFrame({'p_new': p_new, 'mu_new': mu_new})


#data given to researcher
binom_mixture_data = [1, 1, 1, 1, 0, 0, 0, 0, 0, 0]


#initial guesses for E step 
mu_fixed = [0.5, 0.5] #mixture probabilities are fixed so that convergence does not occur in one trial
p = [0.6, 0.8] #success probabilities

#E step
responsibilities = e_step(data = binom_mixture_data, mu = mu_fixed, p = p)

#M step 
estimates = m_step(responsibilities = responsibilities)

#amount that lower bound increased by after M step 
optimized_lower_bound = compute_lower_bound(responsibilities = responsibilities, mu = mu_fixed, p = estimates['p_new']) 
expectation_lower_bound = compute_lower_bound(responsibilities = responsibilities, mu = mu_fixed, p = p)

lower_bound_increase = np.round(optimized_lower_bound - expectation_lower_bound, 5)
  

#amount that expected complete-data log-likelihood increases 
optimized_expected_complete = responsibilities.apply(compute_expected_complete_like, mu = mu_fixed, p = estimates['p_new'], axis=1).sum()     
expectation_expected_complete = responsibilities.apply(compute_expected_complete_like, mu = mu_fixed, p = p, axis=1).sum()   

expected_complete_increase =  np.round(optimized_expected_complete - expectation_expected_complete, 5)

print('Increase in evidence lower bound = ', lower_bound_increase)
print('Increase in expected complete-data log-likelihood:', expected_complete_increase)

 Increase in evidence lower bound =  1.81803
Increase in expected complete-data log-likelihood: 1.81803

Point 2: The Incomplete-Data Log-Likelihood Increases by at Least as Much as the Evidence Lower Bound

In computing new parameter estimates in the M step, $𝜽 𝑛 𝑒 𝑤$ , the incomplete-data log-likelihood increases by at least as much as the evidence lower bound increases. To show the inequality between the increase in the incomplete-data log-likelihood and the evidence lower bound, I first show below that the incomplete-data log-likelihood can be decomposed as the sum of the lower bound and a KL divergence (see Equation $3.24$ . Importantly, I use probability notation in the beginning to highlight that $𝐿 (𝜽 | 𝐱, 𝐳) = 𝑃 (𝐱, 𝐳 | 𝜽)$ , which can be decomposed into $𝑃 (𝐳 | 𝐱, 𝜽) 𝑃 (𝐱 | 𝜽) = 𝑃 (𝐳 | 𝐱, 𝜽) 𝐿 (𝜽 | 𝐱)$ . I denote $𝑃 (𝐳 | 𝐱, 𝜽)$ with probability notation because it is a true probability distribution, but I denote $𝑃 (𝐱 | 𝜽)$ with likelihood notation, $𝐿 (𝜽 | 𝐱)$ , because, when fixing the data and varying the parameters, values become likelihoods (see my previous post on likelihood and probability). As mentioned before, because the M step computes new parameter estimates, $𝜽 𝑛 𝑒 𝑤$ , I distinguish them from the current estimates, $𝜽 𝑜 𝑙 𝑑$ .

$L (𝑞, 𝜽 𝑜 𝑙 𝑑) = \sum 𝐳 𝑞 (𝐳) l o g (𝑃 (𝐱, 𝐳 | 𝜽 𝑜 𝑙 𝑑) 𝑞 (𝐳)) = \sum 𝐳 𝑞 (𝐳) l o g (𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑) 𝑃 (𝐱 | 𝜽 𝑜 𝑙 𝑑) 𝑞 (𝐳)) = \sum 𝐳 𝑞 (𝐳) l o g (𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑) 𝑞 (𝐳)) ⏟_⏟_⏟ -ve (reverse) KL divergence + (\sum 𝐳 𝑞 (𝐳) ⏟ = 1 l o g (𝑃 (𝐱 | 𝜽 𝑜 𝑙 𝑑)) = - 𝔼 𝑞 (𝐳) l o g (𝑞 (𝐳) 𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑)) ⏟_⏟_⏟ KL divergence + (l o g (𝑃 (𝐱 | 𝜽 𝑜 𝑙 𝑑)) ⏟_⏟_⏟ Incomplete-data log-likelihood l o g (𝐿 (𝜽 𝑜 𝑙 𝑑 | 𝐱)) = L (𝑞, 𝜽 𝑜 𝑙 𝑑) + 𝔼 𝑞 (𝐳) l o g (𝑞 (𝐳) 𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑)) ⏟_⏟_⏟ KL divergence l o g (𝐿 (𝜽 𝑜 𝑙 𝑑 | 𝐱)) = L (𝑞, 𝜽 𝑜 𝑙 𝑑) + 𝐾 𝐿 (𝑞 (𝐳) ‖ 𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑)) (3.24)$

As a brief aside, we can see that, by setting $𝑞 (𝐳) = 𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑)$ in the E step, the KL divergence goes to zero and the incomplete-data log-likelihood becomes equal to the evidence lower bound (Equation $3.11$ is repeated below, with $𝜽 = 𝜽 𝑜 𝑙 𝑑$ ).

$l o g 𝐿 (𝜽 𝑜 𝑙 𝑑 | 𝐱) = L (𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑), 𝜽 𝑜 𝑙 𝑑) (3.11 revisited)$

After computing new parameter estimates in the M step, the value of the incomplete-data log-likelihood increases at these new values $𝜽$ such that it it the sum of the old evidence lower bound maximized with respect to $𝜽 𝑛 𝑒 𝑤$ and the KL divergence between the old responsibilities and some new distribution of the responsibilities, $𝑞 (𝐳 𝑛 𝑒 𝑤)$ (Equation $3.25$ ). From Equation $3.10$ , we know that we can use Bayes’ theorem to compute the new distribution of the latent variables, and so I set $𝑞 (𝐳 𝑛 𝑒 𝑤 = 𝑃 (𝐳 | 𝐱, 𝜽 𝑛 𝑒 𝑤)$ ) in Equation $3.26$ .

$l o g 𝐿 (𝜽 𝑛 𝑒 𝑤 | 𝐱) = L (𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑), 𝜽 𝑛 𝑒 𝑤) + 𝐾 𝐿 (𝑞 (𝐳 𝑛 𝑒 𝑤) | 𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑)) = L (𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑), 𝜽 𝑛 𝑒 𝑤) + 𝐾 𝐿 (𝑃 (𝐳 | 𝐱, 𝜽 𝑛 𝑒 𝑤) | 𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑)) (3.25) (3.26)$ In looking at the difference between the incomplete-data log-likelihood at the old parameter values with the old responsibilities and the new parameter values, the incomplete-data log-likelihood increases by at least as much as the evidence lower bound because it absorbs the KL divergence between the old responsibilities and the new responsibilities. Therefore, we arrive at Equation $3.15$ presented at the beginning of this section.

Note that Equation $3.15$ can also be represented using auxiliary function and entropy notation, whereby $𝐻 (𝜽 𝑜 𝑙 𝑑 | 𝜽 𝑜 𝑙 𝑑)$ is the entropy of $𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑)$ and $𝐻 (𝜽 𝑛 𝑒 𝑤 | 𝜽 𝑜 𝑙 𝑑)$ is the cross-entropy of $𝑃 (𝐳 | 𝐱, 𝜽 𝑛 𝑒 𝑤)$ with respect to $𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑)$ (Equation $3.30$ ). Because Gibbs’ inequality states that cross-entropy is greater than or equal to entropy (for a proof, see Appendix C), $𝐻 (𝜽 𝑛 𝑒 𝑤 | 𝜽 𝑜 𝑙 𝑑) \geq 𝐻 (𝜽 𝑜 𝑙 𝑑 | 𝜽 𝑜 𝑙 𝑑)$ , and turns Equation $3.29$ into an inequality when removed.

$l o g 𝐿 (𝜽 𝑛 𝑒 𝑤 | 𝐱) - l o g 𝐿 (𝜽 𝑜 𝑙 𝑑 | 𝐱) = 𝑄 (𝜽 𝑛 𝑒 𝑤 | 𝜽 𝑜 𝑙 𝑑) - 𝑄 (𝜽 𝑜 𝑙 𝑑 | 𝜽 𝑜 𝑙 𝑑) + 𝐻 (𝜽 𝑛 𝑒 𝑤 | 𝜽 𝑜 𝑙 𝑑) - 𝐻 (𝜽 𝑜 𝑙 𝑑 | 𝜽 𝑜 𝑙 𝑑) \geq 𝑄 (𝜽 𝑛 𝑒 𝑤 | 𝜽 𝑜 𝑙 𝑑) - 𝑄 (𝜽 𝑜 𝑙 𝑑 | 𝜽 𝑜 𝑙 𝑑) (3.29) (3.30)$ The Python code block below (see lines 161–210) shows that, after the M step, the incomplete-data log-likelihood indeed increases by as much as the evidence lower bound and the KL divergence between the new and old responsibilities. Note that I fix the values of $𝑚 𝑢 𝑘$ so that convergence does not occur in one trial.

 def compute_cross_entropy(old_responsibilities, new_responsibilities):
  
  ##extract responsibility columns and then compute cross-entropy, \sum_x P(x)*log(Q(x)), for each cell value
  resp_colummns_old = old_responsibilities.filter(like = 'resp_mixture')
  resp_colummns_new = new_responsibilities.filter(like = 'resp_mixture')

  entropy = -np.sum(resp_colummns_old.values * np.log(resp_colummns_new.values))
  
  return entropy


#data given to researcher
binom_mixture_data = [1, 1, 1, 1, 0, 0, 0, 0, 0, 0]

#initial guesses for E step 
mu_fixed = [0.5, 0.5] #mixture probabilities 
p = [0.6, 0.8] #success probabilities

#First iteration of EM algorithm 
old_responsibilities = e_step(data = binom_mixture_data, mu = mu_fixed, p = p) #E step
estimates = m_step(responsibilities = old_responsibilities) # M step 

#Second iteration of EM algorithm  (only E step)
new_responsibilities = e_step(data = binom_mixture_data, mu = mu_fixed, p = estimates['p_new'])


#amount that lower bound increased by after M step: increase in lower bound  + KL divergence 
##1) Increase in lower bound 
optimized_lower_bound = compute_lower_bound(responsibilities = old_responsibilities, mu = mu_fixed, p = estimates['p_new']) 
expectation_lower_bound = compute_lower_bound(responsibilities = old_responsibilities, mu = mu_fixed, p = p)
lower_bound_increase = optimized_lower_bound - expectation_lower_bound

##2) KL divergence between new and old responsibilities 
entropy = compute_entropy(responsibilities = old_responsibilities)
cross_entropy = compute_cross_entropy(old_responsibilities = old_responsibilities, 
                                      new_responsibilities = new_responsibilities)
                                      
KL_divergence = -entropy + cross_entropy
inc_lower_bound_plus_KL = np.round(lower_bound_increase + KL_divergence, 5)


#amount that incomplete-data log-likelihood increases 
new_incomplete_like = compute_incomplete_log_like(data = binom_mixture_data, mu = mu_fixed, p = estimates['p_new']) 
old_incomplete_like = compute_incomplete_log_like(data = binom_mixture_data, mu = mu_fixed, p = p)

incomplete_like_increase = np.round(new_incomplete_like - old_incomplete_like, 5)


print('Increase in evidence lower bound + cross-entropy = ', inc_lower_bound_plus_KL)
print('Increase in incomplete-data log-likelihood:', incomplete_like_increase)

 Increase in evidence lower bound + cross-entropy =  1.91468
Increase in incomplete-data log-likelihood: 1.91468

Computation of New Parameter Estimates

Ending with the computation of the new parameter estimates, I will now show how the M step obtains new parameter values, $𝜽 𝑛 𝑒 𝑤$ . As shown in Equation $3.21$ , the new parameter estimates are obtained by maximizing the expected complete-data log-likelihood, $𝔼 𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑) l o g (𝐿 (𝜽 𝑜 𝑙 𝑑 | 𝐱, 𝐳))$ . Therefore, to obtain new parameter values, I will compute the partial derivatives of $𝔼 𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑) l o g (𝐿 (𝜽 𝑜 𝑙 𝑑 | 𝐱, 𝐳))$ with respect to each parameter in $𝜽 𝑜 𝑙 𝑑 = [𝜇 𝑘, 𝑝 𝑘]$ . As a reference, I repeat the expanded version of the expected complete-data log-likelihood in Equation $3.13$ below:

$𝔼 𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑) l o g (𝐿 (𝜽 𝑜 𝑙 𝑑 | 𝐱, 𝐳)) = 10 \sum 𝑛 = 1 2 \sum 𝑘 = 1 𝛾 (𝑧 𝑛 𝑘) (l o g (𝜇 𝑜 𝑙 𝑑 𝑘) + 𝑥 𝑛 l o g (𝑝 𝑜 𝑙 𝑑 𝑘) + (1 - 𝑥 𝑛) l o g (1 - 𝑝 𝑜 𝑙 𝑑 𝑘)) . (3.14 𝑖 revisited)$

For clarity, I have superscripted each parameter with $𝑜 𝑙 𝑑$ because new parameter estimates are always obtained with the old values.

Beginning with the probability of selecting each $𝑘$ coin, $𝜇 𝑘$ , I compute the corresponding partial derivative of the expected complete-data log-likelihood below.

$𝜕 𝔼 𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑) l o g (𝐿 (𝜽 𝑜 𝑙 𝑑 | 𝐱, 𝐳)) 𝜕 𝜇 𝑜 𝑙 𝑑 𝑘 = \sum 10 𝑛 = 1 𝛾 (𝑧 𝑛 𝑘) 𝜇 𝑜 𝑙 𝑑 𝑘 Set partial derivative = 0 0 = \sum 10 𝑛 = 1 𝛾 (𝑧 𝑛 𝑘) 𝜇 𝑜 𝑙 𝑑 𝑘 (3.31)$

In looking at Equation $3.31$ , it is impossible to obtain a solution for $𝜇 𝑘$ . Thankfully, this problem can be fixed by constraining the optimization problem with a Lagrange multiplier, $𝜆$ (for an explanation, see Lagrange multipliers). Specifically, we constrain the $𝑚 𝑢 𝑘$ values to sum to one and then set the derivatives equal to each and include the Lagrange multiplier in the partial derivative of the constraint.

$2 \sum 𝑘 𝜇 𝑜 𝑙 𝑑 𝑘 = 1 2 \sum 𝑘 𝜇 𝑜 𝑙 𝑑 𝑘 - 1 = 0 𝜕 \sum 2 𝑘 𝜇 𝑜 𝑙 𝑑 𝑘 - 1 𝜕 𝜇 𝑜 𝑙 𝑑 𝑘 = 1 \sum 10 𝑛 = 1 𝛾 (𝑧 𝑛 𝑘) 𝜇 𝑜 𝑙 𝑑 𝑘 = 𝜆 (1) Include summation over 𝑘 and solve for 𝜆 \sum 10 𝑛 = 1 \sum 2 𝑘 = 1 𝛾 (𝑧 𝑛 𝑘) \sum 2 𝑘 = 1 𝜇 𝑜 𝑙 𝑑 𝑘 = 𝜆 10 \sum 𝑛 = 1 2 \sum 𝑘 = 1 𝛾 (𝑧 𝑛 𝑘) = 2 \sum 𝑘 = 1 𝜇 𝑜 𝑙 𝑑 𝑘 ⏟ = 1 𝜆 𝜆 = 10 \sum 𝑛 = 1 2 \sum 𝑘 = 1 𝛾 (𝑧 𝑛 𝑘) = 𝑁 (3.32) (3.33) (3.34)$ Note that adding the responsibilities across all $𝑛$ data points and all $𝑘$ mixtures gives the total number of data points, $\sum 10 𝑛 = 1 \sum 2 𝑘 = 1 𝛾 (𝑧 𝑛 𝑘) = 𝑁$ .

We can then plug in the the solution for $𝜆$ (Equation $3.34$ ) into Equation $3.32$ to obtain the solution for $𝜇 𝑘$ shown below in Equation $3.35$ . Intuitively, the probability of each $𝑘$ coin is obtained by dividing the effective number of data points, $𝑁 𝑘$ , by the total number of data points, $𝑁$ .

$\sum 10 𝑛 = 1 𝛾 (𝑧 𝑛 𝑘) 𝜇 𝑘 = 𝑁 𝜇 𝑘 = \sum 10 𝑛 = 1 𝛾 (𝑧 𝑛 𝑘) 𝑁 = 𝑁 𝑘 𝑁 (3.35)$

Ending with the probability of heads for each $𝑘$ coin, $𝑝 𝑘$ , I compute the corresponding partial derivative of the expected complete-data log-likelihood below. As shown in Equation $3.36$ , the value of $𝑝 𝑘$ that optimizes the lower bound and the incomplete-data log-likelihood is obtained by dividing the sum of flips that result in heads weighted by the corresponding responsibility by the effective number of data points associated with the $𝑘$ coin, $𝑁 𝑘$ .

$𝜕 𝔼 𝑃 (𝐳 | 𝐱, 𝜽 𝑜 𝑙 𝑑) l o g (𝐿 (𝜽 𝑜 𝑙 𝑑 | 𝐱, 𝐳)) 𝜕 𝑝 𝑜 𝑙 𝑑 𝑘 = 10 \sum 𝑛 = 1 𝛾 (𝑧 𝑛 𝑘) (𝑥 𝑛 𝑝 𝑜 𝑙 𝑑 𝑘 - 1 - 𝑥 𝑛 1 - 𝑝 𝑜 𝑙 𝑑 𝑘) = 10 \sum 𝑛 = 1 𝛾 (𝑧 𝑛 𝑘) (𝑥 𝑛 (1 - 𝑝 𝑜 𝑙 𝑑 𝑘 - (𝑝 𝑜 𝑙 𝑑 𝑘 (1 - 𝑥 𝑛)) 𝑝 𝑜 𝑙 𝑑 𝑘 (1 - 𝑝 𝑜 𝑙 𝑑 𝑘)) = 10 \sum 𝑛 = 1 𝛾 (𝑧 𝑛 𝑘) (𝑥 𝑛 - 𝑥 𝑛 𝑝 𝑜 𝑙 𝑑 𝑘 - 𝑝 𝑜 𝑙 𝑑 𝑘 + 𝑥 𝑛 𝑝 𝑜 𝑙 𝑑 𝑘) 𝑝 𝑜 𝑙 𝑑 𝑘 (1 - 𝑝 𝑜 𝑙 𝑑 𝑘)) Set partial derivative 𝑠 𝑝 = 0 0 = 10 \sum 𝑛 = 1 𝛾 (𝑧 𝑛 𝑘) (𝑥 𝑛 - 𝑝 𝑜 𝑙 𝑑 𝑘 𝑝 𝑜 𝑙 𝑑 𝑘 (1 - 𝑝 𝑜 𝑙 𝑑 𝑘)) = 10 \sum 𝑛 = 1 𝛾 (𝑧 𝑛 𝑘) (𝑥 𝑛 - 𝑝 𝑜 𝑙 𝑑 𝑘) 𝑝 𝑜 𝑙 𝑑 𝑘 = \sum 10 𝑛 = 1 𝑥 𝑛 𝛾 (𝑧 𝑛 𝑘) \sum 10 𝑛 = 1 𝛾 (𝑧 𝑛 𝑘) = \sum 10 𝑛 = 1 𝑥 𝑛 𝛾 (𝑧 𝑛 𝑘) 𝑁 𝑘 = 𝑝 𝑛 𝑒 𝑤 𝑘 (3.36)$

Conclusion

In conclusion, the expectation-maximization (EM) algorithm provides a method for indirectly optimizing the incomplete-data log-likelihood. To work with the incomplete-data log-likelihood, the EM algorithm uses a lower-bounding function in the E step that can be used to obtain responsibilities. The responsibilities are then used in the M step to optimize the lower-bounding function and obtain new parameter estimates. By optimizing the lower-bounding function, the incomplete-data log-likelihood increases by at least as much as the lower-bounding function, thus necessitating another E step. The E and M steps repeat until the parameter values stop updating.

References

Bishop, C. (2006). Pattern recognition and machine learning. Springer New York. Retrieved from bit.ly/411YnEq

Do, C. & Batzoglou, S. (2008). What is the expectation maximization algorithm?. Nature Biotechnology, 26(8), 897–899. https://doi.org/10.1038/nbt1406

Appendix A: Proof of Jensen’s Inequality

Jensen’s inequality applies to convex and concave functions. Beginning with convex functions, we can use the first-order condition to prove Jensen’s inequality, which is shown below in Equation $A.1$ :

$𝑓 (𝑥) \geq 𝑓 (𝑐) + 𝑓 (𝑐) (𝑥 - 𝑐) . (A.1)$ Given that this definition holds for all $𝑐$ , we can set $𝑐 = 𝔼 [𝑓 (𝑥)]$ , which gives

$𝑓 (𝑥) \geq 𝑓 (𝔼 [𝑥]) + 𝑓 (𝔼 [𝑥]) (𝑥 - 𝔼 [𝑥]) . (A.2)$

We can then take the expectation of Equation $A.2$ , which results in Jensen’s inequality in Equation $A.4$ .

$𝔼 [𝑓 (𝑥)] \geq 𝔼 [𝑓 (𝔼 [𝑥])] ⏟ = 𝑓 (𝔼 [𝑥] + 𝔼 [𝑓 (𝔼 [𝑥])] (𝔼 [𝑥] - 𝔼 [𝑥]) ⏟_⏟_⏟ = 0 . 𝔼 [𝑓 (𝑥)] \geq 𝑓 (𝔼 [𝑥]) (A.3) (A.4)$

For concave functions, we simply change the direction of the inequality in Equation $A.1$ to obtain Jensen’s inequality for concave functions below (Equation $A.5$ ).

$𝑓 (𝔼 [𝑥]) \geq 𝔼 [𝑓 (𝑥)] (A.5)$

Appendix B: Proof that KL Divergence is Always Non-Negative

As explained in section on the KL divergence and entropy, the KL divergence measures the difference between two probability distributions, $𝑃 (𝐱)$ and $𝑄 (𝐱)$ , by computing the sum of weighted logarithmic ratios, as shown below in a repetition of Equation $3.16$ . Note that I flip the fraction inside the logarithm to obtain the negative (or reverse) KL divergence in Equation $B.1$

$𝐾 𝐿 (𝑃 ‖ 𝑄) = \sum 𝐱 𝑃 (𝐱) l o g (𝑃 (𝐱) 𝑄 (𝐱)) . - 𝐾 𝐿 (𝑃 ‖ 𝑄) = - \sum 𝐱 𝑃 (𝐱) l o g (𝑄 (𝐱) 𝑃 (𝐱)) . (3.16 𝑠 revisited) (B.1)$

Because the negative logarithm causes the KL divergence to become convex, we can use Jensen’s inequality to transform the KL divergence into an inequality by moving the summation inside the logarithm (Equation $B.2$ ). From this point onward, it becomes clear that the KL diverence is always greater than or equal to zero (Equation $B.3$ ).

$- 𝐾 𝐿 (𝑃 ‖ 𝑄) \geq - l o g (\sum 𝐱 𝑃 (𝐱) 𝑄 (𝐱) 𝑃 (𝐱)) \geq l o g (𝑄 (𝐱) ⏟ = 1) \geq l o g (1) \cdot - 1 𝐾 𝐿 (𝑃 ‖ 𝑄) \geq 0 (B.2) (B.3)$

Appendix C: Proof of Gibbs’ Inequality

Gibbs’ inequality states that the cross-entropy, $- 𝑃 (𝐱) l o g (𝑄 (𝐱))$ , is greater than or equal to the entropy, $- 𝑃 (𝐱) l o g (𝑃 (𝐱))$ (see Equation $C.1$ ). Gibbs’ inequality can be easily proven using the fact that the KL divergence is always greater than or equal to zero.

$𝐾 𝐿 (𝑃 ‖ 𝑄) \geq 0 \sum 𝐱 𝑃 (𝐱) l o g (𝑃 (𝐱)) - \sum 𝐱 𝑃 (𝐱) l o g (𝑄 (𝐱)) \geq 0 - \sum 𝐱 𝑃 (𝐱) l o g (𝑄 (𝐱)) \geq - \sum 𝐱 𝑃 (𝐱) l o g (𝑄 (𝐱)) (C.1)$

It should also be noted that maximum likelihood estimation can result in singularities with Gaussian mixture models. In short, if the estimate for a $𝑘 𝑡 ℎ$ mixture’s mean, $𝜇 𝑘$ , happens to exactly match the value of an individual $𝑛$ value, $𝑥 𝑛$ , then the mixture in question can become ‘stuck’ on this data point with all the other data points being modelled by the other mixtures. With the mixture fixed on the one data point, the variance of the mixture will decrease to zero and the log-likelihood will increase to infinity. ↩︎
The first value in Equation $3.17$ is technically the negative entropy. Because entropy (and cross-entropy) compute information/surprise, it makes conceptual sense to represent them as positive values. Unfortunately, the term $\sum 𝐱 𝑃 (𝐱) l o g (𝑃 (𝐱))$ returns negative values. To reflect the conceptualization that entropy computes information, a negative sign is included to multiply the negative value returned by $\sum 𝐱 𝑃 (𝐱) l o g (𝑃 (𝐱))$ into a positive value. ↩︎

	import numpy as np
	import pandas as pd
	from scipy.stats import binom

	def e_step(data, mu, p, n = 1):
	"""
	Compute expectations (i.e., responsibilities) for each data point's membership to each mixture
	Parameters:
	- data: data set
	- mu: Probability of each component
	- p: Probabilities of success for each binomial distribution
	Returns:
	- pandas dataframe
	"""

	assert len(mu) == len(p), "Number of estimates in mu is equal to the number of sucsess probabilities"
	assert sum(mu) == 1, "Sum of mu should be equal to 1"

	#unnormalized responsibilities for each data point for each mixture
	unnormalized_responsibilities = [mu * binom.pmf(x, n=n, p= np.array(p)) for x in data]

	#normalized responsibilities (i.e., probabilities)
	normalized_responsibilities = [rp / np.sum(rp) for rp in unnormalized_responsibilities]

	column_names = ['resp_mixture_{}'.format(mix+1) for mix in range(len(normalized_responsibilities[0]))]

	df_responsibilities = pd.DataFrame(np.vstack(normalized_responsibilities),
	columns = column_names)

	#insert data column as the first one
	df_responsibilities.insert(0, 'data', data)

	return(df_responsibilities)


	#incomplete/complete-data log-likelihood
	def compute_incomplete_log_like(data, mu, p):
	"""
	Compute incomplete-data log-likelihood
	Parameters:
	- data: data set
	- mu: Probability of each component
	- p: Probability of success for each binomial distribution
	"""

	#probability of each data point coming from each distribution
	mixture_sums = [np.sum(mu * binom.pmf(flip_result, n=1, p= np.array(p))) for flip_result in binom_mixture_data]

	#log of mixture_sums
	log_mixture_sums = np.log(mixture_sums)

	#sums of log of mixture_sums
	incomplete_like = np.sum(log_mixture_sums)

	return(incomplete_like)


	#lower bound = expected complete-data log-likelihood + entropy
	def compute_lower_bound(responsibilities, mu, p):

	#expected complete-data log-likelihood
	expected_complete_data_like = responsibilities.apply(compute_expected_complete_like, mu = mu, p = p, axis=1).sum()

	##compute entropy
	entropy = compute_entropy(responsibilities = responsibilities)

	return expected_complete_data_like + entropy


	#entropy: sum of rs*log(rs) for all rs (responsibilities)
	def compute_entropy(responsibilities):

	##extract responsibility columns and then compute entropy
	resp_colummns = responsibilities.filter(like = 'resp_mixture')

	##take sum of x*log(x) for each responsibility
	entropy = -np.sum(resp_colummns.values * np.log(resp_colummns.values))

	return entropy


	#expected complete-data log-likelihood
	def compute_expected_complete_like(row, mu, p):
	resp_columns = [col for col in row.index if 'resp_mixture' in col]
	resp_values = [row[col] for col in resp_columns]

	return np.sum(
	[resp_values * (np.log(mu) +
	row['data'] * np.log(p) + #non-zero if flip result is success (i.e., 'heads')
	(1 - row['data']) * np.log(1 - np.array(p)) #non-zero if flip result is failure (i.e., 'tails')
	)]
	)


	#data given to researcher
	binom_mixture_data = [1, 1, 1, 1, 0, 0, 0, 0, 0, 0]

	#initial guesses for E step
	mu = [0.3, 0.7] #mixture probabilities
	p = [0.6, 0.8] #success probabilities

	#E step
	responsibilities = e_step(data = binom_mixture_data, mu = mu, p = p)

	evidence_lower_bound = np.round(compute_incomplete_log_like(data = binom_mixture_data, mu = mu, p = p), 5)
	incomplete_log_likelihood = np.round(compute_lower_bound(responsibilities = responsibilities, mu = mu, p = p), 5)

	print('Incomplete-data log-likelihood:', evidence_lower_bound)
	print('Lower bound:', incomplete_log_likelihood)

	Incomplete-data log-likelihood: -9.28686
	Lower bound: -9.28686

	def m_step(responsibilities, n = 1):

	#isolate columns that contain responsibilities
	resp_cols = responsibilities.filter(like = 'resp_mixture')

	#New estimate for the probability of heads
	"""Weigh each data point by the corresponding responsibility and divide by the sum
	of responsibilities for each coin (N_k)."""
	#specify axis=1 so that operations are conducted along rows
	p_new = np.sum(responsibilities.filter(regex='^resp_mixture').mul(responsibilities['data'],
	axis=0))/(np.sum(resp_cols)*n)

	#new mixture probabilities
	mu_new = resp_cols.sum()/resp_cols.sum().sum()

	return pd.DataFrame({'p_new': p_new, 'mu_new': mu_new})


	#data given to researcher
	binom_mixture_data = [1, 1, 1, 1, 0, 0, 0, 0, 0, 0]


	#initial guesses for E step
	mu_fixed = [0.5, 0.5] #mixture probabilities are fixed so that convergence does not occur in one trial
	p = [0.6, 0.8] #success probabilities

	#E step
	responsibilities = e_step(data = binom_mixture_data, mu = mu_fixed, p = p)

	#M step
	estimates = m_step(responsibilities = responsibilities)

	#amount that lower bound increased by after M step
	optimized_lower_bound = compute_lower_bound(responsibilities = responsibilities, mu = mu_fixed, p = estimates['p_new'])
	expectation_lower_bound = compute_lower_bound(responsibilities = responsibilities, mu = mu_fixed, p = p)

	lower_bound_increase = np.round(optimized_lower_bound - expectation_lower_bound, 5)


	#amount that expected complete-data log-likelihood increases
	optimized_expected_complete = responsibilities.apply(compute_expected_complete_like, mu = mu_fixed, p = estimates['p_new'], axis=1).sum()
	expectation_expected_complete = responsibilities.apply(compute_expected_complete_like, mu = mu_fixed, p = p, axis=1).sum()

	expected_complete_increase = np.round(optimized_expected_complete - expectation_expected_complete, 5)

	print('Increase in evidence lower bound = ', lower_bound_increase)
	print('Increase in expected complete-data log-likelihood:', expected_complete_increase)

	Increase in evidence lower bound = 1.81803
	Increase in expected complete-data log-likelihood: 1.81803

	def compute_cross_entropy(old_responsibilities, new_responsibilities):

	##extract responsibility columns and then compute cross-entropy, \sum_x P(x)*log(Q(x)), for each cell value
	resp_colummns_old = old_responsibilities.filter(like = 'resp_mixture')
	resp_colummns_new = new_responsibilities.filter(like = 'resp_mixture')

	entropy = -np.sum(resp_colummns_old.values * np.log(resp_colummns_new.values))

	return entropy


	#data given to researcher
	binom_mixture_data = [1, 1, 1, 1, 0, 0, 0, 0, 0, 0]

	#initial guesses for E step
	mu_fixed = [0.5, 0.5] #mixture probabilities
	p = [0.6, 0.8] #success probabilities

	#First iteration of EM algorithm
	old_responsibilities = e_step(data = binom_mixture_data, mu = mu_fixed, p = p) #E step
	estimates = m_step(responsibilities = old_responsibilities) # M step

	#Second iteration of EM algorithm (only E step)
	new_responsibilities = e_step(data = binom_mixture_data, mu = mu_fixed, p = estimates['p_new'])


	#amount that lower bound increased by after M step: increase in lower bound + KL divergence
	##1) Increase in lower bound
	optimized_lower_bound = compute_lower_bound(responsibilities = old_responsibilities, mu = mu_fixed, p = estimates['p_new'])
	expectation_lower_bound = compute_lower_bound(responsibilities = old_responsibilities, mu = mu_fixed, p = p)
	lower_bound_increase = optimized_lower_bound - expectation_lower_bound

	##2) KL divergence between new and old responsibilities
	entropy = compute_entropy(responsibilities = old_responsibilities)
	cross_entropy = compute_cross_entropy(old_responsibilities = old_responsibilities,
	new_responsibilities = new_responsibilities)

	KL_divergence = -entropy + cross_entropy
	inc_lower_bound_plus_KL = np.round(lower_bound_increase + KL_divergence, 5)


	#amount that incomplete-data log-likelihood increases
	new_incomplete_like = compute_incomplete_log_like(data = binom_mixture_data, mu = mu_fixed, p = estimates['p_new'])
	old_incomplete_like = compute_incomplete_log_like(data = binom_mixture_data, mu = mu_fixed, p = p)

	incomplete_like_increase = np.round(new_incomplete_like - old_incomplete_like, 5)


	print('Increase in evidence lower bound + cross-entropy = ', inc_lower_bound_plus_KL)
	print('Increase in incomplete-data log-likelihood:', incomplete_like_increase)

	Increase in evidence lower bound + cross-entropy = 1.91468
	Increase in incomplete-data log-likelihood: 1.91468

Whitepapers Directory