fixes

loelschlaeger · Dec 13, 2023 · 9b94ae9 · 9b94ae9
1 parent cc43425
commit 9b94ae9
Show file tree

Hide file tree

Showing 18 changed files with 413 additions and 182 deletions.
diff --git a/NAMESPACE b/NAMESPACE
@@ -31,11 +31,17 @@ export(fHMM_parameters)
 export(fit_model)
 export(ll_hmm)
 export(npar)
+export(par2parCon)
 export(par2parUncon)
+export(parCon2par)
+export(parCon2parUncon)
+export(parUncon2par)
+export(parUncon2parCon)
 export(prepare_data)
 export(reorder_states)
 export(set_controls)
 export(simulate_hmm)
+export(viterbi)
 importFrom(Rcpp,evalCpp)
 importFrom(checkmate,assert_integerish)
 importFrom(checkmate,assert_number)

diff --git a/NEWS.md b/NEWS.md
@@ -6,6 +6,8 @@
 
 * Two more state-dependent distributions were added: `normal` and `poisson`.
 
+* The Viterbi algorithm can be directly accessed via `viterbi()`.
+
 * Renamed `simulate_data()` -> `simulate_hmm()` to make the functionality clearer. Furthermore, this function is now exported and can be used outside of the package to simulate HMM data.
 
 * `download_data()` no longer saves a .csv-file but returns the data as a `data.frame`. Its `verbose` argument is removed because the function no longer prints any messages.
@@ -19,46 +21,63 @@
 # fHMM 1.1.0
 
 * Extended the time horizon of saved data and updated models for demonstration.
+
 * The `download_data()` function now returns the data as a `data.frame` by default. However, specifying argument `file` still allows for saving the data as a .csv file.
+
 * The `plot.fHMM_model()` function now has the additional argument `ll_relative` (default is `TRUE`) to plot the relative log-likelihood values when `plot_type = "ll"`.
+
 * Significantly increased the test coverage and fixed minor bugs.
+
 * Changed color of time series plot from `"lightgray"` to `"black"` for better readability.
+
 * Added a title to the time series plot when calling `plot.fHMM_model(plot_type = "ts")`. Additionally, a time interval with arguments `from` and `to` can be selected to zoom into the data.
 
 # fHMM 1.0.3
 
 * Added the following methods for an `fHMM_model` object: `AIC()`, `BIC()`, `logLik()`, `nobs()`, `npar()`, `residuals()`.
+
 * The log-normal distribution can now be estimated by setting `sdds = "lnorm"` in the `controls` object.
 
 # fHMM 1.0.2
 
 * Fixed bug in `reorder_states()` that did not order the fine-scale parameter sets when the coarse-scale order was changed.
+
 * Fixed bug in `parameter_labels()` that returned the wrong order of parameter labels.
+
 * Changed plot type of simulated data to lines.
 
 # fHMM 1.0.1
 
 * In the vignette on controls, in the section about example specifications for `controls`, corrected `sdds = "gamma(mu = -1|1)"` to `sdds = "gamma(mu = 0.5|2)"` because mean of the Gamma distribution must be positive.
+
 * Added `digits` argument to `print.fHMM_predict()`.
+
 * Fixed bug in `reorder_states()` that allowed for misspecification of `state_order`.
