# How to plot the log likelihood associated with each iteration of EM algorithm and the density of each component?

I am new to R. I am working with mixture model and would like to get the plot the of the log likelihood associated with each iteration of the EM algorithm and the density of each component. I would like to do that manually.

I'm assuming you are able to calculate the log-likelihood and density in each iteration of your EM algorithm.

The general way to plot a function in R (of which a log-likelihood and density are) is to (a) create inputs for the function over the domain your interested in, (b) generate outputs from your function and save in a data frame, then (c) plot the function using the data frame. Here's an example of plotting a mixture model density using the 2D mixture model given by:

library(mvtnorm)
# example 2D mixture density
mixture_model_density <- function(x){
a <- 0.3
a * dmvnorm(x = x, mean = c(3, 2)) + (1 - a) * dmvnorm(x = x, mean = c(2.5,5.5))
}


(a) Create inputs:

grid_space <- 0.2
x <- seq(from = 0, to = 10, by = grid_space)
y <- seq(from = 0, to = 10, by = grid_space)
grid_mat <- expand.grid(x=x,y=y)


(b) Generate outputs:

 grid_mat\$mm_dens <- mixture_model_density(grid_mat)


(c) Plot results:

ggplot(grid_mat) + stat_contour(aes(x = x, y = y, z = mm_dens))


• Thanks a lot. Not sure how to find log likelihood for each iteration but I will work and see. I just replicate one code (because I am new to R) and tried to get the plot. Thanks a gain.
– user130885
Oct 13, 2016 at 3:22
• @Silver what's the package/code you're using for the EM algorithm? Maybe you can update your question with some more details on that? Oct 13, 2016 at 3:29
• I replicated this code "stats.stackexchange.com/questions/55132/…" (in the answer). But, I was confused why the author of the code used "mysum*(pi1(p1*dnorm) so here the author repeat pi1 two times. And I could not produce the plot.
– user130885
Oct 13, 2016 at 3:41
• In the initialisation line loglik[2]<-mysum(pi1*(log(pi1)+..., the first pi1 here is not tau1 because pi1 and pi2 each start at 0.5 therefore they sum to 1 and tau1 can be simplified. Look later in the code where loglik[k+1]<-mysum(tau1*(log(pi1)+.... As for the plot I need some more information. Maybe you can post a new question with a minimum working example with some example data and code of your first attempt? Oct 17, 2016 at 5:30
• Ok I will send my code as soon as I can. Thanks a lot for your replay.
– user130885
Oct 17, 2016 at 8:34