**--Jean Arreola--**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

## Mixture model

A Gaussian mixture model is a probabilistic way of representing

subpopulations within an overall population. We only observe the data,

not the subpopulation from which observation belongs.

We have $N$ random variables observed, each distributed according to a

mixture of K gaussian components. Each gaussian has its own

parameters, and we should be able to estimate the category using

Expectation Maximization, as we are using a latent variables model.

Now, in a bayesian scenario, each parameter of each gaussian is also a

random variable, as well as the mixture weights. To estimate the

distributions we use Variational Inference, which can be seen as a

generalization of the EM algorithm. Be sure to check this

book to learn all the theory behind gaussian

mixtures and variational inference.

Here is my implementation for Variational Gaussian Mixture Model.

```
#Variational Gaussian Mixture Model
#Constant for Dirichlet Distribution
dirConstant <- function(alpha){
res <- 1
for(i in 1:length(alpha)){
res <- res * gamma(alpha[i])
}
return(gamma(sum(alpha))/res)
}
BWishart <- function(W, v){
D <- ncol(W)
elem1 <- (det(W))^(-v/2)
elem2 <- (2^(v*D/2)) * (pi^(D*(D-1)/4))
elem3 <- 1
for(i in 1:D){
elem3 <- elem3 * gamma((v+1-i)/2)
}
return(elem1 / (elem2 * elem3))
}
#Log precision expected value
espLnPres <- function(W, v){
res <- 0
D <- ncol(W)
for(i in 1:D){
res <- res + digamma((v+1-i)/2)
}
res <- res + D*log(2) + log(det(W))
return(res)
}
#Wishart distribution entropy
entropyWishart <- function(W, v){
D <- ncol(W)
return(-log(BWishart(W,v)) - ((v-D-1)/2) * espLnPres(W,v) + (v*D)/2)
}
# Estimating mixture parameters
vgmm <- function(X, K, iter = 100, eps = 0.001){
D <- ncol(X)
N <- nrow(X)
#Hyperparameters initialization
m0 <- rep(0, D) # mean
W0 <- diag(D) # precision
v0 <- D # degrees of freedom: n > p-1
alpha0 <- 1/K # Dirichlet parameter
beta0 <- 1 # Variance for mean
#For each category
#Initialize the means with centroids from k-means
mk <- kmeans(X,K)$centers
#Initialize presicions with diagonal matrix
Wk <- array(0, c(D, D, K))
for(i in 1:K)
Wk[,,i] <- W0
vk <- rep(v0, K)
#Initialize hyperparameters
betak <- rep(beta0, K)
alphak <- rep(alpha0,K)
# Necessary terms for calculate responsabilities
ln_pres <- rep(0,K)
ln_pi <- rep(0,K)
E_mu_pres <- matrix(0, N, K)
# Iterate
for(it in 1:iter){
#Responsabilities
r <- matrix(0,N, K)
##################### Variational E-Step ##########################33
for(i in 1:K){
#Log precision
ln_pres[i] <- 0
for(j in 1:D){
ln_pres[i] <- ln_pres[i] + digamma((vk[i] + 1 - j) /2)
}
ln_pres[i] <- ln_pres[i] + D * log(2) + log(det(Wk[,,i]))
alpha <- sum(alphak)
ln_pi[i] <- digamma(alphak[i]) - digamma(alpha)
#E[mu,pres] (expected value of joint distribution of mu and pres)
for(k in 1:N){
E_mu_pres[k,i] <- (D / betak[i]) + vk[i] * t(X[k,] - mk[i,]) %*%
Wk[,,i] %*% (X[k,] - mk[i,]) #10.64
r[k,i] <- ln_pi[i] + 0.5 * ln_pres[i] - (D/2) *log(2*pi) -
0.5 * E_mu_pres[k,i]
}
}
# Exp-log-sum trick for numerical stability
rho <- apply(r, 1, function(x){
offset <- max(x)
y <- x - offset
return(exp(y)/sum(exp(y)))
})
rho <- t(rho)
########################### Variational M-Step ##################################
# Auxiliary statistics
Nk <- apply(rho, 2, sum)
# Update means
xBark <- matrix(0, K, D)
for(i in 1:K){
xBark[i,] <- colSums(rho[,i] * X) / Nk[i]
}
# Update covariances
Sk <- array(0, c(D,D,K))
for(i in 1:K){
sum_sk <- 0
for(j in 1:N){
sum_sk <- sum_sk + rho[j,i] * (X[j,] - xBark[i,]) %*% t(X[j,] - xBark[i,])
}
Sk[,,i] <- sum_sk / Nk[i]
}
# Update hyperparameters
for(i in 1:K){
betak[i] <- beta0 + Nk[i]
mk[i,] <- (1/betak[i]) * (beta0 * m0 + Nk[i] * xBark[i,])
Wk[,,i] <- solve(solve(W0) + Nk[i] * Sk[,,i] +
((beta0 * Nk[i]) / (beta0 + Nk[i])) *
(xBark[i,] - m0) %*% t(xBark[i,] - m0))
vk[i] <- v0 + Nk[i]
}
#ELBO (Evidence Lower Bound)
# ELBO is a sum of seven terms
term1 <- 0 #10.71
for(i in 1:K){
term1 <- term1 + Nk[i] * (ln_pres[i] - (D / betak[i]) -
vk[i] * sum(diag(Sk[,,i] %*% Wk[,,i])) -
vk[i] * ( t(xBark[i,] - mk[i,]) %*% Wk[,,i] %*%
(xBark[i,] - mk[i,])) -
D * log(2 * pi))
}
term1 <- 0.5 * term1
term2 <- 0 #10.72
for(i in 1:N){
for(j in 1:K){
term2 <- term2 + (rho[i,j] * ln_pi[j])
}
}
term3 <- 0 #10.73
for(i in 1:K){
term3 <- term3 + ln_pi[i]
}
term3 <- term3 * (alpha0 -1) + log(dirConstant(alpha0))
term4 <- 0 #10.74
sub <- 0
for(i in 1:K){
term4 <- term4 + D * log(beta0 / (2 * pi)) + ln_pres[i]-
((D * beta0)/betak[i]) - beta0 * vk[i] *
t(mk[i,]-m0) %*% Wk[,,i] %*% (mk[i,]-m0)
}
term4 <- 0.5 * term4 + K * log(BWishart(W0,v0))
for(i in 1:K){
sub <- sub + vk[i] * sum(diag(solve(W0) %*% Wk[,,i]))
}
term4 <- term4 + sum(ln_pres) * ((v0-D-1)/2) - 0.5 * sub
term5 <- 0 #10.75
for(i in 1:N){
for(j in 1:K){
stand <- rho[i,j] * log(rho[i,j])
if(!is.finite(stand))
stand <- 0
term5 <- term5 + stand
}
}
term6 <- 0 #10.76
for(i in 1:K){
term6 <- term6 + (alphak[i]-1) * ln_pi[i]
}
term6 <- term6 + log(dirConstant(alphak))
term7 <- 0 #10.77
for(i in 1:K){
term7 <- term7 + 0.5 * ln_pres[i] + (D/2) * log(betak[i]/(2 * pi)) -
(D/2) - entropyWishart(Wk[,,i], vk[i])
}
if(it > 1){
prevELBO <- ELBO
}
ELBO <- term1 + term2 + term3 + term4 - term5 - term6 - term7
# Convergence criteria
if(it > 1 && is.finite(ELBO)){
if(abs(ELBO - prevELBO) < eps){
break
}
}
}
# Return responsabilities, ELBO, covariances and means
# (You can add whatever parameters (or hyperparameters) you need)
lista <- list("rho" = rho, "ELBO" = ELBO, "Wk" = Wk, "mk" = mk)
return(lista)
}
```

