R/validHMM.R

In jiangrongo/HMM:

#' validMeans
#'
#' \code{validMeans} performs simulations and inference on means of HMM.
#'
#'
#' @export
#' @param nrep  the size of sample generated
#' @param trans a matrix of transition probability
#' @param M1    the number of hidden states
#' @return A violin plot of estimated means
#'
#' @examples
#' set.seed(1221)
#' validMeans()
#'
validMeans <- function(nrep=50,
                     trans=matrix(c(0, 3, 7, 5, 2, 3, 5, 5, 0),
                                  3, 3, byrow = TRUE)/10,
                     M1=3){
  # M1=3  don't change. If change, also change section of violin plot
  theta <- vector("list", length=nrep)
  constr <- trans>0  # >0 then 1, otherwise 0

  for (i in 1:nrep){
    # length of Markov chain set to multiple of M1^2+2*M1, the number of parameters
    df <- generateHMM(num=10, n=10*(M1^2+2*M1), trans=trans)$X
    theta[[i]] <- estimateHMM(df, M=M1, constr=constr, tol = 0.001)
  }

  # initialised empty matrix for estimated means
  Est <- matrix(0, M1*nrep, 3)
  for (i in 1:nrep){
    pos <- ((i-1)*M1+1): (i*M1)
    Est[pos,1] <- theta[[i]]$u[,  theta[[i]]$iter+1]
    Est[pos,2] <- theta[[i]]$sig[,theta[[i]]$iter+1]
    # access 3 elements in trans: 1->2, 2->3, 3->1
    Est[pos,3] <- theta[[i]]$trans[[theta[[i]]$iter+1]][c(4,8,3)]
  }

  # create data frame
  df_Est <- data.frame(Est)
  colnames(df_Est) <- c("u", "sig", "trans")
  df_Est$state <- factor(rep(1:M1, nrep))

  # violin plot
  # to make the plot on the same level and comparable to its own true mean
  # we subtract the estimate by their true mean (Note: we set mean as state label)
  df_Est$modi_u <- df_Est$u - c(-5,0,5)
  ggplot2::ggplot(df_Est, ggplot2::aes(x=state, y=modi_u, fill=state)) +
    ggplot2::geom_violin() +
    ggplot2::geom_abline(slope=0, intercept=0, lwd=.8, color=2)
}