Variance-Covariance Matrix for Hidden Markov Models

Description

Computes an approximate variance-covariance matrix of the parameter estimates from a fitted Hidden Markov Model (HMM) using a numerical Hessian of the log-likelihood. This function reformulates the HMM as a special case of a Hidden Semi-Markov Model (HSMM) with geometric dwell-time distributions, and calls HSMMVarianceMatrix.

Usage

HMMVarianceMatrix(x, HMM, obsdist, h = 1e-5, method = "central",
                  verbose = TRUE)

Arguments

Details

This function is a wrapper around HSMMVarianceMatrix, treating an HMM as a special case of an HSMM with geometric dwell-time distributions. The Hessian is computed numerically, inverted, and adjusted for positive-definiteness to provide a variance-covariance matrix for the estimated parameters.

Value

A variance-covariance matrix of the HMM parameter estimates, computed as the inverse of the observed information matrix (negative Hessian of the log-likelihood), adjusted to be positive-definite.

Author(s)

Aimee Cody

See Also

generateHMM, for simulating HMM data. HSMMVarianceMatrix, for the general hidden semi-Markov case. confintHMM, for bootstrap confidence intervals of HMM parameters.
findmleHMM, for fitting Hidden Markov Models.

Examples

# Example with 2-state Normal HMM
J <- 2
Pi <- matrix(c(0.9, 0.1,
               0.2, 0.8), nrow = 2, byrow = TRUE)
delta <- c(0.5, 0.5)
obspar <- list(mean = c(0, 3), sd = c(1, 1))
# Simulate data
sim_data <- generateHMM(n = 100, J = J, obsdist = "norm",
                        obspar = obspar, Pi = Pi, delta = delta)
# Fit HMM
HMM_fit <- findmleHMM(x = sim_data$x, J = J,
                      obsdist = "norm", obspar = obspar,
                      Pi = Pi, delta = delta)
# Compute variance-covariance matrix
vcov_matrix <- HMMVarianceMatrix(x = sim_data$x, HMM = HMM_fit,
                                 obsdist = "norm")
vcov_matrix

Variance-Covariance Matrix for Hidden Semi-Markov Models

Description

Computes an approximate variance-covariance matrix of the parameter estimates from a fitted Hidden Semi-Markov Model (HSMM) using a numerical Hessian of the log-likelihood. The Hessian is inverted and adjusted to ensure positive-definiteness, providing an estimate of parameter uncertainty.

Usage

HSMMVarianceMatrix(x, HSMM, M = NA, obsdist, dwelldist,
                   h = 1e-5, method = "central", verbose = TRUE)

Arguments

Details

This function approximates the observed Fisher information for an HSMM by numerically differentiating the log-likelihood with respect to the model parameters. The resulting Hessian is inverted to obtain a variance-covariance matrix. Since numerical Hessians can yield non-positive-definite matrices due to rounding error or flat likelihood surfaces, the eigenvalues are adjusted to enforce positive-definiteness.

Value

A variance-covariance matrix of the HSMM parameter estimates, computed as the inverse of the observed information matrix (negative Hessian of the log-likelihood), adjusted to be positive-definite.

Author(s)

Aimee Cody

See Also

generateHSMM, for simulating HSMM data. findmleHSMM, for fitting Hidden Semi-Markov Models. HMMVarianceMatrix, for numerically computing the variance-covariance matrix of a fitted hidden Markov model.

Examples

# Example with 2-state Normal HSMM
J <- 2
Pi <- matrix(c(0, 1,
               1, 0), nrow = 2, byrow = TRUE)
obspar <- list(mean = c(0, 3), sd = c(1, 1))
dwellpar <- list(lambda = c(5, 7))
# Simulate data
sim_data <- generateHSMM(n = 100, J = J, obsdist = "norm",
                         dwelldist = "pois", obspar = obspar,
                         dwellpar = dwellpar, Pi = Pi)
# Fit HSMM
HSMM_fit <- findmleHSMM(x = sim_data$x, J = J, M = 50,
                        obsdist = "norm", dwelldist = "pois",
                        obspar = obspar, dwellpar = dwellpar, Pi = Pi)
# Compute variance-covariance matrix
vcov_matrix <- HSMMVarianceMatrix(x = sim_data$x, HSMM = HSMM_fit,
                                  obsdist = "norm", dwelldist = "pois",
                                  M = 50)
vcov_matrix

Calculate Information Criteria for HMM/HSMM Models

Description

Computes information criteria (AIC3, AICu, HQC, BICadj) for Hidden Markov Models (HMM) or Hidden Semi-Markov Models (HSMM).

Usage

IC(data, model, type)

Arguments

Details

Automatically detects HMM vs HSMM models based on dwellparameters presence.

Information criteria formulas:

• AIC3: -2 \log L + 3k

• AICu: Unbiased AIC with finite sample correction

• HQC: -2 \log L + 2k \log(\log n)

• BICadj: -2 \log L + k \log n - k \log(2\pi)

where L is likelihood, k is number of parameters, n is sample size.

Value

Numeric value of the specified information criterion.

Author(s)

Aimee Cody

Examples

J <- 3
Pi <- matrix(c(0.0, 0.6, 0.4,
               0.5, 0.0, 0.5,
               0.3, 0.7, 0.0), nrow = J, byrow = TRUE)

obspar <- list(
  mean = c(-2, 0, 3),
  sd   = c(1, 1.5, 2)
)

dwellpar <- list(
  lambda = c(3, 5, 4)
)

sim <- generateHSMM(n = 1000, J = J, obsdist = "norm",
                    dwelldist = "pois", obspar = obspar,
                    dwellpar = dwellpar, Pi = Pi)

fit <- findmleHSMM(x = sim$x, J = J, M = 100,
                   obsdist = "norm", dwelldist = "pois",
                   obspar = obspar, dwellpar = dwellpar,
                   Pi = Pi)

AIC3 <- IC(sim, fit, "AIC3")
AICu <- IC(sim, fit, "AICu")
HQC <- IC(sim, fit, "HQC")
BICadj <- IC(sim, fit, "BICadj")

Plot Conditional Return Levels from GEV-HMM

Description

Computes and plots conditional return levels from a fitted Hidden Markov Model (HMM) with generalized extreme value (GEV) state-dependent distributions. Bootstrap confidence intervals are included. Return levels are aggregated over user-defined time periods.

Usage

conditionalreturnsHMMgev(x, HMM, return_periods = c(50, 100, 200, 500),
                         varcov = NULL, B = 10000, level = 0.95,
                         time_structure = NULL,
                         plot_title = "Return Levels Over Time",
                         save_plot = FALSE, filename = NULL,
                         width = 12, height = 8, dpi = 300,
                         verbose = TRUE)

Arguments

Details

For each state, the conditional return level for a given return period is computed from the fitted GEV parameters. Uncertainty is estimated by parametric bootstrap using draws from a multivariate normal approximation of parameter estimates. Observation-level return levels are aggregated over time periods defined by time_structure (e.g., years, weeks, months).

Value

Invisibly returns a list containing:

A time-series plot is also produced, showing aggregated conditional return levels with bootstrap confidence intervals.

Author(s)

Aimee Cody

See Also

findmleHMM, for fitting HMMs. generateHMM, for simulating HMM data. globaldecodeHMM, for decoding the most likely state sequence. HMMVarianceMatrix, for variance–covariance estimation.

Examples

set.seed(123)
J <- 2
Pi <- matrix(c(0.9, 0.1,
               0.2, 0.8), nrow = J, byrow = TRUE)
delta <- c(0.5, 0.5)
obspar <- list(loc = c(0, 5), scale = c(1, 2), shape = c(0.1, -0.1))

sim <- generateHMM(n = 200, J = J, obsdist = "gev",
                   obspar = obspar, Pi = Pi, delta = delta)

HMM_fit <- findmleHMM(x = sim$x, J = J, obsdist = "gev",
                      obspar = obspar, Pi = Pi, delta = delta)

time_struct <- list(unit = "week", observations_per_unit = 10, start_point = 1)


result <- conditionalreturnsHMMgev(x = sim$x, HMM = HMM_fit,
                                   return_periods = c(50, 100, 200, 500),
                                   B = 1000,
                                   time_structure = time_struct,
                                   plot_title = "GEV-HMM Conditional Return Levels")

result <- conditionalreturnsHMMgev(x = sim$x, HMM = HMM_fit,
                                   return_periods = c(50, 100),
                                   B = 1000,
                                   time_structure = time_struct,
                                   save_plot = TRUE,
                                   filename = tempfile(fileext = ".png"))

Plot Conditional Return Levels from GEV-HSMM

Description

Computes and plots conditional return levels from a fitted Hidden Semi-Markov Model (HSMM) with generalized extreme value (GEV) state-dependent distributions. Bootstrap confidence intervals are included. Return levels are aggregated over user-defined time periods.

Usage

conditionalreturnsHSMMgev(x, HSMM, dwelldist, M = NA,
                          return_periods = c(50, 100, 200, 500),
                          varcov = NULL, B = 10000, level = 0.95,
                          time_structure = NULL, shift = FALSE,
                          plot_title = "Return Levels Over Time",
                          save_plot = FALSE, filename = NULL,
                          width = 12, height = 8, dpi = 300,
                          verbose = TRUE)

Arguments

Details

For each state, the conditional return level for a given return period is computed from the fitted GEV parameters. Uncertainty is estimated by parametric bootstrap using draws from a multivariate normal approximation of parameter estimates. Observation-level return levels are aggregated over time periods defined by time_structure (e.g., years, weeks, months).

Value

Invisibly returns a list containing:

A time-series plot is also produced, showing aggregated conditional return levels with bootstrap confidence intervals.

Author(s)

Aimee Cody

See Also

findmleHSMM, for fitting HSMMs. generateHSMM, for simulating HSMM data. globaldecodeHSMM, for decoding the most likely state sequence. HSMMVarianceMatrix, for variance–covariance estimation.

Examples

set.seed(123)

J <- 2
Pi <- matrix(c(0, 1,
               1, 0), nrow = J, byrow = TRUE)
delta <- c(0.5, 0.5)

obspar <- list(
  loc   = c(0, 5),
  scale = c(1, 2),
  shape = c(0.1, -0.1)
)

dwellpar <- list(lambda = c(5, 10))

sim <- generateHSMM(
  n = 200, J = J, obsdist = "gev", dwelldist = 'pois',
  obspar = obspar, dwellpar = dwellpar, Pi = Pi, delta = delta
)

HSMM_fit <- findmleHSMM(
  x = sim$x, J = J, obsdist = "gev", dwelldist = 'pois',
  obspar = obspar, dwellpar = dwellpar, Pi = Pi, delta = delta
)

time_struct <- list(
  unit = "week", observations_per_unit = 10, start_point = 1
)


result <- conditionalreturnsHSMMgev(
  x = sim$x, HSMM = HSMM_fit,
  return_periods = c(50, 100, 200, 500),
  dwelldist = "pois", B = 1000,
  time_structure = time_struct,
  plot_title = "GEV-HSMM Conditional Return Levels"
)

result <- conditionalreturnsHSMMgev(
  x = sim$x, HSMM = HSMM_fit,
  return_periods = c(50, 100),
  dwelldist = "pois", B = 1000,
  time_structure = time_struct,
  save_plot = TRUE,
  filename = tempfile(fileext = ".png")
)

Bootstrap Confidence Intervals for Hidden Markov Models

Description

Computes bootstrap confidence intervals for the parameters of a fitted Hidden Markov Model (HMM). The function resamples data by simulating from the fitted model, refits the HMM repeatedly to the simulated data, and constructs confidence intervals for the state-dependent parameters and the transition probability matrix.

Usage

confintHMM(x, HMM, obsdist, level = 0.95, B = 100, EM = FALSE,
           verbose = TRUE, seed = NULL, ...)

Arguments

Details

The function uses a parametric bootstrap procedure. Data are simulated from the fitted HMM using generateHMM, and parameters are re-estimated with findmleHMM. Confidence intervals are then obtained from the empirical quantiles of the bootstrap replicates. When verbose = TRUE, the function reports progress. The function will stop if too many failed bootstrap attempts occur.

Value

A list with three elements:

Author(s)

Aimee Cody

See Also

findmleHMM, for estimating HMM parameters from data. generateHMM, for simulating HMM data. confintHSMM, computing bootstrapped confidence intervals for a hidden semi-Markov model.

Examples

J <- 2
Pi <- matrix(c(0.9, 0.1,
               0.2, 0.8), nrow = 2, byrow = TRUE)
obspar <- list(lambda = c(2, 7))
sim_data <- generateHMM(n = 1000, J = J, obsdist = "pois", obspar = obspar, Pi = Pi)
HMM_fit <- findmleHMM(J = J, x = sim_data$x, obsdist = "pois",
                      obspar = obspar, Pi = Pi, EM = FALSE)

ci <- confintHMM(x = sim_data$x, HMM = HMM_fit, obsdist = "pois",
                 B = 50, level = 0.9)
ci$obspar_CI
ci$Pi_CI
ci_silent <- confintHMM(x = sim_data$x, HMM = HMM_fit, obsdist = "pois",
                        B = 50, level = 0.9, verbose = FALSE)

Bootstrap Confidence Intervals for Hidden Semi-Markov Models

Description

Computes bootstrap confidence intervals for the parameters of a fitted Hidden Semi-Markov Model (HSMM). The function resamples data by simulation from the fitted model, refits the HSMM repeatedly, and constructs confidence intervals for the state-dependent observation parameters, dwell-time distribution parameters, and the transition probability matrix.

Usage

confintHSMM(x, HSMM, obsdist, dwelldist, level = 0.95, B = 100, M = NA,
            shift = FALSE, maxiter = 100, tol = 1e-5, verbose = TRUE, seed = NULL)

Arguments

Details

The function uses a parametric bootstrap procedure. Data are simulated from the fitted HSMM using generateHSMM, and parameters are re-estimated with findmleHSMM. Confidence intervals are then obtained from the empirical quantiles of the bootstrap replicates. When verbose = TRUE, progress is reported. Replicates that fail to converge are retried until B successful fits are obtained.

Value

A list with four elements:

Author(s)

Aimee Cody

See Also

findmleHSMM, for estimating HSMM parameters from data. generateHSMM, for simulating HSMM data. confintHMM, computing bootstrapped confidence intervals for a hidden Markov model.

Examples

J <- 2
Pi <- matrix(c(0, 1,
               1, 0), nrow = 2, byrow = TRUE)

obspar <- list(mean = c(4, 3), sd = c(1, 1))

dwellpar <- list(lambda = c(5, 7))

sim_data <- generateHSMM(n = 1000, J = J, obsdist = "norm",
                         dwelldist = "pois", obspar = obspar,
                         dwellpar = dwellpar, Pi = Pi)

HSMM_fit <- findmleHSMM(x = sim_data$x, J = J, M = 100,
                        obsdist = "norm", dwelldist = "pois",
                        obspar = obspar, dwellpar = dwellpar, Pi = Pi)


ci <- confintHSMM(x = sim_data$x, HSMM = HSMM_fit,
                  obsdist = "norm", dwelldist = "pois",
                  B = 50, M = 50, level = 0.9)

ci$obspar_CI
ci$dwellpar_CI
ci$Pi_CI

ci_silent <- confintHSMM(x = sim_data$x, HSMM = HSMM_fit,
                         obsdist = "norm", dwelldist = "pois",
                         B = 50, M = 50, level = 0.9, verbose = FALSE)

Daily maxima of geomagnetic rate-of-change data from Eskdalemuir Magnetic Observatory

Description

A dataset containing daily block maxima of the horizontal geomagnetic field rate-of-change at the Eskdalemuir Magnetic Observatory, Scotland, United Kingdom. The data reflect the daily extremes of geomagnetic variability, serving as indicators of geomagnetic disturbance over the period 1999–2022.

Usage

data(daily_maxima)

Format

A numeric vector of daily maxima of

R_1(t) = \sqrt{(X_t - X_{t-1})^2 + (Y_t - Y_{t-1})^2},

where X_t and Y_t denote the horizontal north and east components of the magnetic field, respectively.

Details

These daily maxima capture short-term geomagnetic disturbances while maintaining a sufficiently fine temporal resolution to observe rapid changes. Like the weekly maxima, they are naturally modelled using the Generalised Extreme Value (GEV) distribution, with location, scale, and shape parameters \mu, \sigma, and \xi.

Source

Eskdalemuir Magnetic Observatory, British Geological Survey, United Kingdom.

Derived and aggregated for use in HMMHSMM.

Examples

data(daily_maxima)
plot(daily_maxima, type = "l", main = "Daily maxima of geomagnetic variability",
     ylab = "Daily maximum rate of change (nT/min)", xlab = "Day index")

Plot Exceedance Probabilities from GEV-HMM

Description

Computes and plots exceedance probabilities of a given threshold from a fitted Hidden Markov Model (HMM) with generalized extreme value (GEV) state-dependent distributions. Bootstrap confidence intervals are included. Probabilities are aggregated over user-defined time periods.

Usage

exceedanceplotHMMgev(x, HMM, threshold, varcov = NULL, B = 10000, level = 0.95,
                     time_structure = NULL,
                     plot_title = "Exceedance Probabilities Over Time",
                     save_plot = FALSE, filename = NULL,
                     width = 12, height = 8, dpi = 300,
                     verbose = TRUE)

Arguments

Details

For each state, the exceedance probability of threshold is computed from the fitted GEV parameters. Uncertainty is estimated by parametric bootstrap using draws from a multivariate normal approximation of parameter estimates. Observation-level probabilities are aggregated over time periods defined by time_structure (e.g., years, months, weeks). When verbose = TRUE, progress messages are displayed during computation.

Value

Invisibly returns a list containing:

A time-series plot is also produced, showing aggregated exceedance probabilities with confidence intervals.

Author(s)

Aimee Cody

See Also

findmleHMM, for fitting HMMs. generateHMM, for simulating HMM data. globaldecodeHMM, for decoding the most likely state sequence. HMMVarianceMatrix, for variance–covariance estimation.

Examples

set.seed(99)
J <- 2
Pi <- matrix(c(0.9, 0.1,
               0.2, 0.8), nrow = J, byrow = TRUE)
delta <- c(0.5, 0.5)
obspar <- list(loc = c(0, 5), scale = c(1, 2), shape = c(0.1, -0.1))

sim <- generateHMM(n = 200, J = J, obsdist = "gev",
                   obspar = obspar, Pi = Pi, delta = delta)

HMM_fit <- findmleHMM(x = sim$x, J = J, obsdist = "gev",
                      obspar = obspar, Pi = Pi, delta = delta)

time_struct <- list(unit = "week", observations_per_unit = 10, start_point = 1)


result <- exceedanceplotHMMgev(x = sim$x, HMM = HMM_fit,
                               threshold = 20, B = 10000,
                               time_structure = time_struct,
                               plot_title = "GEV-HMM Exceedance Probabilities")

result <- exceedanceplotHMMgev(x = sim$x, HMM = HMM_fit,
                               threshold = 20, B = 10000,
                               time_structure = time_struct,
                               save_plot = TRUE,
                               filename = tempfile(fileext = ".png"),
                               verbose = FALSE)

Plot Exceedance Probabilities from GEV-HSMM

Description

Computes and plots exceedance probabilities of a given threshold from a fitted Hidden Semi-Markov Model (HSMM) with generalized extreme value (GEV) state-dependent distributions. Bootstrap confidence intervals are included. Probabilities are aggregated over user-defined time periods.

Usage

exceedanceplotHSMMgev(x, HSMM, threshold, dwelldist, M = NA,
                      varcov = NULL, B = 10000, level = 0.95,
                      time_structure = NULL, shift = FALSE,
                      plot_title = "Exceedance Probabilities Over Time",
                      save_plot = FALSE, filename = NULL,
                      width = 12, height = 8, dpi = 300,
                      verbose = TRUE)

Arguments

Details

For each state, the exceedance probability of threshold is computed from the fitted GEV parameters. Uncertainty is estimated by parametric bootstrap using draws from a multivariate normal approximation of parameter estimates. Observation-level probabilities are aggregated over time periods defined by time_structure (e.g., years, months, weeks). When verbose = TRUE, progress messages are displayed during computation.

Value

Invisibly returns a list containing:

A time-series plot is also produced, showing aggregated exceedance probabilities with confidence intervals.

Author(s)

Aimee Cody

See Also

findmleHSMM, for fitting HSMMs. generateHSMM, for simulating HSMM data. globaldecodeHSMM, for decoding the most likely state sequence. HSMMVarianceMatrix, for variance–covariance estimation.

Examples

set.seed(123)

J <- 2
Pi <- matrix(c(0, 1,
               1, 0), nrow = J, byrow = TRUE)
delta <- c(0.5, 0.5)

obspar <- list(
  loc   = c(0, 5),
  scale = c(1, 2),
  shape = c(0.1, -0.1)
)

dwellpar <- list(lambda = c(5, 10))

sim <- generateHSMM(
  n = 200, J = J, obsdist = "gev", dwelldist = "pois",
  obspar = obspar, dwellpar = dwellpar, Pi = Pi, delta = delta
)

HSMM_fit <- findmleHSMM(
  x = sim$x, J = J, obsdist = "gev", dwelldist = "pois",
  obspar = obspar, dwellpar = dwellpar, Pi = Pi, delta = delta
)

time_struct <- list(
  unit = "week", observations_per_unit = 10, start_point = 1
)


result <- exceedanceplotHSMMgev(
  x = sim$x, HSMM = HSMM_fit,
  threshold = 20, dwelldist = "pois",
  B = 10000, time_structure = time_struct,
  plot_title = "GEV-HSMM Exceedance Probabilities"
)

result <- exceedanceplotHSMMgev(
  x = sim$x, HSMM = HSMM_fit,
  threshold = 20, dwelldist = "pois",
  B = 10000, time_structure = time_struct,
  save_plot = TRUE,
  filename = tempfile(fileext = ".png"),
  verbose = FALSE
)

Maximum Likelihood Estimation for Hidden Markov Models

Description

Estimates the parameters of a Hidden Markov Model (HMM) using either the EM-based semi-Markov parameterisation or direct numerical maximization of the likelihood.

Usage

findmleHMM(J, x, obsdist, obspar, Pi, EM = FALSE, verbose = TRUE, seed = NULL, ...)

Arguments

Details

This function fits a Hidden Markov Model to a sequence of observations using either:

• An EM-based approximation via a semi-Markov model when EM = TRUE

• Direct numerical maximization of the likelihood using nlm when EM = FALSE

The EM-based approach initializes dwell-time probabilities from the diagonal of Pi and estimates the HMM via the corresponding semi-Markov model. After fitting, transition probabilities are adjusted to incorporate self-transition probabilities. Supported observation distributions include normal, Poisson, Weibull, zero-inflated Poisson (ZIP), negative binomial, zero-inflated negative binomial (ZINB), exponential, gamma, log-normal, generalized extreme value (GEV), zero-inflated normal, and zero-inflated gamma. Parameter requirements vary by distribution. When verbose = TRUE and EM = TRUE, progress messages from the EM algorithm are displayed.

Value

A list containing:

Author(s)

Aimee Cody

See Also

generateHMM for simulating HMM data. findmleHSMM for EM-based semi-Markov estimation.

Examples

J <- 3
Pi <- matrix(c(0.8, 0.15, 0.05,
               0.1, 0.7, 0.2,
               0.2, 0.3, 0.5), nrow = 3, byrow = TRUE)

obspar <- list(
  mean = c(-2, 0, 3),
  sd = c(0.5, 1, 1.5)
)

x <- generateHMM(n = 200, J = J, obsdist = "norm",
                       obspar = obspar, Pi = Pi)$x

fit <- findmleHMM(J = J, x = x, obsdist = "norm",
                   obspar = obspar, Pi = Pi, EM = FALSE)

fit$estimate
fit$loglik
fit$AIC

fit_em <- findmleHMM(J = J, x = x, obsdist = "norm",
                      obspar = obspar, Pi = Pi, EM = TRUE, verbose = FALSE)

Multiple Initialization Maximum Likelihood Estimation for Hidden Markov Models

Description

Fits a Hidden Markov Model (HMM) by repeatedly initializing observation and transition parameters and selecting the fit with the highest log-likelihood. This approach helps avoid convergence to poor local optima. For the generalized extreme value (GEV) distribution, starting values are generated from repeated maximum likelihood fits on random data subsets.

Usage

findmleHMMnostarting(J, x, obsdist, no.initials = 50, EM = FALSE,
                     verbose = TRUE, seed = NULL, ...)

Arguments

Details

This function automates multiple trials of findmleHMM with randomized starting values, returning the fit that achieves the highest log-likelihood.

• For most observation distributions, starting values are generated via clusterHMM.

• For the GEV distribution, starting values are drawn from repeated fits of evd::fgev on random data segments. Up to 20,000 attempts are made, and a warning is issued if fewer than 1000 valid estimates are obtained.

During each iteration:

1. Observation parameters are perturbed slightly to encourage exploration.

2. A transition matrix Pi is drawn from a random uniform distribution with added self-transition bias.

3. The HMM is estimated via findmleHMM.

4. If the resulting log-likelihood exceeds the current best, the model is updated.

At the end of all iterations, the best-fitting model is returned. When verbose = TRUE, iteration numbers and error messages are displayed during the fitting process.

Value

A list corresponding to the best fit across all initializations, containing:

Author(s)

Aimee Cody

See Also

findmleHMM for fitting an HMM with user-supplied starting values. generateHMM for simulating HMM data. findmleHSMMnostarting for the non-initialised estimation of hidden semi-Markov models.

Examples

set.seed(123)
J <- 3
Pi <- matrix(c(0.7, 0.2, 0.1,
               0.1, 0.8, 0.1,
               0.2, 0.3, 0.5), nrow = 3, byrow = TRUE)
obspar <- list(mean = c(-2, 0, 3),
               sd   = c(0.5, 1, 1.5))
x <- generateHMM(n = 200, J = J, Pi = Pi, obsdist = "norm", obspar = obspar)$x


fit <- findmleHMMnostarting(J = J, x = x, obsdist = "norm",
                            no.initials = 30)

fit$loglik
fit$estimate

fit_silent <- findmleHMMnostarting(J = J, x = x, obsdist = "norm",
                                   no.initials = 30, verbose = FALSE)

Maximum Likelihood Estimation for Hidden Semi-Markov Models

Description

Estimates the parameters of a Hidden Semi-Markov Model (HSMM) using the EM algorithm with specified observation and dwell-time distributions.

Usage

findmleHSMM(x, J, M = NA, obsdist, dwelldist, obspar, dwellpar,
            Pi, delta = NULL, maxiter = 100, tol = 1e-05, shift = FALSE,
            verbose = TRUE, seed = NULL)

Arguments

Details

This function fits a Hidden Semi-Markov Model to a sequence of observations using an EM algorithm. At each iteration, the algorithm:

• Computes observation probabilities for each state

• Calculates dwell-time probabilities and survival functions

• Performs a forward-backward-like procedure to estimate expected sufficient statistics

• Re-estimates observation and dwell-time parameters

Supported observation distributions include normal, Poisson, Weibull, zero-inflated Poisson (ZIP), negative binomial, zero-inflated negative binomial (ZINB), exponential, gamma, log-normal, generalized extreme value (GEV), zero-inflated normal, and zero-inflated gamma.

Supported dwell-time distributions include Poisson, negative binomial, beta-binomial, and geometric.

When verbose = TRUE, the function displays log-likelihood values at each iteration and reports whether convergence was achieved.

Value

A list containing:

Author(s)

Aimee Cody

See Also

findmleHMM for Hidden Markov Model estimation. generateHSMM for simulating HSMM data.

Examples

J <- 3
Pi <- matrix(c(0.0, 0.6, 0.4,
               0.5, 0.0, 0.5,
               0.3, 0.7, 0.0), nrow = J, byrow = TRUE)

obspar <- list(
  mean = c(-2, 0, 3),
  sd   = c(1, 1.5, 2)
)

dwellpar <- list(
  lambda = c(3, 5, 4)
)

sim <- generateHSMM(n = 1000, J = J, obsdist = "norm",
                    dwelldist = "pois", obspar = obspar,
                    dwellpar = dwellpar, Pi = Pi)

fit <- findmleHSMM(x = sim$x, J = J, M = 100,
                   obsdist = "norm", dwelldist = "pois",
                   obspar = obspar, dwellpar = dwellpar,
                   Pi = Pi)

fit$observationparameters
fit$dwellparameters
fit$Pi
fit$delta
fit$loglikelihoods

fit_silent <- findmleHSMM(x = sim$x, J = J, M = 100,
                          obsdist = "norm", dwelldist = "pois",
                          obspar = obspar, dwellpar = dwellpar,
                          Pi = Pi, verbose = FALSE)

Fit Hidden Semi-Markov Model (HSMM) Without User-Provided Starting Values

Description

Fits a hidden semi-Markov model (HSMM) to a univariate time series by maximum likelihood, using multiple randomly generated initializations instead of user-provided starting parameters. The best-fitting model is selected by log-likelihood.

Usage

findmleHSMMnostarting(J, x, obsdist, dwelldist, M = NA, no.initials = 50,
                      verbose = TRUE, seed = NULL, ...)

Arguments

Details

This function automates parameter initialization by sampling plausible random values for state-dependent observation parameters, dwell-time parameters, initial distribution \delta, and transition matrix \Pi. It repeatedly calls findmleHSMM with these initializations, and retains the model with the highest log-likelihood. When obsdist = "gev", initial values are obtained via repeated calls to evd::fgev on random data segments until a sufficient number of estimates are collected.

Value

An object containing:

Author(s)

Aimee Cody

See Also

findmleHSMM, for fitting an HSMM with user-supplied initial values. generateHSMM, for simulating HSMM data. findmleHMMnostarting, for the non-initialised estimation of hidden Markov models.

Examples

set.seed(321)
J <- 3
Pi <- matrix(c(0.0, 0.6, 0.4,
               0.5, 0.0, 0.5,
               0.3, 0.7, 0.0), nrow = J, byrow = TRUE)
obspar <- list(
  mean = c(-2, 0, 3),
  sd   = c(1, 1.5, 2)
)
dwellpar <- list(
  lambda = c(3, 5, 4)
)
sim <- generateHSMM(n = 200, J = J, obsdist = "norm",
                    dwelldist = "pois", obspar = obspar,
                    dwellpar = dwellpar, Pi = Pi)
fit <- findmleHSMMnostarting(J = J, x = sim$x,
                                       obsdist = "norm", dwelldist = "pois",
                                       M = 100, no.initials = 30)
fit$observationparameters
fit$dwellparameters
fit$Pi
fit$delta
fit$loglikelihoods

Generate Data from a Hidden Markov Model

Description

Simulates observations and hidden states from a Hidden Markov Model (HMM) with specified observation distributions.

Usage

generateHMM(n, J, obsdist, obspar, Pi, delta = NULL, seed = NULL)

Arguments

Details

This function simulates data from a Hidden Markov Model where:

• Hidden states follow a discrete-time Markov chain with transition matrix Pi

• At each time step, observations are generated from state-dependent distributions

The function supports multiple observation distributions including normal, Poisson, Weibull, zero-inflated Poisson (ZIP), negative binomial, zero-inflated negative binomial (ZINB), exponential, gamma, log-normal, generalized extreme value (GEV), zero-inflated normal, and zero-inflated gamma. The initial state is sampled according to the initial distribution delta. If delta is not provided, it is computed as the stationary distribution of the Markov chain defined by Pi.

Value

A list containing:

Author(s)

[Aimee Cody]

See Also

generateHSMM for Hidden Semi-Markov Models.

Examples

# Example with 3 states, normal observations
J <- 3
# HMM transition matrix
Pi <- matrix(c(0.8, 0.15, 0.05,
               0.1, 0.7,  0.2,
               0.2, 0.3,  0.5), nrow = 3, byrow = TRUE)
# Observation parameters (normal distribution)
obspar <- list(
  mean = c(-2, 0, 3),
  sd = c(0.5, 1, 1.5)
)
# Generate 200 observations
sim_data <- generateHMM(n = 200, J = J, obsdist = "norm",
                       obspar = obspar, Pi = Pi)
# View the results
head(sim_data$x)      # observations
head(sim_data$state)  # hidden states
length(sim_data$x)    # number of observations

Generate Data from a Hidden Semi-Markov Model

Description

Simulates observations and hidden states from a Hidden Semi-Markov Model (HSMM) with specified observation and dwell time distributions.

Usage

generateHSMM(n, J, obsdist, dwelldist, obspar, dwellpar, Pi,
             delta = NULL, simtype = "nobs", shift = FALSE, seed = NULL)

Arguments

Details

This function simulates data from a Hidden Semi-Markov Model where:

• Hidden states follow a Markov chain with transition matrix Pi

• Each state has an associated dwell time distribution that determines how long the process remains in that state

• Observations are generated from state-dependent distributions

The function supports multiple observation distributions including normal, Poisson, Weibull, zero-inflated Poisson (ZIP), negative binomial, zero-inflated negative binomial (ZINB), exponential, gamma, log-normal, generalized extreme value (GEV), zero-inflated normal and zero-inflated gamma.

Dwell time distributions include Poisson, negative binomial, and beta-binomial, all with optional shift parameters to ensure minimum dwell times.

When simtype = "nobs", the function generates exactly n observations. When simtype = "nseq", it generates n state sequences and the total number of observations depends on the realized dwell times.

Value

A list containing:

Author(s)

[Aimee Cody]

See Also

generateHMM for Hidden Markov Models.

Examples

# Example with 3 states, normal observations, Poisson dwell times
J <- 3
# HSMM transition matrix
Pi <- matrix(c(0.0, 0.6, 0.4,
               0.5, 0.0, 0.5,
               0.3, 0.7, 0.0), nrow = 3, byrow = TRUE)

# Observation parameters (normal distribution)
obspar <- list(
  mean = c(-2, 0, 3),
  sd = c(1, 1.5, 2)
)

# Dwell time parameters (Poisson distribution)
dwellpar <- list(
  lambda = c(3, 5, 4)
)

# Generate 100 observations
sim_data <- generateHSMM(n = 100, J = J, obsdist = "norm", dwelldist = "pois",
                        obspar = obspar, dwellpar = dwellpar, Pi = Pi)

# View the results
head(sim_data$x)        # observations
head(sim_data$states)   # hidden states
sim_data$N              # total number of observations

Global Decoding for Hidden Markov Models

Description

Computes the most likely hidden state sequence for a univariate time series using a fitted Hidden Markov Model (HMM) via global decoding (Viterbi-style).

Usage

globaldecodeHMM(x, HMM, obsdist)

Arguments

Details

This function implements global decoding for an HMM using the estimated parameters from a fitted model. It calculates the forward probabilities and recursively determines the most probable sequence of hidden states. Global decoding provides the single most likely state sequence, unlike local decoding which gives the most probable state at each time point independently.

Value

Author(s)

Aimee Cody

See Also

findmleHMM, for estimating HMM parameters from data. generateHMM, for simulating HMM data. localdecodeHMM, for local decoding of an HMM. globaldecodeHSMM, for global decoding of a hidden semi-Markov model.

Examples

# Example with 3 states, normal observations
J <- 3

# HMM transition matrix
Pi <- matrix(c(0.8, 0.15, 0.05,
               0.1, 0.7, 0.2,
               0.2, 0.3, 0.5), nrow = 3, byrow = TRUE)

# Observation parameters (normal distribution)
obspar <- list(
  mean = c(-2, 0, 3),
  sd = c(0.5, 1, 1.5)
)

# Simulate HMM data
sim_data <- generateHMM(n = 200, J = J, obsdist = "norm", obspar = obspar, Pi = Pi)

# Fit HMM to simulated data
HMM_fit <- findmleHMM(J = J, x = sim_data$x, obsdist = "norm",
                       obspar = obspar, Pi = Pi, EM = FALSE)

# Compute global decoding
decoded_states <- globaldecodeHMM(x = sim_data$x, HMM = HMM_fit, obsdist = "norm")

decoded_states

Global Decoding of Hidden Semi-Markov Models

Description

Computes the most likely hidden state sequence for a univariate time series using a fitted Hidden Semi-Markov Model (HSMM) via the Viterbi algorithm.

Usage

globaldecodeHSMM(x, M = NA, HSMM, obsdist, dwelldist, shift = FALSE)

Arguments

Details

This function computes the global most probable state sequence (Viterbi path) for a sequence of observations under a Hidden Semi-Markov Model using the parameters contained in a fitted HSMM object. The function calculates observation log-probabilities, dwell-time log-probabilities, and forward log-probabilities to determine the Viterbi path.

Value

A numeric vector of length length(x), giving the most likely hidden state at each observation.

Author(s)

Aimee Cody

See Also

findmleHSMM, for estimating HSMM parameters from data. generateHSMM, for simulating HSMM data. localdecodeHSMM, for local decoding of HSMM states. globaldecodeHMM, for global decoding hidden Markov models.

Examples

J <- 3
Pi <- matrix(c(0.0, 0.6, 0.4,
               0.5, 0.0, 0.5,
               0.3, 0.7, 0.0), nrow = J, byrow = TRUE)

obspar <- list(
  mean = c(4, 2, 6),
  sd   = c(1, 1, 2)
)

dwellpar <- list(
  lambda = c(3, 5, 4)
)

# Simulate HSMM data
sim <- generateHSMM(n = 1000, J = J, obsdist = "norm",
                    dwelldist = "pois", obspar = obspar,
                    dwellpar = dwellpar, Pi = Pi)

# Fit HSMM using the true parameters
HSMM_true <- findmleHSMM(x = sim$x, J = J, M = 100,
                          obsdist = "norm", dwelldist = "pois",
                          obspar = obspar, dwellpar = dwellpar,
                          Pi = Pi)

# Decode states using globaldecodeHSMM
states <- globaldecodeHSMM(x = sim$x, M = 100, HSMM = HSMM_true,
                           obsdist = "norm", dwelldist = "pois")

head(states)

Local Decoding for Hidden Markov Models

Description

Computes the most likely hidden state sequence for a univariate time series using a fitted Hidden Markov Model (HMM) object via the forward-backward algorithm.

Usage

localdecodeHMM(x, HMM, obsdist)

Arguments

Details

This function uses the forward-backward algorithm to compute local state probabilities from a fitted Hidden Markov Model (HMM). The most probable state at each observation is returned. Initial state probabilities are derived from the stationary distribution of the transition matrix Pi if not explicitly provided.

Supported observation distributions include normal, Poisson, Weibull, zero-inflated Poisson (ZIP), negative binomial, zero-inflated negative binomial (ZINB), exponential, gamma, log-normal, generalized extreme value (GEV), zero-inflated normal, and zero-inflated gamma.

Value

Numeric vector of length length(x), containing the most likely hidden state at each observation, computed via local decoding using the forward-backward algorithm.

Author(s)

Aimee Cody

See Also

findmleHMM, for estimating HMM parameters from data. generateHMM, for simulating HMM data. localdecodeHSMM, for local decoding of Hidden Semi-Markov Models. globaldecodeHMM, for global decoding in HMMs.

Examples

# Example with 3 states, normal observations
J <- 3

# HMM transition matrix
Pi <- matrix(c(0.8, 0.15, 0.05,
               0.1, 0.7, 0.2,
               0.2, 0.3, 0.5), nrow = 3, byrow = TRUE)

# Observation parameters (normal distribution)
obspar <- list(
  mean = c(-2, 0, 3),
  sd = c(0.5, 1, 1.5)
)

# Simulate HMM data
sim_data <- generateHMM(n = 200, J = J, obsdist = "norm", obspar = obspar, Pi = Pi)

# Fit HMM to simulated data
HMM_fit <- findmleHMM(J = J, x = sim_data$x, obsdist = "norm",
                       obspar = obspar, Pi = Pi, EM = FALSE)

# Compute local decoding
decoded_states <- localdecodeHMM(x = sim_data$x, HMM = HMM_fit, obsdist = "norm")

decoded_states

Local Decoding for Hidden Semi-Markov Models

Description

Computes the most likely hidden state sequence for a univariate time series using a fitted Hidden Semi-Markov Model (HSMM).

Usage

localdecodeHSMM(x, M = NA, HSMM, obsdist, dwelldist)

Arguments

Details

This function uses the forward-backward algorithm for HSMMs to compute local state probabilities and returns the most probable state at each time point. It supports a wide range of observation distributions and dwell-time distributions. The user can supply a maximum dwell-time truncation parameter M and optionally use a shift parameter for dwell times.

Value

Author(s)

Aimee Cody

See Also

findmleHSMM, for estimating HSMM parameters from data. generateHSMM, for simulating HSMM data. localdecodeHMM, for local decoding of a hidden Markov model. globaldecodeHSMM, for finding the most likely sequence of states using global decoding.

Examples

J <- 3
Pi <- matrix(c(0.0, 0.6, 0.4,
               0.5, 0.0, 0.5,
               0.3, 0.7, 0.0), nrow = J, byrow = TRUE)

obspar <- list(
  mean = c(-4, 0, 4),
  sd   = c(1, 0.5, 1.5)
)

dwellpar <- list(
  lambda = c(3, 5, 4)
)

# Simulate HSMM data
sim <- generateHSMM(n = 1000, J = J, obsdist = "norm",
                    dwelldist = "pois", obspar = obspar,
                    dwellpar = dwellpar, Pi = Pi)

# Fit HSMM using the true parameters
HSMM_true <- findmleHSMM(x = sim$x, J = J, M = 100,
                          obsdist = "norm", dwelldist = "pois",
                          obspar = obspar, dwellpar = dwellpar,
                          Pi = Pi)

# Decode states using localdecodeHSMM
decoded <- localdecodeHSMM(x = sim$x, HSMM = HSMM_true,
                            obsdist = "norm", dwelldist = "pois",
                            M = 100)
decoded

Plot Hidden Markov Model Parameters Over Time

Description

Plots the estimated state-dependent parameters of a fitted Hidden Markov Model (HMM) over time. Confidence intervals can be included, along with optional overlay of observed data for visual comparison.

Usage

plotHMMParameters(x, HMM, obsdist, confint_result = NULL,
                  level = 0.95, B = 100,
                  time_structure = NULL,
                  plot_title = "HMM Parameters Over Time",
                  overlay_data = NULL, overlay_label = "Data",
                  colors = c("black", "red", "blue", "green"),
                  save_plot = FALSE, filename = NULL,
                  width = 12, height = 8, dpi = 300,
                  verbose = TRUE, seed=NULL)

Arguments

Details

For each observation, the most likely hidden state is decoded using globaldecodeHMM. The corresponding state-dependent parameters are extracted and plotted over time. If confint_result is not supplied, bootstrap confidence intervals are computed via confintHMM. Custom time scaling can be applied using the time_structure argument. Overlay data (such as covariates or raw series) can be plotted on a secondary axis for comparison. When verbose = TRUE, progress messages are displayed during confidence interval computation.

Value

Invisibly returns a list containing:

The function also produces parameter trajectory plots with confidence intervals and optional overlays.

Author(s)

Aimee Cody

See Also

findmleHMM, for fitting HMMs. generateHMM, for simulating HMM data. globaldecodeHMM, for decoding the most likely state sequence. confintHMM, for bootstrap confidence intervals of HMM parameters.

Examples

set.seed(123)
J <- 2
Pi <- matrix(c(0.9, 0.1,
               0.2, 0.8), nrow = 2, byrow = TRUE)
delta <- c(0.5, 0.5)
obspar <- list(lambda = c(2, 7))

sim_data <- generateHMM(n = 120, J = J, obsdist = "pois",
                        obspar = obspar, Pi = Pi, delta = delta)

HMM_fit <- findmleHMM(x = sim_data$x, J = J, obsdist = "pois",
                      obspar = obspar, Pi = Pi, delta = delta)

overlay_series <- sim_data$x

time_struct <- list(unit = "day", observations_per_unit = 12, start_point = 1)


result <- plotHMMParameters(x = sim_data$x, HMM = HMM_fit, obsdist = "pois",
                            overlay_data = overlay_series,
                            overlay_label = "Observed counts",
                            time_structure = time_struct,
                            plot_title = "Poisson HMM Parameters with Overlay")

result <- plotHMMParameters(x = sim_data$x, HMM = HMM_fit, obsdist = "pois",
                            overlay_data = overlay_series,
                            overlay_label = "Observed counts",
                            time_structure = time_struct,
                            save_plot = TRUE,
                            filename = tempfile(fileext = ".png"),
                            verbose = FALSE)

Plot Hidden Semi-Markov Model Parameters Over Time

Description

Plots state-dependent observation parameters and dwell-time parameters of a fitted Hidden Semi-Markov Model (HSMM) over time. Confidence intervals can be included, with optional overlay of observed data for comparison.

Usage

plotHSMMParameters(x, HSMM, obsdist, dwelldist, confint_result = NULL,
                   level = 0.95, B = 100, M = NA, include_dwell = FALSE,
                   shift = FALSE, time_structure = NULL,
                   plot_title = "HSMM Parameters Over Time",
                   overlay_data = NULL, overlay_label = "Data",
                   colors = c("black", "red", "blue", "darkgreen"),
                   save_plot = FALSE, filename = NULL,
                   width = 12, height = 8, dpi = 300,
                   verbose = TRUE, seed = NULL)

Arguments

Details

The most likely hidden state sequence is decoded using globaldecodeHSMM. For each decoded state, observation and dwell-time parameters are extracted and plotted over time. If confint_result is not provided, bootstrap intervals are computed via confintHSMM. Plots include observation parameter trajectories, dwell-time parameter trajectories (if include_dwell = TRUE), and optional overlay data on secondary axes. Time scaling is configurable using time_structure. When verbose = TRUE, progress messages are displayed during confidence interval computation. If save_plot = TRUE, the plot is saved as a PNG file to the specified filename with the given dimensions and resolution.

Value

Invisibly returns a list containing:

The function also generates multi-panel plots of parameter trajectories with confidence intervals and optional overlays.

Author(s)

Aimee Cody

See Also

findmleHSMM, for fitting HSMMs. generateHSMM, for simulating HSMM data. globaldecodeHSMM, for decoding the most likely state sequence. confintHSMM, for bootstrap confidence intervals of HSMM parameters.

Examples

set.seed(42)
J <- 2
Pi <- matrix(c(0.0, 1.0,
               1.0, 0.0), nrow = J, byrow = TRUE)

obspar <- list(mean = c(0, 3), sd = c(1, 1.2))
dwellpar <- list(lambda = c(5, 8))

sim <- generateHSMM(n = 200, J = J,
                    obsdist = "norm", dwelldist = "pois",
                    obspar = obspar, dwellpar = dwellpar,
                    Pi = Pi)

HSMM_fit <- findmleHSMM(x = sim$x, J = J, M = 100,
                        obsdist = "norm", dwelldist = "pois",
                        obspar = obspar, dwellpar = dwellpar,
                        Pi = Pi)

overlay_series <- sim$x
time_struct <- list(unit = "day", observations_per_unit = 10, start_point = 1)


result <- plotHSMMParameters(x = sim$x, HSMM = HSMM_fit,
                             obsdist = "norm", dwelldist = "pois",
                             overlay_data = overlay_series,
                             overlay_label = "Observed values",
                             time_structure = time_struct,
                             plot_title = "HSMM Parameters with Overlay")

result_silent <- plotHSMMParameters(x = sim$x, HSMM = HSMM_fit,
                                    obsdist = "norm", dwelldist = "pois",
                                    include_dwell = TRUE,
                                    overlay_data = overlay_series,
                                    time_structure = time_struct,
                                    verbose = FALSE)

Ordinary Residuals for Hidden Markov Models

Description

Computes ordinary (probability integral transform) residuals for a fitted Hidden Markov Model (HMM). Generates quantile-quantile and autocorrelation plots for assessing model fit.

Usage

residualsHMM(x, HMM, obsdist, lag.max = 50, verbose = TRUE)

Arguments

Details

The function computes probability integral transform (PIT) residuals for HMM observations by combining forward and backward probabilities with the cumulative distribution functions of the state-dependent observations. The residuals are assessed via:

• QQ plots comparing observed residuals to the expected standard normal distribution.

• Autocorrelation function plots to detect remaining temporal dependence.

This provides a visual check of model adequacy, where ideally residuals should behave as independent standard normal variables. When verbose = TRUE, the proportion of residuals outside the 95% interval is printed.

Value

Invisibly returns a list containing:

The function also produces diagnostic plots:

• Generates a QQ plot of ordinary residuals against expected standard normal quantiles, with 95% confidence intervals.

• Generates an ACF plot of the ordinary residuals with approximate 95% confidence bands.

Author(s)

Aimee Cody

See Also

findmleHMM, for fitting HMMs. generateHMM, for simulating HMM data. residualsHSMM, for residual analysis of hidden semi-Markov models.

Examples

J <- 2
Pi <- matrix(c(0.9, 0.1,
               0.2, 0.8), nrow = 2, byrow = TRUE)
delta <- c(0.5, 0.5)
obspar <- list(lambda = c(2, 7))

sim_data <- generateHMM(n = 200, J = J, obsdist = "pois", obspar = obspar, Pi = Pi, delta = delta)

HMM_fit <- findmleHMM(x = sim_data$x, J = J, obsdist = "pois",
                      obspar = obspar, Pi = Pi, delta = delta)

result <- residualsHMM(x = sim_data$x, HMM = HMM_fit, obsdist = "pois", lag.max = 30)

result_silent <- residualsHMM(x = sim_data$x, HMM = HMM_fit, obsdist = "pois",
                              lag.max = 30, verbose = FALSE)

Ordinary Residuals for Hidden Semi-Markov Models

Description

Computes ordinary (probability integral transform) residuals for a fitted Hidden Semi-Markov Model (HSMM). Generates quantile-quantile and autocorrelation plots for assessing model fit via simulation-based envelopes.

Usage

residualsHSMM(x, HSMM, obsdist, dwelldist, M = NA, lag.max = 50, nsim = 100,
              use.theoretical.acf = FALSE, verbose = TRUE, seed = NULL)

Arguments

Details

The function computes probability integral transform (PIT) residuals for HSMM observations by combining forward-backward probabilities with cumulative distribution functions of the state-dependent observations and dwell-time distributions. Residuals are assessed via:

• QQ plots comparing observed residuals to the expected standard normal distribution with a 95% simulation-based envelope.

• Autocorrelation function plots to detect remaining temporal dependence with simulation-based or theoretical 95% bands.

This provides a visual check of model adequacy, where ideally residuals behave as independent standard normal variables.

Value

The function generates QQ and ACF plots and prints messages about residuals outside the 95% envelope. An invisible list is also returned containing diagnostic information.

Author(s)

Aimee Cody

See Also

findmleHSMM, for estimating HSMM parameters from data. generateHSMM, for simulating HSMM data. residualsHMM, for residual analysis of hidden Markov models.

Examples

J <- 3
Pi <- matrix(c(0.0, 0.6, 0.4,
               0.5, 0.0, 0.5,
               0.3, 0.7, 0.0), nrow = J, byrow = TRUE)
obspar <- list(
  mean = c(-2, 0, 3),
  sd   = c(1, 1.5, 2)
)
dwellpar <- list(
  lambda = c(3, 5, 4)
)
# Simulate HSMM data
sim <- generateHSMM(n = 200, J = J, obsdist = "norm",
                    dwelldist = "pois", obspar = obspar,
                    dwellpar = dwellpar, Pi = Pi)
# Fit HSMM using the true parameters
HSMM_true <- findmleHSMM(x = sim$x, J = J, M = 100,
                          obsdist = "norm", dwelldist = "pois",
                          obspar = obspar, dwellpar = dwellpar,
                          Pi = Pi)
# Compute residuals and diagnostic plots
residualsHSMM(x = sim$x, HSMM = HSMM_true,
              obsdist = "norm", dwelldist = "pois",
              M = 100, nsim = 50, lag.max = 30)

Weekly maxima of geomagnetic rate-of-change data from Eskdalemuir Magnetic Observatory

Description

A dataset containing weekly block maxima of the horizontal geomagnetic field rate-of-change at the Eskdalemuir Magnetic Observatory, Scotland, United Kingdom. These data quantify the variability of the Earth's magnetic field and capture periods of extreme geomagnetic disturbance from 1999 to 2022.

Usage

data(weekly_maxima)

Format

A numeric vector of weekly maxima of

R_1(t) = \sqrt{(X_t - X_{t-1})^2 + (Y_t - Y_{t-1})^2},

where X_t and Y_t are the horizontal north and east components of the Earth's magnetic field, respectively, measured at one-minute intervals.

Details

Aggregating to a weekly scale smooths out short-lived fluctuations while retaining the most significant geomagnetic disturbances, as well as reducing short-term serial dependencies. The extreme nature of the data motivates modelling via the Generalised Extreme Value (GEV) distribution, with density

f(x) = \frac{1}{\sigma} \left[1+\xi\left(\frac{x-\mu}{\sigma}\right)\right]^{-\frac{1}{\xi}-1} \exp\left\{-\left[1+\xi\left(\frac{x-\mu}{\sigma}\right)\right]^{-\frac{1}{\xi}}\right\},

where \mu, \sigma, and \xi are the location, scale, and shape parameters, respectively.

Source

Eskdalemuir Magnetic Observatory, British Geological Survey, United Kingdom.

Derived and aggregated for use in HMMHSMM.

Examples

data(weekly_maxima)
hist(weekly_maxima, breaks = 30, main = "Weekly maxima of geomagnetic variability",
     xlab = "Weekly maximum rate of change (nT/min)")