R/p_em.R

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

#' A mytoolbox function to compute the p function of an EM algorithm
#'
#' This function allows you to compute the p function for a mixture (EM)
#' @param theta parameters of the mixture that need estimation
#' @param y observations
#' @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.
#' @param p.0 absolute probability for latent variable.
#' @keywords mixture model
#' @export
#' @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.
#' @examples
#' p_em()
p_em = function(theta, y, p.0, mean_func, ...){

  ID = levels(y$ID)

  m = mean_func(theta, ...)

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

  p.old = matrix(0, nrow = length(ID), ncol = ncol(m))

  for(i in 1:nrow(p.old)){
    for(j in 1:ncol(p.old)){

      y.p = y[y$ID == ID[i], ]

      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)
      p.old[i,j] = sum(py.cx, na.rm=T) + log(pc.x)
    }

    mem.p = max(p.old[i,])
    p.temp = exp(p.old[i,] - mem.p)

    p.old[i,] = p.temp / sum(p.temp)
  }
  return(p.old)
}