User:Timothee Flutre/Notebook/Postdoc/2011/12/14
From OpenWetWare
(→Learn about mixture models and the EM algorithm: add R code loglik) 
m (→Learn about mixture models and the EM algorithm) 

Line 10:  Line 10:  
''(Caution, this is my own quickanddirty tutorial, see the references at the end for presentations by professional statisticians.)''  ''(Caution, this is my own quickanddirty 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. It seems therefore natural to say that the samples come from a mixture of clusters. The aim is thus to recover from the data, ie. to infer, (i) how many clusters there are, (ii) what are the features of these clusters, and (  +  * '''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. It seems therefore natural to say that the samples come from a mixture of clusters. The aim is thus to recover from the data, ie. to infer, (i) how many clusters there are, (ii) what are the features of these clusters, and (iii) from which cluster each sample comes from. In the following, I will focus on points (ii) and (iii). 
* '''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. 
Revision as of 16:12, 27 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.)
The constraints are: and
We can now write the complete likelihood, ie. the likelihood of the augmented model, assuming all observations are independent: . And also the incomplete (or marginal) likelihood:
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 because of the constraints on the w_{k}. But we can handle it by writing them in terms of unconstrained variables γ_{k} (softmax function):
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 #' @param gap difference between all component means GetUnivariateSimulatedData < function(K=2, N=100, gap=6){ mus < seq(0, gap*(K1), gap) 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.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=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] }) }
GetLogLikelihood < function(data, mus, sigmas, mix.weights){ loglik < sum(sapply(data, function(x){ log(sum(unlist(Map(function(mu, sigma, mix.weight){ mix.weight * GetUnivariateNormalDensity(x, mu, sigma) }, mus, sigmas, mix.weights)))) })) return(loglik) }
EMalgo < function(data, params, threshold.convergence=10^(2), nb.iter=10, verbose=1){ logliks < vector() i < 1 if(verbose > 0) cat(paste("iter ", i, "\n", sep="")) membership.probas < Estep(data, params) params < Mstep(data, params, membership.probas) loglik < GetLogLikelihood(data, params$mus, params$sigmas, params$mix.weights) logliks < append(logliks, loglik) while(i < nb.iter){ i < i + 1 if(verbose > 0) cat(paste("iter ", i, "\n", sep="")) membership.probas < Estep(data, params) params < Mstep(data, params, membership.probas) loglik < GetLogLikelihood(data, params$mus, params$sigmas, params$mix.weights) if(loglik < logliks[length(logliks)]){ msg < paste("the loglikelihood is decreasing:", loglik, "<", logliks[length(logliks)]) stop(msg, call.=FALSE) } logliks < append(logliks, loglik) if(abs(logliks[i]  logliks[i1]) <= threshold.convergence) break } return(list(params=params, membership.probas=membership.probas, logliks=logliks, nb.iters=i)) }
## 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="loglikelihood", 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() It seems to work well, which was expected as the clusters are well separated from each other... The classification of each observation can be obtained via the following command: ## get the classification of the observations memberships < apply(res$membership.probas, 1, function(x){which(x > 0.5)}) table(memberships)
