# Expectation-Maximization with a MLE package in R

As a follow up to one answer of the topic Expectation-Maximization with a coin toss: One of the user posted an R-code with MLE example almost a year ago (and his last online time here was 3 months ago, so I doubt he answers my question regarding his code).

In maximization step he multiplies likelihood of an A coin (row$weight_A) by a log-likelihood (llf_A). I do not understand why it is allowed to multiply a likelihood with a log-likelihood. I mean why wouldn't we multiple a likelihood with a likelihood or log-likelihood with a log-likelihood? And in the code he mixes it. require("stats4"); ## sample data from Do and Batzoglou ds<-data.frame(heads=c(5,9,8,4,7),n=c(10,10,10,10,10), coin=c("B","A","A","B","A"),weight_A=1:5*0) ## "baby likelihood" for a single observation llf <- function(heads, n, theta) { comb <- function(n, x) { #nCr function return(factorial(n) / (factorial(x) * factorial(n-x))) } if (theta<0 || theta >1) { # probabilities should be in [0,1] return(-Inf); } z<-comb(n,heads)* theta^heads * (1-theta)^(n-heads); return (log(z)) } ## the "E-M" likelihood function em <- function(theta_A,theta_B) { # expectation step: given current parameters, what is the likelihood # an observation is the result of tossing coin A (vs coin B)? ds$weight_A <<- by(ds, 1:nrow(ds), function(row) {
llf_A <- llf(row$heads,row$n, theta_A);
llf_B <- llf(row$heads,row$n, theta_B);

return(exp(llf_A)/(exp(llf_A)+exp(llf_B)));
})

# maximisation step: given params and weights, calculate likelihood of the sample
return(- sum(by(ds, 1:nrow(ds), function(row) {
llf_A <- llf(row$heads,row$n, theta_A);
llf_B <- llf(row$heads,row$n, theta_B);

return(row$weight_A*llf_A + (1-row$weight_A)*llf_B);
})))
}

est<-mle(em,start = list(theta_A=0.6,theta_B=0.5), nobs=NROW(ds))

• If we will exponentiate the last line, we will get likelihood(A) to the power row$weight_A multiply by likelihood(B) to the power 1 - row$weight_A. The final function returns actually log likelihood of the data. It was done in order not to multiply super small numbers (but I wonder why it was not done for the return of llf function itself, it is much safer to do z calculation on log scale) Sep 21, 2016 at 15:15
• row$weight_A is a likelihood of the data being produced by coin A, so it is exactly likelihood(A). llf_A is the log probability of tossing ones/zeros under a theta_A . If we exponentiate, we will have probability(theta_A) to the power likelihood(A) multiply by ... . But I still do not get why one would multiply likelihood(A) by a log probability (different scales of variables). ( logarithmizing avoids underflow but I asked why one mixes the log probabilities with non-log probabilities while multiplying) And why the whole term row$weight_A*llf_A is in log? Sep 21, 2016 at 17:27
• It looks like row$weight_A is a weight of the cluster and written as exp(llf_A)/(exp(llf_A)+exp(llf_B)). Which gives the likelihood exactly as it is defined as$\prod p^x (1-p)^{(1-x)}$. May be I am wrong, but it looks exactly like this. Sep 21, 2016 at 19:02 • Now I know where I had a mistake, I was always thinking about$\prod (p^{heads})(1-p)^{tails}$and therefore could not relate row$weight_A to the likelihood estimation. Sep 22, 2016 at 7:41
• Since your question is closed here, can you perhaps change/remove your comment on the original question stats.stackexchange.com/a/131479/65037 Nov 15, 2016 at 8:56