## Applications

Gaussian Mixture Models can be seen as a form of clustering, but each

observation will belong to all clusters simultaneously, as we are

estimating the probabilities for belonging to each gaussian

distribution. This is called “soft clustering”, as opposed to other

algorithms like k-means, which is a “hard clustering technique” (each

observation belongs to only one cluster). As a matter of fact, k-means

is a special case of a gaussian mixture when the variances all are the

same, and there aren’t covariances (so all the clusters will have a

circular shape).

A consequence of this is that gaussian mixture are more flexible than

k-means because the clusters can have an “elliptical form”. In

particular, in image segmentation, gaussian mixture are the prefered

algorithm. For example, in image matting (segment an image by

background and foreground pixels), GMM are a natural choice because

each pixel will have a probability for belongin to the foreground and

the background.

## Eigenfaces

In this post, we will use variational GMM to do face detection. We

will use the faces94

dataset, and

choose the most probable category for each face.

The representation that I choose for the images are the Eigenfaces,

which are the eigenvectors of the matrix of faces (each column is an

image and each row has all the pixels values of the image). It’s

important to note that the images have to be centered (sustract the mean).

To reduce dimensionality, we will work with the eigenvectors of the

matrix X’X, so we will have instead a matrix of N x N.

## Results

The first five eigenfaces:

Now the results of the classification:

We can see that the algorithm only misclassified one point. Notice

that the groups are almost linearly separable, so eigenfaces was an

extremely helpful representation.

# Final thoughts

A gaussian mixture model is a powerful technique for unsupervised

learning. With Variational Inference, we can give more abilities to

the mixture, like working with missing values, or adding additional

levels to the hierarchical model. GMM are also the principles for

learning advances models like Hidden Markov Models.

**leave a comment**for the author, please follow the link and comment on their blog:

**--Jean Arreola--**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.