R/computeH.R

In jiangrongo/HMM:

#' computeH
#'
#' \code{computeH} computes matrix, sumH, where sumH[k,j] = Sum over t (1->T-1) H_{k,j}(t)
#'
#'
#' H_{k,j}(t) = P(Z[t]=k , Z[t+1]=j | whole X = x).
#' We omitted the normalising constant 1/P(whole X = x)).
#'
#' @export
#' @param alpha a matrix, alpha, from Forward Algorithm
#' @param beta  a matrix, beta, from Backward Algorithm
#' @param X     a vector of observed states
#' @param trans a matrix of transition probability
#' @param u     a vector of means of Normal distribution for each state (emission probability)
#' @param sig   a vector of standard deviations of Normal distribution for each state (emission probability)
#' @return      A matrix, sumH, where sumH[k,j] = Sum over t (1->T-1) P(Z[t]=k , Z[t+1]=j | whole X = x)
#'
#' @examples
#' set.seed(1221)
#' df <- generateHMM(num=2,n=10)
#' alpha <- forwardAlg(df$X[,1], df$trans, df$u, df$sig)
#' beta <- backwardAlg(df$X[,1], df$trans, df$u, df$sig)
#' computeH(alpha, beta, df$X[,1], df$trans, df$u, df$sig)

computeH <- function(alpha, beta, X, trans, u, sig){
  # compute sumH directly
  M <- dim(trans)[1]
  n <- dim(alpha)[2]
  sumH <- matrix(0, M, M)
  for (t in 1:(n-1)){
    bf <- as.matrix(dnorm(X[t+1], u, sig) * beta[,t+1])
    sumH <- sumH + (as.matrix(alpha[, t]) %*% t(bf)) * trans
  }
  return(sumH)
}