Nothing
rehmm_obj_fun <- function(paras, Y, N, K, quad_out) {
para_mat <- inflate_paras_rehmm(paras, N, K)
para_P1 <- para_mat$para_P1
para_a <- para_mat$para_a
para_b <- para_mat$para_b
para_alpha <- para_mat$para_alpha
para_beta <- para_mat$para_beta
-compute_llh_rehmm(Y, para_a, para_b, para_alpha, para_beta, para_P1, quad_out$nodes, quad_out$weights)
}
rehmm_gr_fun <- function(paras, Y, N, K, quad_out) {
para_mat <- inflate_paras_rehmm(paras, N, K)
para_P1 <- para_mat$para_P1
para_a <- para_mat$para_a
para_b <- para_mat$para_b
para_alpha <- para_mat$para_alpha
para_beta <- para_mat$para_beta
n <- length(Y)
para_a_gr <- matrix(0, K, K-1)
para_b_gr <- matrix(0, K, K-1)
para_alpha_gr <- matrix(0, K, N-1)
para_beta_gr <- matrix(0, K, N-1)
para_P1_gr <- rep(0, K-1)
for (index_seq in 1:n) {
seq_gr_res <- seq2llh_gr_rehmm(Y[[index_seq]], para_a, para_b, para_alpha, para_beta, para_P1, quad_out$nodes, quad_out$weights)
para_a_gr <- para_a_gr + seq_gr_res$para_a_gr
para_b_gr <- para_b_gr + seq_gr_res$para_b_gr
para_alpha_gr <- para_alpha_gr + seq_gr_res$para_alpha_gr
para_beta_gr <- para_beta_gr + seq_gr_res$para_beta_gr
para_P1_gr <- para_P1_gr + seq_gr_res$para_P1_gr
}
-deflate_paras_rehmm(para_a_gr, para_b_gr, para_alpha_gr, para_beta_gr, para_P1_gr)
}
#' Estimate latent traits in LHMM
#'
#' Compute MAP estimates of latent traits given LHMM parameters
#'
#' @param int_seqs a list of \code{n} action sequences where actions are coded as integers 0, ..., N-1
#' @param para_a \code{K} by \code{K-1} matrix. discrimination parameters of state transition probability matrix
#' @param para_b \code{K} by \code{K-1} matrix. location parameters of state transition probability matrix
#' @param para_alpha \code{K} by \code{N-1} matrix. discrimination parameters of state-action (emission) probability matrix
#' @param para_beta \code{K} by \code{N-1} matrix. location parameters of state-action (emission) probability matrix
#' @param para_P1 a vector of length \code{K-1}. parameters of initial state probability vector
#' @param n_pts number of quadrature points
#'
#' @return a vector of length \code{n}. Estimated latent traits.
#' @export
compute_theta <- function(int_seqs, para_a, para_b, para_alpha, para_beta, para_P1, n_pts) {
n_seq <- length(int_seqs)
N <- ncol(para_beta) + 1
K <- nrow(para_a)
out <- numeric(n_seq)
quad_out <- statmod::gauss.quad(n_pts, kind="hermite")
theta_list <- quad_out$nodes
weight_list <- quad_out$weights
P1 <- compute_P1_lhmm(para_P1)
# this loop is computationally intensive
for (index_seq in 1:n_seq) {
d_list <- rep(0, n_pts)
for (index_pt in 1:n_pts) {
paras <- compute_PQ_lhmm(sqrt(2)*theta_list[index_pt], para_a, para_b, para_alpha, para_beta)
d_list[index_pt] <- seq2llh(int_seqs[[index_seq]], P1, paras$P, paras$Q)
}
d_list <- d_list + log(weight_list)
max_d <- max(d_list)
theta_d_log <- max_d + log(sum(exp(d_list - max_d)))
out[index_seq] <- sum(exp(d_list - theta_d_log) * sqrt(2) * theta_list)
}
out
}
#' MMLE of LHMM
#'
#' Maximum marginalized likelihood estimation of LHMM.
#' Marginalization over latent trait is computed numerically using Guassian-Hermite quadratures from \code{\link{statmod}}.
#' Optimization is performed through \code{\link{optim}}.
#'
#'
#' @param action_seqs a list of \code{n} action sequences
#' @param K number of hidden states
#' @param paras a list of elements named \code{para_a}, \code{para_b}, \code{para_alpha}, \code{para_beta}, and \code{para_P1},
#' providing initial values of model parameters
#' @param n_pts number of quadrature points
#' @param verbose logical. If \code{TRUE}, progress messages are printed.
#' @param ... additional arguments passed to \code{\link{optim}}
#'
#' @return A list containing the following elements
#' \tabular{ll}{
#' \code{seqs} \tab action sequences coded in integers \cr
#' \tab \cr
#' \code{K} \tab number of hidden states \cr
#' \tab \cr
#' \code{N} \tab number of distinct actions \cr
#' \tab \cr
#' \code{paras_init} \tab a list containing initial values of parameters \cr
#' \tab \cr
#' \code{paras_est} \tab a list containing parameter estimates \cr
#' \tab \cr
#' \code{theta_est} \tab a vector of length \code{n}. estimated latent traits \cr
#' \tab \cr
#' \code{init_mllh} \tab initial value of the marginalized likelihood function \cr
#' \tab \cr
#' \code{opt_mllh} \tab maximized marginalized likelihood function \cr
#' \tab \cr
#' \code{opt_res} \tab object returned by \code{\link{optim}} \cr
#' }
#'
#' @examples
#' # generate data
#' paras_true <- sim_lhmm_paras(5, 2)
#' sim_data <- sim_lhmm(10, paras_true, 3, 5)
#' # randomly initialize parameters
#' paras_init <- sim_lhmm_paras(5, 2)
#' # fit model
#' lhmm_res <- lhmm(sim_data$seqs, 2, paras_init)
#'
#' @export
#'
lhmm <- function(action_seqs, K, paras, n_pts = 100, verbose = TRUE, ...) {
n <- length(action_seqs) # sample size
# create action set
action_set <- unique(unlist(action_seqs))
N <- length(action_set)
action2id <- 1:N
names(action2id) <- action_set
# convert action sequences to integer sequences consisting of 0, 1, ..., N-1
int_seqs <- list(n)
for (i in 1:n) {
seqt <- action_seqs[[i]]
int_seqs[[i]] <- action2id[seqt] - 1
}
quad_out <- statmod::gauss.quad(n_pts, kind="hermite")
# flatten parameters for optimization
paras_init_rehmm <- deflate_paras_rehmm(paras$para_a, paras$para_b, paras$para_alpha, paras$para_beta, paras$para_P1)
# objective function value at intial parameter values
obj_val_init_rehmm <- rehmm_obj_fun(paras_init_rehmm, int_seqs, N, K, quad_out)
if (verbose) cat("Optimizing obj fun...\n")
# minimize the marginalized maximum likelihood function
opt_res_rehmm <- optim(paras_init_rehmm, fn = rehmm_obj_fun, gr = rehmm_gr_fun, Y = int_seqs, N = N, K = K, quad_out = quad_out, method="BFGS", ...)
# objective function value at optimized parameter values
obj_val_opt_rehmm <- opt_res_rehmm$value
# optimized parameter values
opt_paras_rehmm <- inflate_paras_rehmm(opt_res_rehmm$par, N, K)
rownames(opt_paras_rehmm$para_a) <- paste0("state", 1:K)
rownames(opt_paras_rehmm$para_b) <- paste0("state", 1:K)
rownames(opt_paras_rehmm$para_alpha) <- paste0("state", 1:K)
rownames(opt_paras_rehmm$para_beta) <- paste0("state", 1:K)
names(opt_paras_rehmm$para_P1) <- paste0("state", 2:K)
colnames(opt_paras_rehmm$para_a) <- paste0("state", 2:K)
colnames(opt_paras_rehmm$para_b) <- paste0("state", 2:K)
colnames(opt_paras_rehmm$para_alpha) <- action_set[-1]
colnames(opt_paras_rehmm$para_beta) <- action_set[-1]
# compute MAE of theta
if (verbose) cat("Computing theta...\n")
theta_est <- compute_theta(int_seqs, opt_paras_rehmm$para_a, opt_paras_rehmm$para_b, opt_paras_rehmm$para_alpha, opt_paras_rehmm$para_beta, opt_paras_rehmm$para_P1, n_pts)
out <- list(seqs = int_seqs, K = K, N = N,
paras_init = inflate_paras_rehmm(paras_init_rehmm, N, K),
paras_est = opt_paras_rehmm,
theta_est = theta_est,
init_mllh = -obj_val_init_rehmm,
opt_mllh = -obj_val_opt_rehmm,
opt_res = opt_res_rehmm)
return(out)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.