Using EM To Estimate A Probability Density With A Mixture Of Gaussians

Aaron A. D’Souza

adsouza@usc.edu

1 Introduction

The problem we are trying to address in this note is simple. Given a set of data points X = {x_{i}}_{1}^{N}, we wish to determine the underlying probability distribution p(x), that generates this data. In general, the distribution p(x) could be any real-valued, scalar function with the following constraints:

p(x) ≥ 0, ∀x and ∫_{−∞}^{∞} p(x)dx = 1

The framework of mixture modeling using Gaussians makes the following assumptions:

• The data was generated using a set of K probability distributions.

• Each of the individual probability distributions is a Gaussian:

x; θ_{k} ∼ Normal (x; μ_{k}, Σ_{k})^{∗}

i.e., the probability of generating a data point x under the m^{th} model is given according to a Gaussian distribution with mean μ_{k} and covariance Σ_{k}.

• Each data point is generated according to the following algorithm:

1: for i = 1 to N do

2: m ← index of one of the K models randomly selected

according to the prior probability vector π 3: Randomly generate x_{i} according to the distribution Normal (x_{i}; μ_{k}, Σ_{k}) 4: end for

Representing probabilistic systems as graphical models is rapidly becoming a useful tool in Bayesian analysis. The graphical model corresponding to our formulation of the data generation process is shown in fig. 1. We adopt the convention that circular nodes correspond to random variables in the model, while rectangular nodes correspond to variables that parameterize the distributions of these random variables. Circular nodes with double borders indicate observed variables in our model (the data).

∗We use the notation Normal (x; μ, Σ) to denote the multivariate Gaussian distribution which is mathematically defined as:

Normal (x; μ, Σ) = (2π)^{−d/2} |Σ|^{−1/2} exp \{− \frac{1}{2} (x − μ)^{T} Σ^{−1}(x − μ)\} where d is the dimensionality of x.

μ_{m}  π x_{i}  z_{i}^{Σ_{m}}

Figure 1: Graphical model for maximum likelihood density estimation using a mixture of Gaussians

In this model we have introduced an additional variable z_{i} associated with each x_{i}. The z_{i} variables are indicator variables that are multinomially distributed according to the parameter vector π, and indicate which component generates the corresponding x_{i}. It is easiest to think of each z_{i} as an K dimensional vector with a 1 in the element corresponding to the selected mixture component, and 0’s in all other elements. The probability of the m^{th} element being 1, is w_{k}.

z_{i} = \underbrace{[0 0 ⋯ 0 1 0 ⋯ 0 0]^{T}}_{K elements}

Since we do not know the corresponding z_{i} for each x_{i} (if we did, then we simply group the x_{i} according to their z_{i}, and fit a single Gaussian to each group), these variables are called hidden variables.

Our problem is now reduced to finding the values of the model parameters μ_{k} and Σ_{k} for each of the K models, as well as the prior probability vector π, which when plugged into the generative model, is most likely to generate the observed data distribution. In other words we are interested in maximizing the likelihood L(θ) = p(X; θ) of generating the observed data given the model parameters θ = {μ_{k}, Σ_{k}, w_{k}}_{1}^{K}. This approach is called the Maximum Likelihood (ML) framework since it finds the parameter settings that maximize the likelihood of observing the data.

Although the ML approach is an intuitively appealing solution, we often find that maximizing the expressions for likelihood w.r.t. the parameters θ are often analytically intractable. The Expectation Maximization (EM) algorithm can be used to simplify the math considerably.

Figure 2: Mixtures equations

2 Estimating the Model Parameters using EM

Instead of attempting to maximize the likelihood of the observed data p(X; θ), we attempt instead to maximize the likelihood of the joint distribution of X and Z = {z_{i}}_{1}^{N}, p(X, Z; θ). For the purposes of maximization we can also work with the logarithm of this quantity, log p(X, Z; θ). This quantity is also known as the complete log-likelihood. Since we cannot observe the values of the random variables z_{i} we must work with the expectation of this quantity w.r.t. some distribution Q(Z).

The log of the complete data likelihood can be written as follows:

l_{c}(θ) = log p(X, Z; θ)

= log ∏_{i}^{N} p(x_{i}, z_{i}; θ)

= log ∏_{i}^{N} ∏_{k}^{K} [p(x_{i}|π_{ik} = 1; θ)p(π_{ik} = 1)]^{π_{ik}}

