rkHMM:

compare_to_depmix <- function(trn_mtx) {
  p_state_0 <- matrix(c(0.5, 0.5), 1, 2)

  s1 <- c(0.9, 0.9, 0.1, 0.9, 0.9)
  s2 <- c(0.2, 0.2, 0.8, 0.2, 0.2)
  p_obs_per_state <- array(c(s1, rep(1, length(s1)), s2, rep(1, length(s2))), c(5, 2, 2))
  joint_p_obs_per_state <-  apply(p_obs_per_state, c(1, 3), prod)

  # We transpose because depmixS4 goes from columns to rows
  trDens <- array(t(trn_mtx), dim = c(1, 2, 2))
  depmix_results <- depmixS4::fb(
    init = p_state_0,
    A = trDens,
    B = p_obs_per_state,
    useC = FALSE,
    homogeneous = TRUE,
    na.allow = FALSE,
    return.all = TRUE
  )

  f_p <- forward(p_state_0[1, ], trn_mtx, p_obs = joint_p_obs_per_state)
  f_b_p <- forward_backward(p_state_0[1, ], trn_mtx, p_obs = joint_p_obs_per_state)

  expect_true(all(abs(depmix_results$alpha - f_p$norm_forward_p) < 1e-8))
  expect_true(all(abs(f_b_p$alpha - f_p$norm_forward_p) < 1e-8))

  b_p <- backward(trn_mtx, p_obs = joint_p_obs_per_state, scale_factor = 1 / rowSums(f_p$forward_p))
  expect_true(all(abs(depmix_results$beta - b_p$norm_backward_p) < 1e-8))
  expect_true(all(abs(f_b_p$beta - b_p$norm_backward_p) < 1e-8))

  likelhood <- -sum(log(1 / rowSums(f_p$forward_p)))
  expect_true(all(abs(depmix_results$logLike - likelhood) < 1e-8))
  expect_true(all(abs(f_b_p$logLike - likelhood) < 1e-8))

  gamma_calc <- f_p$norm_forward_p * b_p$norm_backward_p / (1 / rowSums(f_p$forward_p))
  expect_true(all(abs(depmix_results$gamma - gamma_calc) < 1e-8))
  expect_true(all(abs(f_b_p$gamma - gamma_calc) < 1e-8))

  expect_true(all(abs(depmix_results$xi[-5 , , 2] - f_b_p$xi[-5 , , 2]) < 1e-8))
}

test_that("forward_backward", {
  trn_mtx <- matrix(c(0.7, 0.3, 0.3, 0.7), 2, 2, byrow = TRUE)
  compare_to_depmix(trn_mtx)

  trn_mtx <- matrix(c(0.5, 0.5, 0.3, 0.7), 2, 2, byrow = TRUE)
  compare_to_depmix(trn_mtx)
})