2
$\begingroup$

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.

Any help please?

$\endgroup$

1 Answer 1

1
$\begingroup$

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))

Contour plot

$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – user130885
    Oct 13, 2016 at 3:22
  • $\begingroup$ @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? $\endgroup$
    – BonStats
    Oct 13, 2016 at 3:29
  • $\begingroup$ 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. $\endgroup$
    – user130885
    Oct 13, 2016 at 3:41
  • $\begingroup$ 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? $\endgroup$
    – BonStats
    Oct 17, 2016 at 5:30
  • $\begingroup$ Ok I will send my code as soon as I can. Thanks a lot for your replay. $\endgroup$
    – user130885
    Oct 17, 2016 at 8:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy