User:Timothee Flutre/Notebook/Postdoc/2011/12/14

From OpenWetWare

(Difference between revisions)
Jump to: navigation, search
m (Learn about mixture models and the EM algorithm)
(Learn about mixture models and the EM algorithm: add code for example)
Line 10: Line 10:
''(Caution, this is my own quick-and-dirty tutorial, see the references at the end for presentations by professional statisticians.)''
''(Caution, this is my own quick-and-dirty tutorial, see the references at the end for presentations by professional statisticians.)''
-
* '''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.
+
* '''Motivation''': a large part of any scientific activity is about measuring things, in other words collecting data, and it is not infrequent 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 <math>X = (x_1, x_2, ..., x_N)</math>. For the moment, we suppose that each observation <math>x_i</math> is univariate, ie. each corresponds to only one number.
* '''Data''': we have N observations, noted <math>X = (x_1, x_2, ..., x_N)</math>. For the moment, we suppose that each observation <math>x_i</math> 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 <math>K</math> clusters (see examples above). This means that we expect a subset of the observations to come from cluster <math>k=1</math>, another subset to come from cluster <math>k=2</math>, and so on.
+
* '''Hypothesis''': let's assume that the data are heterogeneous and that they can be partitioned into <math>K</math> clusters (see examples above). This means that we expect a subset of the observations to come from cluster <math>k=1</math>, another subset to come from cluster <math>k=2</math>, and so on.
* '''Model''': technically, we say that the observations were generated according to a [http://en.wikipedia.org/wiki/Probability_density_function density function] <math>f</math>. More precisely, this density is itself a mixture of densities, one per cluster. In our case, we will assume that each cluster <math>k</math> corresponds to a Normal distribution, which density is here noted <math>g</math>, with mean <math>\mu_k</math> and standard deviation <math>\sigma_k</math>. Moreover, as we don't know for sure from which cluster a given observation comes from, we define the mixture weight <math>w_k</math> to be the probability that any given observation comes from cluster <math>k</math>. As a result, we have the following list of parameters: <math>\theta=(w_1,...,w_K,\mu_1,...\mu_K,\sigma_1,...,\sigma_K)</math>. Finally, for a given observation <math>x_i</math>, we can write the model <math>f(x_i/\theta) = \sum_{k=1}^{K} w_k g(x_i/\mu_k,\sigma_k)</math> , with <math>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}</math>.
* '''Model''': technically, we say that the observations were generated according to a [http://en.wikipedia.org/wiki/Probability_density_function density function] <math>f</math>. More precisely, this density is itself a mixture of densities, one per cluster. In our case, we will assume that each cluster <math>k</math> corresponds to a Normal distribution, which density is here noted <math>g</math>, with mean <math>\mu_k</math> and standard deviation <math>\sigma_k</math>. Moreover, as we don't know for sure from which cluster a given observation comes from, we define the mixture weight <math>w_k</math> to be the probability that any given observation comes from cluster <math>k</math>. As a result, we have the following list of parameters: <math>\theta=(w_1,...,w_K,\mu_1,...\mu_K,\sigma_1,...,\sigma_K)</math>. Finally, for a given observation <math>x_i</math>, we can write the model <math>f(x_i/\theta) = \sum_{k=1}^{K} w_k g(x_i/\mu_k,\sigma_k)</math> , with <math>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}</math>.
Line 67: Line 67:
* '''EM algorithm''': ... <TO DO> ...
* '''EM algorithm''': ... <TO DO> ...
-
* '''Simulate data''':
+
* '''R code to simulate data''':
  #' Generate univariate observations from a mixture of Normals
  #' Generate univariate observations from a mixture of Normals
Line 91: Line 91:
  }
  }
-
* '''Implement the E step''':
+
* '''R code for the E step''':
  #' Return probas of latent variables given data and parameters from previous iteration
  #' Return probas of latent variables given data and parameters from previous iteration
Line 129: Line 129:
  }
  }
-
* '''Implement the M step''':
+
* '''R code for the M step''':
  #' Return ML estimates of parameters
  #' Return ML estimates of parameters
Line 193: Line 193:
   })
   })
  }
  }
 +
 +
* '''R code for the EM loop''':
 +
 +
... <TO DO> ...
 +
 +
* '''Example''': using the code above, I simulated data, ran the EM algorithm and plotted the results. It seems to work well, which was expected as the clusters are well separated from each other.
 +
 +
## simulate data
 +
K <- 3
 +
N <- 300
 +
simul <- GetUnivariateSimulatedData(K, N)
 +
data <- simul$obs
 +
 +
## run the EM algorithm
 +
params0 <- list(mus=runif(n=K, min=min(data), max=max(data)),
 +
                sigmas=rep(1, K),
 +
                mix.weights=rep(1/K, K))
 +
res <- EMalgo(data, params0, 10^(-3), 1000, 1)
 +
 +
## check its convergence
 +
plot(res$logliks, xlab="iterations", ylab="log-likelihood",
 +
      main="Convergence of the EM algorithm", type="b")
 +
 +
## plot the data along with the inferred densities
 +
png("mixture_univar_em.png")
 +
hist(data, breaks=30, freq=FALSE, col="grey", border="white", ylim=c(0,0.15),
 +
      main="Histogram of data overlaid with densities inferred by EM")
 +
rx <- seq(from=min(data), to=max(data), by=0.1)
 +
ds <- lapply(1:K, function(k){dnorm(x=rx, mean=res$params$mus[k], sd=res$params$sigmas[k])})
 +
f <- sapply(1:length(rx), function(i){
 +
  res$params$mix.weights[1] * ds[[1]][i] + res$params$mix.weights[2] * ds[[2]][i] + res$params$mix.weights[3] * ds[[3]][i]
 +
})
 +
lines(rx, f, col="red", lwd=2)
 +
dev.off()
* '''References''':
* '''References''':
-
** tutorial: document from Carlo Tomasi (Duke University)
+
** introduction (ch.1) of the PhD thesis from Matthew Stephens (Oxford, 2000)
-
** introduction to mixture models: PhD thesis from Matthew Stephens (Oxford, 2000)
+
** tutorial from Carlo Tomasi (Duke University)
-
** articles on the Bayesian approach: Diebolt and Robert (1994); Richardson and Green (1997); Jasra, Holmes and Stephens (2005)
+
<!-- ##### DO NOT edit below this line unless you know what you are doing. ##### -->
<!-- ##### DO NOT edit below this line unless you know what you are doing. ##### -->

Revision as of 16:50, 3 January 2012

Project name Main project page
Previous entry      Next entry

Learn about mixture models and the EM algorithm

(Caution, this is my own quick-and-dirty tutorial, see the references at the end for presentations by professional statisticians.)

  • Motivation: a large part of any scientific activity is about measuring things, in other words collecting data, and it is not infrequent 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 X = (x1,x2,...,xN). For the moment, we suppose that each observation xi is univariate, ie. each corresponds to only one number.
  • Hypothesis: let's assume that the data are heterogeneous and that they can be partitioned into K clusters (see examples above). This means that we expect a subset of the observations to come from cluster k = 1, another subset to come from cluster k = 2, and so on.
  • Model: technically, we say that the observations were generated according to a density function f. More precisely, this density is itself a mixture of densities, one per cluster. In our case, we will assume that each cluster k corresponds to a Normal distribution, which density is here noted g, with mean μk and standard deviation σk. Moreover, as we don't know for sure from which cluster a given observation comes from, we define the mixture weight wk to be the probability that any given observation comes from cluster k. As a result, we have the following list of parameters: θ = (w1,...,wK1,...μK1,...,σK). Finally, for a given observation xi, we can write the model f(x_i/\theta) = \sum_{k=1}^{K} w_k g(x_i/\mu_k,\sigma_k) , with 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: L(θ) = P(X / θ). 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: 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 l(θ). It gives the same solution because the log is monotonically increasing, but it's easier to derive the log-likelihood than the likelihood. Here is the whole formula:

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 Z1,...,Zi,...,ZN, one for each observation, such that Zi = k means that xi belongs to cluster k. Thanks to this, we can reinterpret the mixture weights: wk = P(Zi = k / θ). Moreover, we can now define the membership probabilities, one for each observation: 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 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. Introducing the latent variables corresponds to what is called the "missing data formulation" of the mixture problem.
  • MLE analytical formulae: a few important rules are required, but only from a high-school level in maths (see here). Let's start by finding the maximum-likelihood estimates of the mean of each cluster:

\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 μk, all the others means μl with l \ne k are constant, and thus disappear:

\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:

\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:

\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) = \sum_{i=1}^N \frac{1}{\sigma^2} p(k/i) (\mu_k - x_i)

By convention, we note \hat{\mu_k} the maximum-likelihood estimate of μk:

\frac{\partial l(\theta)}{\partial \mu_k}_{\mu_k=\hat{\mu_k}} = 0

Therefore, we finally obtain:

\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 derive the log-likelihood w.r.t. σk:

\frac{\partial l(\theta)}{\partial \sigma_k} = \sum_{i=1}^N p(k/i) (\frac{-1}{\sigma_k} + \frac{(x_i - \mu_k)^2}{\sigma_k^3})

And then we obtain the ML estimates for the standard deviation of each cluster:

\hat{\sigma_k} = \sqrt{\frac{\sum_{i=1}^N p(k/i) (x_i - \mu_k)^2}{\sum_{i=1}^N p(k/i)}}

The partial derivative of l(θ) w.r.t. wk is tricky. ... <TO DO> ...

\frac{\partial l(\theta)}{\partial w_k} = \sum_{i=1}^N (p(k/i) - w_k)

Finally, here are the ML estimates for the mixture weights:

\hat{w}_k = \frac{1}{N} \sum_{i=1}^N p(k/i)

  • EM algorithm: ... <TO DO> ...
  • R code to simulate data:
#' 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*(K-1), 6)
  sigmas <- runif(n=K, min=0.5, max=1.5)
  tmp <- floor(rnorm(n=K-1, 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.weights=ns/N))
}
  • R code for the E step:
#' 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.weights
Estep <- function(data, params){
  GetMembershipProbas(data, params$mus, params$sigmas, params$mix.weights)
}
#' 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.weights Kx1 vector of mixture weights w_k=P(zi=k/theta)
#' @return NxK matrix of membership probas
GetMembershipProbas <- function(data, mus, sigmas, mix.weights){
  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.weight){
      mix.weight * GetUnivariateNormalDensity(x, mu, sigma)}, mus, sigmas, mix.weights)))
    unlist(Map(function(mu, sigma, mix.weight){
      mix.weight * GetUnivariateNormalDensity(x, mu, sigma) / norm.const
    }, mus[-K], sigmas[-K], mix.weights[-K]))
  })), ncol=K-1, 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)*(x-mu)^2) )
}
  • R code for the M step:
#' Return ML estimates of parameters
#'
#' @param data Nx1 vector of observations
#' @param params list which components are mus, sigmas and mix.weights
#' @param membership.probas NxK matrix with entry i,k being P(zi=k/xi,theta)
Mstep <- function(data, params, membership.probas){
  params.new <- list()
  sum.membership.probas <- apply(membership.probas, 2, sum)
  params.new$mus <- GetMlEstimMeans(data, membership.probas,
                                    sum.membership.probas)
  params.new$sigmas <- GetMlEstimStdDevs(data, params.new$mus,
                                         membership.probas,
                                         sum.membership.probas)
  params.new$mix.weights <- GetMlEstimMixWeights(data, membership.probas,
                                                 sum.membership.probas)
  return(params.new)
}
#' Return ML estimates of the means (1 per cluster)
#'
#' @param data Nx1 vector of observations
#' @param membership.probas NxK matrix with entry i,k being P(zi=k/xi,theta)
#' @param sum.membership.probas Kx1 vector of sum per column of matrix above
#' @return Kx1 vector of means
GetMlEstimMeans <- function(data, membership.probas, sum.membership.probas){
  K <- ncol(membership.probas)
  sapply(1:K, function(k){
    sum(unlist(Map("*", membership.probas[,k], data))) /
      sum.membership.probas[k]
  })
}
#' Return ML estimates of the std deviations (1 per cluster)
#'
#' @param data Nx1 vector of observations
#' @param membership.probas NxK matrix with entry i,k being P(zi=k/xi,theta)
#' @param sum.membership.probas Kx1 vector of sum per column of matrix above
#' @return Kx1 vector of std deviations
GetMlEstimStdDevs <- function(data, means, membership.probas,
                              sum.membership.probas){
  K <- ncol(membership.probas)
  sapply(1:K, function(k){
    sqrt(sum(unlist(Map(function(p_ki, x_i){
      p_ki * (x_i - means[k])^2
    }, membership.probas[,k], data))) /
         sum.membership.probas[k])
  })
}
#' Return ML estimates of the mixture weights
#'
#' @param data Nx1 vector of observations
#' @param membership.probas NxK matrix with entry i,k being P(zi=k/xi,theta)
#' @param sum.membership.probas Kx1 vector of sum per column of matrix above
#' @return Kx1 vector of mixture weights
GetMlEstimMixWeights <- function(data, membership.probas,
                                 sum.membership.probas){
  K <- ncol(membership.probas)
  sapply(1:K, function(k){
    1/length(data) * sum.membership.probas[k]
  })
}
  • R code for the EM loop:

... <TO DO> ...

  • Example: using the code above, I simulated data, ran the EM algorithm and plotted the results. It seems to work well, which was expected as the clusters are well separated from each other.
## simulate data
K <- 3
N <- 300
simul <- GetUnivariateSimulatedData(K, N)
data <- simul$obs
## run the EM algorithm
params0 <- list(mus=runif(n=K, min=min(data), max=max(data)),
                sigmas=rep(1, K),
                mix.weights=rep(1/K, K))
res <- EMalgo(data, params0, 10^(-3), 1000, 1)
## check its convergence
plot(res$logliks, xlab="iterations", ylab="log-likelihood",
     main="Convergence of the EM algorithm", type="b")
## plot the data along with the inferred densities
png("mixture_univar_em.png")
hist(data, breaks=30, freq=FALSE, col="grey", border="white", ylim=c(0,0.15),
     main="Histogram of data overlaid with densities inferred by EM")
rx <- seq(from=min(data), to=max(data), by=0.1)
ds <- lapply(1:K, function(k){dnorm(x=rx, mean=res$params$mus[k], sd=res$params$sigmas[k])})
f <- sapply(1:length(rx), function(i){
  res$params$mix.weights[1] * ds1[i] + res$params$mix.weights[2] * ds2[i] + res$params$mix.weights[3] * ds3[i]
})
lines(rx, f, col="red", lwd=2)
dev.off()
  • References:
    • introduction (ch.1) of the PhD thesis from Matthew Stephens (Oxford, 2000)
    • tutorial from Carlo Tomasi (Duke University)


Personal tools