+
 * Added option to `fit_model()` to initialize at the estimates of another model (#73).
 
 # fHMM 1.0.0
 
 * Enhanced the package by S3 classes.
+
 * Added more `controls` specifications.
+
 * Included a prediction function.
+
 * Improved documentations.
 
 # fHMM 0.3.0
 
 * Added vignettes.
+
 * Improved specification of `controls`.
+
 * Fixed minor bugs.
 
 # fHMM 0.2.0
 
 * Improved documentation of functions and README.
+
 * Improved specification of `controls`. (#37 and #38)
 
 # fHMM 0.1.0

diff --git a/R/decode_states.R b/R/decode_states.R
@@ -32,52 +32,6 @@ decode_states <- function(x, verbose = TRUE) {
     stop("'verbose' must be either TRUE or FALSE.", call. = FALSE)
   }
 
-  ### definition of the Viterbi algorithm for state decoding
-  viterbi <- function(observations, nstates, Gamma, mu, sigma, df, sdd) {
-    T <- length(observations)
-    delta <- oeli::stationary_distribution(Gamma, soft_fail = TRUE)
-    allprobs <- matrix(0, nstates, T)
-    for (n in seq_len(nstates)) {
-      if (sdd == "t") {
-        allprobs[n, ] <- (1 / sigma[n]) * stats::dt(
-          (observations - mu[n]) / sigma[n],
-          df[n]
-        )
-      }
-      if (sdd == "gamma") {
-        allprobs[n, ] <- stats::dgamma(observations,
-          shape = mu[n]^2 / sigma[n]^2,
-          scale = sigma[n]^2 / mu[n]
-        )
-      }
-      if (sdd == "normal") {
-        allprobs[n, ] <- stats::dnorm(observations, mean = mu[n], sd = sigma[n])                            
-      }
-      if (sdd == "lognormal") {
-        allprobs[n, ] <- stats::dlnorm(observations, meanlog = mu[n], 
-                                       sdlog = sigma[n])                            
-      }
-      if (sdd == "poisson") {
-        allprobs[n, ] <- stats::dpois(observations, lambda = mu[n])                            
-      }
-    }
-    xi <- matrix(0, nstates, T)
-    for (n in seq_len(nstates)) {
-      xi[n, 1] <- log(delta[n]) + log(allprobs[n, 1])
-    }
-    for (t in seq_len(T)[-1]) {
-      for (n in seq_len(nstates)) {
-        xi[n, t] <- max(xi[, t - 1] + log(Gamma[, n])) + log(allprobs[n, t])
-      }
-    }
-    iv <- numeric(T)
-    iv[T] <- which.max(xi[, T])
-    for (t in rev(seq_len(T - 1))) {
-      iv[t] <- which.max(xi[, t] + log(Gamma[, iv[t + 1]]))
-    }
-    return(iv)
-  }
-
   ### apply Viterbi algorithm
   par <- parUncon2par(x$estimate, x$data$controls)
   if (!x$data$controls$hierarchy) {
@@ -120,3 +74,102 @@ decode_states <- function(x, verbose = TRUE) {
   x$decoding <- decoding
   return(x)
 }
+
+#' @rdname decode_states
+#' @param observations
+#' A \code{numeric} \code{vector} of state-dependent observations.
+#' @param nstates
+#' The number of states.
+#' @param sdd
+#' A \code{character}, specifying the state-dependent distribution. One of 
+#' \itemize{
+#'   \item \code{"normal"} (the normal distribution),
+#'   \item \code{"lognormal"} (the log-normal distribution),
+#'   \item \code{"t"} (the t-distribution),
+#'   \item \code{"gamma"} (the gamma distribution),
+#'   \item \code{"poisson"} (the Poisson distribution).
+#' }
+#' @param Gamma
+#' A transition probability \code{matrix} of dimension \code{nstates}.
+#' @param mu
+#' A \code{numeric} vector of expected values for the state-dependent 
+#' distribution in the different states of length \code{nstates}.
+#' 
+#' For the gamma- or Poisson-distribution, \code{mu} must be positive.
+#' 
+#' @param sigma
+#' A positive \code{numeric} vector of standard deviations for the 
+#' state-dependent distribution in the different states of length \code{nstates}. 
+#' 
+#' Not relevant in case of a state-dependent Poisson distribution.
+#' 
+#' @param df
+#' A positive \code{numeric} vector of degrees of freedom for the 
+#' state-dependent distribution in the different states of length \code{nstates}. 
+#' 
+#' Only relevant in case of a state-dependent t-distribution.
+#' 
+#' @examples
+#' viterbi(
+#'   observations = c(1, 1, 1, 10, 10, 10),
+#'   nstates      = 2,
+#'   sdd          = "poisson",
+#'   Gamma        = matrix(0.5, 2, 2),
+#'   mu           = c(1, 10)
+#' )
+#' 
+#' @export
+
+viterbi <- function(
+    observations, nstates, sdd, Gamma, mu, sigma = NULL, df = NULL
+  ) {
+  T <- length(observations)
+  delta <- oeli::stationary_distribution(Gamma, soft_fail = TRUE)
+  allprobs <- matrix(0, nstates, T)
+  for (n in seq_len(nstates)) {
+    if (sdd == "t") {
+      allprobs[n, ] <- (1 / sigma[n]) * stats::dt(
+        (observations - mu[n]) / sigma[n],
+        df[n]
+      )
+    }
+    if (sdd == "gamma") {
+      allprobs[n, ] <- stats::dgamma(observations,
+                                     shape = mu[n]^2 / sigma[n]^2,
+                                     scale = sigma[n]^2 / mu[n]
+      )
+    }
+    if (sdd == "normal") {
+      allprobs[n, ] <- stats::dnorm(observations, mean = mu[n], sd = sigma[n])                            
+    }
+    if (sdd == "lognormal") {
+      allprobs[n, ] <- stats::dlnorm(observations, meanlog = mu[n], 
+                                     sdlog = sigma[n])                            
+    }
+    if (sdd == "poisson") {
+      allprobs[n, ] <- stats::dpois(observations, lambda = mu[n])                            
+    }
+  }
+  xi <- matrix(0, nstates, T)
+  for (n in seq_len(nstates)) {
+    xi[n, 1] <- log(delta[n]) + log(allprobs[n, 1])
+  }
+  for (t in seq_len(T)[-1]) {
+    for (n in seq_len(nstates)) {
+      xi[n, t] <- max(xi[, t - 1] + log(Gamma[, n])) + log(allprobs[n, t])
+    }
+  }
+  iv <- numeric(T)
+  iv[T] <- which.max(xi[, T])
+  for (t in rev(seq_len(T - 1))) {
+    iv[t] <- which.max(xi[, t] + log(Gamma[, iv[t + 1]]))
+  }
+  return(iv)
+}
+
+
+
+
+
+
+
diff --git a/R/fHMM_parameters.R b/R/fHMM_parameters.R
@@ -612,6 +612,7 @@ par2parUncon <- function(par, controls, use_parameter_labels = TRUE) {
 #' @rdname parameter_transformations
 #' @return
 #' For \code{parUncon2parCon}: a vector of constrained model parameters.
+#' @export
 
 parUncon2parCon <- function(parUncon, controls, use_parameter_labels = TRUE) {
   stopifnot(inherits(parUncon, "parUncon"))
@@ -722,6 +723,7 @@ parUncon2parCon <- function(parUncon, controls, use_parameter_labels = TRUE) {
 #' @rdname parameter_transformations
 #' @return
 #' For \code{parCon2par}: an object of class \code{\link{fHMM_parameters}}.
+#' @export
 
 parCon2par <- function(parCon, controls, use_parameter_labels = TRUE) {
   stopifnot(inherits(parCon, "parCon"))
@@ -816,6 +818,7 @@ parCon2par <- function(parCon, controls, use_parameter_labels = TRUE) {
 #' @rdname parameter_transformations
 #' @return
 #' For \code{par2parCon}: a vector of constrained model parameters.
+#' @export
 
 par2parCon <- function(par, controls, use_parameter_labels = TRUE) {
   stopifnot(inherits(par, "fHMM_parameters"))
@@ -830,6 +833,7 @@ par2parCon <- function(par, controls, use_parameter_labels = TRUE) {
 #' @rdname parameter_transformations
 #' @return
 #' For \code{parCon2parUncon}: a vector of unconstrained model parameters.
+#' @export
 
 parCon2parUncon <- function(parCon, controls, use_parameter_labels = TRUE) { 
   stopifnot(inherits(parCon, "parCon"))
@@ -844,6 +848,7 @@ parCon2parUncon <- function(parCon, controls, use_parameter_labels = TRUE) {
 #' @rdname parameter_transformations
 #' @return
 #' For \code{parUncon2par}: an object of class \code{fHMM_parameters}.
+#' @export
 
 parUncon2par <- function(
     parUncon, controls, use_parameter_labels = TRUE