4
$\begingroup$

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))
$\endgroup$
5
  • $\begingroup$ 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) $\endgroup$ Sep 21, 2016 at 15:15
  • $\begingroup$ 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? $\endgroup$
    – Alina
    Sep 21, 2016 at 17:27
  • $\begingroup$ 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. $\endgroup$ Sep 21, 2016 at 19:02
  • $\begingroup$ 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. $\endgroup$
    – Alina
    Sep 22, 2016 at 7:41
  • $\begingroup$ Since your question is closed here, can you perhaps change/remove your comment on the original question stats.stackexchange.com/a/131479/65037 $\endgroup$ Nov 15, 2016 at 8:56

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.