Learn about mixture models and the EM algorithm
 Motivation and examples: a large part of any scientific activity is about measuring things, in other words collecting data, and it is not unfrequent to collect heterogeneous data. For instance, we measure the height of individuals without recording their gender, we measure the levels of expression of a gene in several individuals without recording which ones are healthy and which ones are sick, etc. It seems therefore natural to say that the samples come from a mixture of clusters. The aim is then to recover from the data, ie. to infer, (i) the values of the parameters of the probability distribution of each cluster, and (ii) from which cluster each sample comes from.
 Data: we have N observations, noted Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle X = (x_1, x_2, ..., x_N)}
. For the moment, we suppose that each observation Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle x_i}
is univariate, ie. each corresponds to only one number.
 Hypotheses and aim: let's assume that the data are heterogeneous and that they can be partitioned into Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle K}
clusters (see examples above). This means that we expect a subset of the observations to come from cluster Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle k=1}
, another subset to come from cluster Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle k=2}
, and so on.
 Model: technically, we say that the observations were generated according to a density function Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle f}
. More precisely, this density is itself a mixture of densities, one per cluster. In our case, we will assume that each cluster Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle k}
corresponds to a Normal distribution, here noted Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle g}
, of mean Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle \mu_k}
and standard deviation Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle \sigma_k}
. Moreover, as we don't know for sure from which cluster a given observation comes from, we define the mixture probability Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle w_k}
to be the probability that any given observation comes from cluster Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle k}
. As a result, we have the following list of parameters: Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle \theta=(w_1,...,w_K,\mu_1,...\mu_K,\sigma_1,...,\sigma_K)}
. Finally, for a given observation Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle x_i}
, we can write the model Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle f(x_i/\theta) = \sum_{k=1}^{K} w_k g(x_i/\mu_k,\sigma_k)}
, with Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle g(x_i/\mu_k,\sigma_k) = \frac{1}{\sqrt{2\pi} \sigma_k} \exp^{\frac{1}{2}(\frac{x_i  \mu_k}{\sigma_k})^2}}
.
 Likelihood: this corresponds to the probability of obtaining the data given the parameters: Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle L(\theta) = P(X/\theta)}
. We assume that the observations are independent, ie. they were generated independently, whether they are from the same cluster or not. Therefore we can write: Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle L(\theta) = \prod_{i=1}^N f(x_i/\theta)}
.
 Estimation: now we want to find the values of the parameters that maximize the likelihood. This reduces to (i) differentiating the likelihood with respect to each parameter, and then (ii) finding the value at which each partial derivative is zero. Instead of maximizing the likelihood, we maximize its logarithm, noted Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle l(\theta)}
. It gives the same solution because the log is monotonically increasing, but it's easier to derive the loglikelihood than the likelihood. Here is the whole formula:
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle l(\theta) = \sum_{i=1}^N log(f(x_i/\theta)) = \sum_{i=1}^N log( \sum_{k=1}^{K} w_k \frac{1}{\sqrt{2\pi} \sigma_k} \exp^{\frac{1}{2}(\frac{x_i  \mu_k}{\sigma_k})^2})}
 Latent variables: here it's worth noting that, although everything seems fine, a big information is lacking, we aim at finding the parameters defining the mixture but we don't know from which cluster each observation is coming... That's why we need to introduce the following N latent variables Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle Z_1,...,Z_i,...,Z_N}