= ∑_{i=1}^{N} ∑_{k}^{K} π_{ik} log p(x_{i}|π_{ik} = 1; θ) + π_{ik} log w_{k}

Since we have assumed that each of the individual models is a Gaussian, the quantity p(x_{i}|m, θ) is simply the conditional probability of generating x_{i} given that the m^{th} model is chosen:

log p(x_{i}|π_{ik} = 1; θ) = \frac{1}{(2π)^{d/2} |Σ_{k}|^{1/2}} exp \{−\frac{1}{2} (x_{i} − μ_{k})^{T} Σ_{k}^{−1} (x_{i} − μ_{k})\}  (1)

Taking expectations w.r.t. Q(Z) we get:

⟨l_{c}(θ)⟩_{Q(Z)} = ∑_{i=1}^{N} ∑_{k}^{K} ⟨π_{ik}⟩ log p(x_{i}|π_{ik} = 1; θ) + ⟨π_{ik}⟩ log w_{k}  (2)

2.1 The M step

The “M” step in EM takes the expected complete log-likelihood as defined in eq. (2) and maximizes it w.r.t. the parameters that are to be estimated; in this case w_{k}, μ_{k}, and Σ_{k}.

Differentiating eq. (2) w.r.t. μ_{k} we get:

\frac{∂ ⟨l_{c}(θ)⟩_{Q(Z)}}{∂μ_{k}} = ∑_{i=1}^{N} ⟨π_{ik}⟩ \frac{∂}{∂μ_{k}} log p(x_{i}|π_{ik} = 1; θ) = 0  (3)

We can compute \frac{∂}{∂μ_{k}} log p(x_{i}|π_{ik} = 1; θ) using eq. (1) as follows:

\frac{∂}{∂μ_{k}} log p(x_{i}|π_{ik} = 1; θ) = \frac{∂}{∂μ_{k}} log \{\frac{1}{(2π)^{d/2} |Σ_{k}|^{1/2}} exp \{−\frac{1}{2} (x_{i} − μ_{k})^{T} Σ_{k}^{−1} (x_{i} − μ_{k})\}\}

= −\frac{1}{2} \frac{∂}{∂μ_{k}} (x_{i} − μ_{k})^{T} Σ_{k}^{−1} (x_{i} − μ_{k})

= (x_{i} − μ_{k})^{T} Σ_{k}^{−1}^{†}

Substituting this result into eq. (3), we get:

∑_{i=1}^{N} ⟨π_{ik}⟩ (x_{i} − μ_{k})^{T} Σ_{k}^{−1} = 0

giving us the update equation:

μ_{k} = \frac{∑_{i=1}^{N} ⟨π_{ik}⟩ x_{i}}{∑_{i=1}^{N} ⟨π_{ik}⟩}  (4)

Differentiating eq. (2) w.r.t. Σ_{k}^{−1} we get:

\frac{∂ ⟨l_{c}(θ)⟩_{Q(Z)}}{∂Σ_{k}^{−1}} = ∑_{i=1}^{N} ⟨π_{ik}⟩ \frac{∂}{∂Σ_{k}^{−1}} log p(x_{i}|π_{ik} = 1; θ) = 0  (5)

†Where we have used the relation \frac{∂}{∂x} x^{T} Ax = x^{T} (A + A^{T})

We can compute _{\frac{∂}{∂Σ_{k}^{−1}}} log p(x_{i}|π_{ik} = 1; θ) using eq. (1) as follows:

\frac{∂}{∂Σ_{k}^{−1}} log p(x_{i}|π_{ik} = 1; θ) = \frac{∂}{∂Σ_{k}^{−1}} log \{\frac{1}{(2π)^{d/2} |Σ_{k}|^{1/2}} exp \{−\frac{1}{2} (x_{i} − μ_{k})^{T} Σ_{k}^{−1} (x_{i} − μ_{k})\}\}

= \frac{∂}{∂Σ_{k}^{−1}} \{\frac{1}{2} log |Σ_{k}^{−1}| − \frac{1}{2} (x_{i} − μ_{k})^{T} Σ_{k}^{−1} (x_{i} − μ_{k})\}

= \frac{1}{2}Σ_{k} − \frac{1}{2} (x_{i} − μ_{k}) (x_{i} − μ_{k})^{T} ^{‡}

Substituting this result into eq. (5), we get:

∑_{i=1}^{N} ⟨π_{ik}⟩ (\frac{1}{2}Σ_{k} − \frac{1}{2} (x_{i} − μ_{k}) (x_{i} − μ_{k})^{T} ) = 0

giving us the update equation:

Σ_{k} = \begin{matrix}∑ & {N \atop i=1} & ⟨π_{ik}⟩ & (x_{i} − μ_{k}) & (x_{i} − μ_{k}) & T \\ ∑ & {N \atop i=1} & ⟨π_{ik}⟩ \end{matrix}  (6)

In order to maximize the expected log-likelihood in eq. (2) w.r.t. w_{k}, we have to keep in mind that the maximization has the constraint that ∑_{k}^{K} w_{k} = 1. In order to enforce this constraint we use the Lagrange multiplier λ, and augment eq. (2) as follows:

L′ (θ) = ⟨l_{c}(θ)⟩_{Q(Z)} − λ (∑_{k}^{K} w_{k} − 1)  (7)

We now differentiate this new expression w.r.t. each w_{k} giving us:

\frac{∂}{∂w_{k}} ⟨l_{c}(θ)⟩_{Q(Z)} − λ = 0 for 1 ≤ m ≤ K

Using eq. (2) we get:

\begin{matrix}\frac{1}{w_{k}} ∑_{i=1}^{N} ⟨π_{ik}⟩ − λ & = 0 \\ or equivalently ∑_{i=1}^{N} ⟨π_{ik}⟩ − λw_{k} & = 0 \end{matrix} \}  for 1 ≤ m ≤ K  (8)

Summing eq. (8) over all K models we get:

∑_{k}^{K} ∑_{i=1}^{N} ⟨π_{ik}⟩ − λ ∑_{k}^{K} w_{k} = 0

But since ∑_{k}^{K} w_{k} = 1 we get:

λ = ∑_{k}^{K} ∑_{i=1}^{N} ⟨π_{ik}⟩ = N  (9)

‡Where we have used the relation \frac{∂}{∂X} log |X| = (X^{−1})^{T} and \frac{∂}{∂A} x^{T} Ax = xx^{T}

Substituting this result back into eq. (8) we get the following update equation:

w_{k} = \frac{∑_{i=1}^{N} ⟨π_{ik}⟩}{N}  (10)

which preserves the constraint that ∑_{k}^{K} w_{k} = 1.

2.2 The E step

Now that we have derived the update equations that maximize the expected complete loglikelihood ⟨log p(X, Z; θ)⟩, we wish to ensure that we are indeed also maximizing the incomplete log-likelihood p(X; θ) (which is the quantity that we are truly interested in maximizing).

As we mentioned earlier in section 2, we are guaranteed to maximize the incomplete loglikelihood only when the expectation is taken w.r.t. the posterior distribution of Z, namely p(Z|X; θ). Hence each of the expectations ⟨π_{ik}⟩ that appear in the update equations derived in the previous section (section 2.1), should be computed as follows:

⟨π_{ik}⟩_{p(Z|X;θ)} = 1 · p(π_{ik} = 1|x_{i}; θ) + 0 · p(π_{ik} = 0|x_{i}; θ)

= p(π_{ik} = 1|x_{i}; θ)

= \frac{p(x_{i}|π_{ik} = 1; θ)p(π_{ik} = 1)}{∑_{k=1}^{K} p(x_{i}|z_{ik} = 1; θ)p(z_{ik} = 1)}

= \frac{p(x_{i}|π_{ik} = 1; θ)w_{k}}{∑_{k=1}^{K} p(x_{i}|π_{ik} = 1; θ)w_{k}}

= \frac{p(x_{i}|π_{ik} = 1; θ)w_{k}}{notmalized over k}

3 Summary

To summarize, given a set of data points X, if we wish to estimate the underlying probability distribution using EM to fit K Gaussians, we apply Algorithm 1, iterating until convergence of the model parameters.

4 Why does this work? (A brief review of EM theory)

Knowing that log p(X; θ) is difficult to maximize analytically, we (seemingly arbitrarily) chose to maximize the expected complete log-likelihood ⟨log p(X, Z; θ)⟩_{Q(Z)} in the hope that this also increases the incomplete log-likelihood log p(X; θ) (the quantity we are really interested in). This section will justify this choice and prove that we are indeed maximizing log p(X; θ)

Let us rewrite log p(X|θ) as follows:

Initialize: all ⟨π_{ik}⟩, w_{k}, μ_{k}, and Σ_{k} 1: repeat 2: for i = 1 to N do //The E step 3: for m = 1 to K do

4:

p(x_{i}|π_{ik} = 1; θ) = (2π)^{−d/2} |Σ_{k}|^{−1/2} exp \{−\frac{1}{2}(x_{i} − μ_{k})^{T} Σ_{k}^{−1}(x_{i} − μ_{k})\}
⟨π_{ik}⟩ = \frac{p(x_{i}|π_{ik} = 1; θ)w_{k}}{normalized over k}

5: end for 6: end for 7: for m = 1 to K do //The M step

8:

Σ_{k} = \begin{matrix}∑ & {N \atop i=1} & ⟨π_{ik}⟩ & (x_{i} − μ_{k}) & (x_{i} − μ_{k}) & T \\ ∑ & {N \atop i=1} & ⟨π_{ik}⟩ \end{matrix}

μ_{k} = \frac{∑_{i=1}^{N} ⟨π_{ik}⟩ x_{i}}{∑_{i=1}^{N} ⟨π_{ik}⟩}

w_{k} = \frac{∑_{i=1}^{N} ⟨π_{ik}⟩}{N}

9: end for 10: until model parameters converge

Algorithm 1: Estimate w_{k}, μ_{k}, Σ_{k} for 1 ≤ m ≤ K

log p(X; θ) = log ∫ p(X, Z; θ)dZ
= log ∫ Q(Z)\frac{p(X, Z; θ)}{Q(Z)}dZ  (11)
≥ ∫ Q(Z) log \frac{p(X, Z; θ)}{Q(Z)}dZ  (Jensen’s inequality)  (12)
= ∫ Q(Z) log p(X, Z; θ)dZ − ∫ Q(Z) log Q(Z)dZ  (13)

Hence we have the following lower-bound to log p(X; θ):

log p(X; θ) ≥ \underbrace{⟨log p(X, Z; θ)⟩_{Q(Z)}}_{exp. comp. log-lik.} + \underbrace{H [Q(Z)]}_{entropy of Q(Z)} = F(Q, θ)  (14)

Since Q(Z) is an arbitrary distribution, it is independent of θ. Hence in order to maximize the functional F(Q, θ) w.r.t. θ, it suffices to simply maximize ⟨log p(X, Z; θ)⟩_{Q(Z)}. (Hence the M-step).

Does this maximization achieve our aim? Eq. (14) shows that the functional F(Q, θ) is a lower bound to the quantity we are interested in. In which case, maximizing F(Q, θ) does not guarantee that we are improving log p(X; θ) at all! If however, we set Q(Z) = p(Z|X; θ) in Eq. (12), then we see that the lower bound in fact becomes an equality.

∫ Q(Z) log \frac{p(X, Z; θ)}{Q(Z)}dZ = ∫ p(Z|X; θ) log \frac{p(X, Z; θ)}{p(Z|X; θ)}dZ

= ∫ p(Z|X; θ) log \frac{p(Z|X; θ)p(X; θ)}{p(Z|X; θ)}dZ

= ∫ p(Z|X; θ) log p(X; θ)dZ

= log p(X; θ) ∫ p(Z|X; θ)dZ

= log p(X; θ)

This means that when computing the expected complete log-likelihood ⟨log p(X, Z; θ)⟩_{Q(Z)}, the expectation should be taken w.r.t. the true posterior p(Z|X; θ) of the hidden variables (the E step).

5 Examples

Figures 3, 4 and 5 show the result of fitting 2, 3, and 4 Gaussians respectively to a set of data points. In each figure the first plot shows the positions of the means and the relative covariances of each Gaussian, while the second shows the resulting estimated distribution obtained by marginalizing over the models as follows:

p(x; θ) = ∑_{k}^{K} p(x|z_{k} = 1; θ)p(z_{k} = 1)

Figure 6 demonstrates that one can model rather arbitrary non-gaussian distributions provided we have a sufficient number of mixture components. The tradeoff here is that too few components will fail to model the structure of the data, while too many will “overfit” the data. The model selection problem in this context is determining an appropriate compromise.

Figure 3: Fitting with 2 Gaussians 9 Figure 4: Fitting with 3 Gaussians

10 Figure 5: Fitting with 4 Gaussians

11 Figure 6: A sufficient number of mixture components can model arbitrary distributions

12