User:Timothee Flutre/Notebook/Postdoc/2011/12/14
From OpenWetWare
(Difference between revisions)
(→Learn about mixture models and the EM algorithm: add code M step) 
(→Learn about mixture models and the EM algorithm: replace mixing proportions by mixture weights) 

Line 16:  Line 16:  
* '''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.  * '''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.  
  * '''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, here noted <math>g</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>. 
* '''Likelihood''': this corresponds to the probability of obtaining the data given the parameters: <math>L(\theta) = P(X/\theta)</math>. 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: <math>L(\theta) = \prod_{i=1}^N f(x_i/\theta)</math>.  * '''Likelihood''': this corresponds to the probability of obtaining the data given the parameters: <math>L(\theta) = P(X/\theta)</math>. 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: <math>L(\theta) = \prod_{i=1}^N f(x_i/\theta)</math>.  
Line 23:  Line 23:  
<math>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})</math>  <math>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})</math>  
  * '''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 <math>Z_1,...,Z_i,...,Z_N</math>, one for each observation, such that <math>Z_i=k</math> means that <math>x_i</math> belongs to cluster <math>k</math>. Thanks to this, we can reinterpret the  +  * '''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 <math>Z_1,...,Z_i,...,Z_N</math>, one for each observation, such that <math>Z_i=k</math> means that <math>x_i</math> belongs to cluster <math>k</math>. Thanks to this, we can reinterpret the mixture weights: <math>w_k = P(Z_i=k/\theta)</math>. Moreover, we can now define the membership probabilities, one for each observation: <math>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)}</math>. We will note these membership probabilities <math>p(k/i)</math> 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 highschool level in maths (see [http://en.wikipedia.org/wiki/Differentiation_%28mathematics%29#Rules_for_finding_the_derivative here]). Let's start by finding the maximumlikelihood estimates of the mean of each cluster: 
<math>\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}</math>  <math>\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}</math>  
Line 61:  Line 61:  
<math>\frac{\partial l(\theta)}{\partial w_k} = \sum_{i=1}^N (p(k/i)  w_k)</math>  <math>\frac{\partial l(\theta)}{\partial w_k} = \sum_{i=1}^N (p(k/i)  w_k)</math>  
  Finally, here are the ML estimates for the mixture  +  Finally, here are the ML estimates for the mixture weights: 
<math>\hat{w}_k = \frac{1}{N} \sum_{i=1}^N p(k/i)</math>  <math>\hat{w}_k = \frac{1}{N} \sum_{i=1}^N p(k/i)</math>  
Line 88:  Line 88:  
clusters < clusters[new.order]  clusters < clusters[new.order]  
return(list(obs=obs, clusters=clusters, mus=mus, sigmas=sigmas,  return(list(obs=obs, clusters=clusters, mus=mus, sigmas=sigmas,  
  mix.  +  mix.weights=ns/N)) 
}  }  
Line 96:  Line 96:  
#'  #'  
#' @param data Nx1 vector of observations  #' @param data Nx1 vector of observations  
  #' @param params list which components are mus, sigmas and mix.  +  #' @param params list which components are mus, sigmas and mix.weights 
Estep < function(data, params){  Estep < function(data, params){  
  GetMembershipProbas(data, params$mus, params$sigmas, params$mix.  +  GetMembershipProbas(data, params$mus, params$sigmas, params$mix.weights) 
}  }  
Line 106:  Line 106:  
#' @param mus Kx1 vector of means  #' @param mus Kx1 vector of means  
#' @param sigmas Kx1 vector of std deviations  #' @param sigmas Kx1 vector of std deviations  
  #' @param mix.  +  #' @param mix.weights Kx1 vector of mixture weights w_k=P(zi=k/theta) 
#' @return NxK matrix of membership probas  #' @return NxK matrix of membership probas  
  GetMembershipProbas < function(data, mus, sigmas, mix.  +  GetMembershipProbas < function(data, mus, sigmas, mix.weights){ 
N < length(data)  N < length(data)  
K < length(mus)  K < length(mus)  
tmp < matrix(unlist(lapply(1:N, function(i){  tmp < matrix(unlist(lapply(1:N, function(i){  
x < data[i]  x < data[i]  
  norm.const < sum(unlist(Map(function(mu, sigma, mix.  +  norm.const < sum(unlist(Map(function(mu, sigma, mix.weight){ 
  mix.proba * GetUnivariateNormalDensity(x, mu, sigma)}, mus, sigmas, mix.  +  mix.proba * GetUnivariateNormalDensity(x, mu, sigma)}, mus, sigmas, mix.weights))) 
  unlist(Map(function(mu, sigma, mix.  +  unlist(Map(function(mu, sigma, mix.weight){ 
mix.proba * GetUnivariateNormalDensity(x, mu, sigma) / norm.const  mix.proba * GetUnivariateNormalDensity(x, mu, sigma) / norm.const  
  }, mus[K], sigmas[K], mix.  +  }, mus[K], sigmas[K], mix.weights[K])) 
})), ncol=K1, byrow=TRUE)  })), ncol=K1, byrow=TRUE)  
membership.probas < cbind(tmp, apply(tmp, 1, function(x){1  sum(x)}))  membership.probas < cbind(tmp, apply(tmp, 1, function(x){1  sum(x)}))  
Line 134:  Line 134:  
#'  #'  
#' @param data Nx1 vector of observations  #' @param data Nx1 vector of observations  
  #' @param params list which components are mus, sigmas and mix.  +  #' @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)  #' @param membership.probas NxK matrix with entry i,k being P(zi=k/xi,theta)  
Mstep < function(data, params, membership.probas){  Mstep < function(data, params, membership.probas){  
Line 144:  Line 144:  
membership.probas,  membership.probas,  
sum.membership.probas)  sum.membership.probas)  
  params.new$mix.  +  params.new$mix.weights < GetMlEstimMixWeights(data, membership.probas, 
  +  sum.membership.probas)  
return(params.new)  return(params.new)  
}  }  
Line 180:  Line 180:  
}  }  
  #' Return ML estimates of the  +  #' Return ML estimates of the mixture weights 
#'  #'  
#' @param data Nx1 vector of observations  #' @param data Nx1 vector of observations  
#' @param membership.probas NxK matrix with entry i,k being P(zi=k/xi,theta)  #' @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  #' @param sum.membership.probas Kx1 vector of sum per column of matrix above  
  #' @return Kx1 vector of  +  #' @return Kx1 vector of mixture weights 
  +  GetMlEstimMixWeights < function(data, membership.probas,  
  +  sum.membership.probas){  
K < ncol(membership.probas)  K < ncol(membership.probas)  
sapply(1:K, function(k){  sapply(1:K, function(k){ 
Revision as of 12: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 quickanddirty tutorial, see the references at the end for presentations by professional statisticians.)
As we derive with respect to μ_{k}, all the others means μ_{l} with are constant, and thus disappear:
And finally:
Once we put all together, we end up with:
By convention, we note the maximumlikelihood estimate of μ_{k}:
Therefore, we finally obtain:
By doing the same kind of algebra, we derive the loglikelihood w.r.t. σ_{k}:
And then we obtain the ML estimates for the standard deviation of each cluster:
The partial derivative of l(θ) w.r.t. w_{k} is tricky. ... <TO DO> ...
Finally, here are the ML estimates for the mixture weights:
#' 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.weights=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.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.proba * GetUnivariateNormalDensity(x, mu, sigma)}, mus, sigmas, mix.weights))) unlist(Map(function(mu, sigma, mix.weight){ mix.proba * GetUnivariateNormalDensity(x, mu, sigma) / norm.const }, mus[K], sigmas[K], mix.weights[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) ) }
#' 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] }) }
