R/Q_em.R

In ick003/SpTMixture: Mixture models with spatio-temporal dependencies

#' A mytoolbox function to compute the Q function of an EM algorithm
#'
#' This function allows you to compute the Q function for a mixture (EM)
#' @param theta parameters of the mixture that need estimation
#' @param y observations
#' @param p.theta result from the expectation step (EM - probabilities)
#' @param mean_func function to calculate the mean for each mixture. The result is a matrix where the number of columns is the number of element in the mixture.
#' @keywords mixture model
#' @references Dempster, A. and Laird, N. and Rubin, D. (1977) Maximum likelihood from incomplete data via the EM algorithm. \emph{JRSS B} \bold{39} 1--38.
#' @export
#' @examples
#' Q_em()
Q_em = function(theta,y, p.theta,p.0, mean_func, ...){

  ID = levels(y$ID)

  m = mean_func(theta, ...)

  sigma = theta[(length(theta) - ncol(p.theta) + 1):length(theta)]

  q = 0
#  q.hist = matrix(NA, nrow = nrow(p.theta), ncol = ncol(p.theta))
  for(i in 1:nrow(p.theta)){
    y.p = y[y$ID == ID[i], ]
    for(j in 1:ncol(p.theta)){

      m.p = m[, j]

      py.cx = dnorm(y.p$obs, m.p, exp(sigma[j]), log = T)

      pc.x = max(p.0[i,j], 10^-290)

      lpc.x.p.theta = log(pc.x)* p.theta[i,j]

      if(is.nan(lpc.x.p.theta)){lpc.x.p.theta = 0}

      q = q + (sum(py.cx, na.rm=T)* p.theta[i,j] + lpc.x.p.theta)

      #q.hist[i,j] = (sum(py.cx, na.rm=T) + lpc.x) * p.theta[i,j]
    }

  }
  if(is.na(q)){q = -Inf}
  return(q)
}