, one for each observation, such that Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle Z_i=k}
means that Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle x_i}
belongs to cluster Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle k}
. Thanks to this, we can reinterpret the mixing probabilities: Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle w_k = P(Z_i=k/\theta)}
. Moreover, we can now define the membership probabilities, one for each observation: Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle P(Z_i=k/x_i,\theta) = \frac{w_k g(x_i/\mu_k,\sigma_k)}{\sum_{l=1}^K w_l g(x_i/\mu_l,\sigma_l)}}
. We will note these membership probabilities Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle p(k/i)}
as they will have a big role in the EM algorithm below. Indeed, we don't know the values taken by the latent variables, so we will have to infer their probabilities from the data.
 Technical details: a few important rules are required, but only from a highschool level in maths (see here). Let's start by finding the maximumlikelihood estimates of the mean of each cluster:
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle \frac{\partial l(\theta)}{\partial \mu_k} = \sum_{i=1}^N \frac{1}{f(x_i/\theta)} \frac{\partial f(x_i/\theta)}{\partial \mu_k}}
As we derive with respect to Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle \mu_k}
, all the others means Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle \mu_l}
with Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle l \ne k}
are constant, and thus disappear:
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle \frac{\partial f(x_i/\theta)}{\partial \mu_k} = w_k \frac{\partial g(x_i/\mu_k,\sigma_k)}{\partial \mu_k}}
And finally:
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle \frac{\partial g(x_i/\mu_k,\sigma_k)}{\partial \mu_k} = \frac{\mu_k  x_i}{\sigma_k^2} g(x_i/\mu_k,\sigma_k)}
Once we put all together, we end up with:
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle \frac{\partial l(\theta)}{\partial \mu_k} = \sum_{i=1}^N \frac{1}{\sigma^2} \frac{w_k g(x_i/\mu_k,\sigma_k)}{\sum_{l=1}^K w_l g(x_i/\mu_l,\sigma_l)} (\mu_k  x_i)}
By convention, we note Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle \hat{\mu_k}}
the maximumlikelihood estimate of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle \mu_k}
:
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle \frac{\partial l(\theta)}{\partial \mu_k}_{\mu_k=\hat{\mu_k}} = 0}
Therefore, we finally obtain:
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle \hat{\mu_k} = \frac{\sum_{i=1}^N p(k/i) x_i}{\sum_{i=1}^N p(k/i)}}
By doing the same kind of algebra, we also obtain the ML estimates for the standard deviation of each cluster:
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle \hat{\sigma_k} = \sqrt{\frac{\sum_{i=1}^N p(k/i) (x_i  \mu_k)^2}{\sum_{i=1}^N p(k/i)}}}
#' Generate univariate observations from a mixture of Normals
#'
#' @param K number of components
#' @param N number of observations
GetUnivariateSimulatedData < function(K=2, N=100){
mus < seq(0, 6*(K1), 6)
sigmas < runif(n=K, min=0.5, max=1.5)
tmp < floor(rnorm(n=K1, mean=floor(N/K), sd=5))
ns < c(tmp, N  sum(tmp))
clusters < as.factor(matrix(unlist(lapply(1:K, function(k){rep(k, ns[k])})),
ncol=1))
obs < matrix(unlist(lapply(1:K, function(k){
rnorm(n=ns[k], mean=mus[k], sd=sigmas[k])
})))
new.order < sample(1:N, N)
obs < obs[new.order]
rownames(obs) < NULL
clusters < clusters[new.order]
return(list(obs=obs, clusters=clusters, mus=mus, sigmas=sigmas,
mix.probas=ns/N))
}
#' Return probas of latent variables given data and parameters from previous iteration
#'
#' @param data Nx1 vector of observations
#' @param params list which components are mus, sigmas and mix.probas
Estep < function(data, params){
GetMembershipProbas(data, params$mus, params$sigmas, params$mix.probas)
}
#' Return the membership probabilities P(zi=k/xi,theta)
#'
#' @param data Nx1 vector of observations
#' @param mus Kx1 vector of means
#' @param sigmas Kx1 vector of std deviations
#' @param mix.probas Kx1 vector of mixing probas P(zi=k/theta)
#' @return NxK matrix of membership probas
GetMembershipProbas < function(data, mus, sigmas, mix.probas){
N < length(data)
K < length(mus)
tmp < matrix(unlist(lapply(1:N, function(i){
x < data[i]
norm.const < sum(unlist(Map(function(mu, sigma, mix.proba){
mix.proba * GetUnivariateNormalDensity(x, mu, sigma)}, mus, sigmas, mix.probas)))
unlist(Map(function(mu, sigma, mix.proba){
mix.proba * GetUnivariateNormalDensity(x, mu, sigma) / norm.const
}, mus[K], sigmas[K], mix.probas[K]))
})), ncol=K1, byrow=TRUE)
membership.probas < cbind(tmp, apply(tmp, 1, function(x){1  sum(x)}))
names(membership.probas) < NULL
return(membership.probas)
}
#' Univariate Normal density
GetUnivariateNormalDensity < function(x, mu, sigma){
return( 1/(sigma * sqrt(2*pi)) * exp(1/(2*sigma^2)*(xmu)^2) )
}
