Starting from:

$25

STAT542 - A3 - Solved

Implement the EM algorithm for a p-dimensional Gaussian mixture model with G components:

G

X

pk · N(x;µk,Σ). k=1

Store the parameters as a list in R with three components

•    prob: a G-dimensional probability vector;

•    mean: a p-by-G matrix with the k-th column being µk, the p-dimensional mean for the k-th Gaussian component;

•    Sigma: a p-by-p covariance matrix shared by all G components.

Your code should have the following structure.

Estep <- function(data, G, para){

# Your Code

# return the n-by-G probability matrix

}

Mstep <- function(data, G, para, post.prob){

# Your Code

# Return the updated parameters

}

myEM <- function(data, T, G, para){ for(t in 1:T){ post.prob <- Estep(data, G, para) para <- Mstep(data, G, para, post.prob)

} return(para)

}
You should test your code on the faithful data from the R package mclust with G = 2. The estimated parameters from your algorithm and the one from mclust after T = 10 iterations should be the same.

library(mclust)

n <- nrow(faithful) Z <- matrix(0, n, 2)

Z[sample(1:n, 120), 1] <- 1 Z[, 2] <- 1 - Z[, 1]

ini0 <- mstep(modelName="EEE", faithful, Z)$parameters

# Output from my EM alg

para0 <- list(prob = ini0$pro, mean = ini0$mean,

Sigma = ini0$variance$Sigma)

myEM(T=10, para=para0)

# Output from mclust

Rout <- em(modelName = "EEE", data = faithful, control = emControl(eps=0, tol=0, itmax = 10), parameters = ini0)$parameters

list(Rout$pro, Rout$mean, Rout$variance$Sigma)

More products