Type:	Package
Title:	Fast Numerical Maximum Likelihood Estimation for Latent Markov Models
Version:	2.0.5
Description:	A variety of latent Markov models, including hidden Markov models, hidden semi-Markov models, state-space models and continuous-time variants can be formulated and estimated within the same framework via directly maximising the likelihood function using the so-called forward algorithm. Applied researchers often need custom models that standard software does not easily support. Writing tailored 'R' code offers flexibility but suffers from slow estimation speeds. We address these issues by providing easy-to-use functions (written in 'C++' for speed) for common tasks like the forward algorithm. These functions can be combined into custom models in a Lego-type approach, offering up to 10-20 times faster estimation via standard numerical optimisers. To aid in building fully custom likelihood functions, several vignettes are included that show how to simulate data from and estimate all the above model classes.
URL:	https://janoleko.github.io/LaMa/
License:	GPL-3
Encoding:	UTF-8
Imports:	Rcpp, mgcv, Matrix, stats, utils, MASS, splines2, methods, circular, sn, numDeriv
LinkingTo:	Rcpp, RcppArmadillo
Depends:	R (≥ 3.5.0), RTMB
RoxygenNote:	7.3.2
Suggests:	knitr, rmarkdown, testthat (≥ 3.0.0), PHSMM, MSwM, scales
VignetteBuilder:	knitr
Config/testthat/edition:	3
LazyData:	true
NeedsCompilation:	yes
Packaged:	2025-06-13 09:16:07 UTC; jan-ole
Author:	Jan-Ole Koslik [aut, cre]
Maintainer:	Jan-Ole Koslik <jan-ole.koslik@uni-bielefeld.de>
Repository:	CRAN
Date/Publication:	2025-06-15 09:10:17 UTC

LaMa: Fast Numerical Maximum Likelihood Estimation for Latent Markov Models

Description

A variety of latent Markov models, including hidden Markov models, hidden semi-Markov models, state-space models and continuous-time variants can be formulated and estimated within the same framework via directly maximising the likelihood function using the so-called forward algorithm. Applied researchers often need custom models that standard software does not easily support. Writing tailored 'R' code offers flexibility but suffers from slow estimation speeds. We address these issues by providing easy-to-use functions (written in 'C++' for speed) for common tasks like the forward algorithm. These functions can be combined into custom models in a Lego-type approach, offering up to 10-20 times faster estimation via standard numerical optimisers. To aid in building fully custom likelihood functions, several vignettes are included that show how to simulate data from and estimate all the above model classes.

Author(s)

Maintainer: Jan-Ole Koslik jan-ole.koslik@uni-bielefeld.de (ORCID)

See Also

Useful links:

• https://janoleko.github.io/LaMa/

Calculate the index of the first observation of each track based on an ID variable

Description

Function to conveniently calculate the trackInd variable that is needed internally when fitting a model to longitudinal data with multiple tracks.

Usage

calc_trackInd(ID)

Arguments

Value

A vector of indices of the first observation of each track which can be passed to the forward and forward_g to sum likelihood contributions of each track

Examples

uniqueID = c("Animal1", "Animal2", "Animal3")
ID = rep(uniqueID, c(100, 200, 300))
trackInd = calc_trackInd(ID)

Evaluate trigonometric basis expansion

Description

This function can be used to evaluate a trigonometric basis expansion for a given periodic variable and period. It can also be used in formulas passed to make_matrices.

Usage

cosinor(x = 1:24, period = 24, eval = TRUE)

Arguments

Details

The returned basis can be used for linear predictors of the form

\eta^{(t)} = \beta_0 + \sum_{k} \bigl( \beta_{1k} \sin(\frac{2 \pi t}{period_k}) + \beta_{2k} \cos(\frac{2 \pi t}{period_k}) \bigr).

This is relevant for modeling e.g. diurnal variation and the flexibility can be increased by adding smaller frequencies (i.e. increasing the length of period).

Value

either a desing matrix with the evaluated cosinor terms (eval = TRUE) or a character vector with the terms as strings (eval = FALSE).

Examples

## Evaluate cosinor terms
# builds design matrix
X = cosinor(1:24, period = 24)
X = cosinor(1:24, period = c(24, 12, 6))

## Usage in model formulas
# e.g. frequencies of 24 and 12 hours + interaction with temperature
form = ~ x + temp * cosinor(hour, c(24, 12)) 
data = data.frame(x = runif(24), temp = rnorm(24,20), hour = 1:24)
modmat = make_matrices(form, data = data)

State dwell-time distributions of periodically inhomogeneous Markov chains

Description

Computes the dwell-time distribution of a periodically inhomogeneous Markov chain for a given transition probability matrix.

Usage

ddwell(x, Gamma, time = NULL, state = NULL)

Arguments

Details

For Markov chains whose transition probabilities vary only periodically, which is achieved for example by expressing the transition probability matrix as a periodic function of the time of day using tpm_p or cosinor, the probability distribution of time spent in a state can be computed analytically. This function computes said distribution, either for a specific time point (conditioning on transitioning into the state at that time point) or for the overall distribution (conditioning on transitioning into the state at any time point).

Value

either time-varying dwell-time distribution(s) if time is specified, or overall dwell-time distribution if time is NULL. If more than one state is specified, a named list over states is returned.

References

Koslik, J. O., Feldmann, C. C., Mews, S., Michels, R., & Langrock, R. (2023). Inference on the state process of periodically inhomogeneous hidden Markov models for animal behavior. arXiv preprint arXiv:2312.14583.

Examples

# setting parameters for trigonometric link
beta = matrix(c(-1, 2, -1, -2, 1, -1), nrow = 2, byrow = TRUE)
Gamma = tpm_p(beta = beta, degree = 1)

# at specific times and for specific state
ddwell(1:20, Gamma, time = 1:4, state = 1)
# results in 4x20 matrix

# or overall distribution for all states
ddwell(1:20, Gamma)
# results in list of length 2, each element is a vector of length 20

Reparametrised multivariate Gaussian distribution

Description

Density function of the multivariate Gaussian distribution reparametrised in terms of its precision matrix (inverse variance). This implementation is particularly useful for defining the joint log-likelihood with penalised splines or i.i.d. random effects that have a multivariate Gaussian distribution with fixed precision/ penalty matrix \lambda S. As S is fixed and only scaled by \lambda, it is more efficient to precompute the determinant of S (for the normalisation constant) and only scale the quadratic form by \lambda when multiple spline parameters/ random effects with different \lambda's but the same penalty matrix S are evaluated.

Usage

dgmrf2(x, mu = 0, S, lambda, logdetS = NULL, log = FALSE)

Arguments

Details

This implementation allows for automatic differentiation with RTMB.

Value

vector of density values

Examples

x = matrix(runif(30), nrow = 3)

# iid random effects
S = diag(10)
sigma = c(1, 2, 3) # random effect standard deviations
lambda = 1 / sigma^2
d = dgmrf2(x, 0, S, lambda)

# P-splines
L = diff(diag(10), diff = 2) # second-order difference matrix
S = t(L) %*% L
lambda = c(1,2,3)
d = dgmrf2(x, 0, S, lambda, log = TRUE)

Dirichlet distribution

Description

Density of the Dirichlet distribution.

Usage

ddirichlet(x, alpha, log = TRUE)

Arguments

Details

This implementation of ddirichlet allows for automatic differentiation with RTMB.

Value

ddirichlet gives the density.

Examples

ddirichlet(c(0.2, 0.3, 0.5), c(1, 2, 3))

Forward algorithm with homogeneous transition probability matrix

Description

Calculates the log-likelihood of a sequence of observations under a homogeneous hidden Markov model using the forward algorithm.

Usage

forward(
  delta,
  Gamma,
  allprobs,
  trackID = NULL,
  ad = NULL,
  report = TRUE,
  logspace = FALSE
)

Arguments

Value

log-likelihood for given data and parameters

See Also

Other forward algorithms: forward_g(), forward_hsmm(), forward_ihsmm(), forward_p(), forward_phsmm()

Examples

## negative log likelihood function
nll = function(par, step) {
 # parameter transformations for unconstrained optimisation
 Gamma = tpm(par[1:2]) # multinomial logit link
 delta = stationary(Gamma) # stationary HMM
 mu = exp(par[3:4])
 sigma = exp(par[5:6])
 # calculate all state-dependent probabilities
 allprobs = matrix(1, length(step), 2)
 ind = which(!is.na(step))
 for(j in 1:2) allprobs[ind,j] = dgamma2(step[ind], mu[j], sigma[j])
 # simple forward algorithm to calculate log-likelihood
 -forward(delta, Gamma, allprobs)
}

## fitting an HMM to the trex data
par = c(-2,-2,            # initial tpm params (logit-scale)
        log(c(0.3, 2.5)), # initial means for step length (log-transformed)
        log(c(0.2, 1.5))) # initial sds for step length (log-transformed)
mod = nlm(nll, par, step = trex$step[1:1000])

General forward algorithm with time-varying transition probability matrix

Description

Calculates the log-likelihood of a sequence of observations under a hidden Markov model with time-varying transition probabilities using the forward algorithm.

Usage

forward_g(
  delta,
  Gamma,
  allprobs,
  trackID = NULL,
  ad = NULL,
  report = TRUE,
  logspace = FALSE
)

Arguments

Value

log-likelihood for given data and parameters

See Also

Other forward algorithms: forward(), forward_hsmm(), forward_ihsmm(), forward_p(), forward_phsmm()

Examples

## Simple usage
Gamma = array(c(0.9, 0.2, 0.1, 0.8), dim = c(2,2,10))
delta = c(0.5, 0.5)
allprobs = matrix(0.5, 10, 2)
forward_g(delta, Gamma, allprobs)


## Full model fitting example
## negative log likelihood function
nll = function(par, step, Z) {
 # parameter transformations for unconstrained optimisation
 beta = matrix(par[1:6], nrow = 2)
 Gamma = tpm_g(Z, beta) # multinomial logit link for each time point
 delta = stationary(Gamma[,,1]) # stationary HMM
 mu = exp(par[7:8])
 sigma = exp(par[9:10])
 # calculate all state-dependent probabilities
 allprobs = matrix(1, length(step), 2)
 ind = which(!is.na(step))
 for(j in 1:2) allprobs[ind,j] = dgamma2(step[ind], mu[j], sigma[j])
 # simple forward algorithm to calculate log-likelihood
 -forward_g(delta, Gamma, allprobs)
}

## fitting an HMM to the trex data
par = c(-1.5,-1.5,        # initial tpm intercepts (logit-scale)
        rep(0, 4),        # initial tpm slopes
        log(c(0.3, 2.5)), # initial means for step length (log-transformed)
        log(c(0.2, 1.5))) # initial sds for step length (log-transformed)
mod = nlm(nll, par, step = trex$step[1:500], Z = trigBasisExp(trex$tod[1:500]))

Forward algorithm for homogeneous hidden semi-Markov models

Description

Calculates the (approximate) log-likelihood of a sequence of observations under a homogeneous hidden semi-Markov model using a modified forward algorithm.

Usage

forward_hsmm(
  dm,
  omega,
  allprobs,
  trackID = NULL,
  delta = NULL,
  eps = 1e-10,
  report = TRUE
)

Arguments

Details

Hidden semi-Markov models (HSMMs) are a flexible extension of HMMs, where the state duration distribution is explicitly modelled by a distribution on the positive integers. For direct numerical maximum likelhood estimation, HSMMs can be represented as HMMs on an enlarged state space (of size M) and with structured transition probabilities.

This function is designed to be used with automatic differentiation based on the R package RTMB. It will be very slow without it!

Value

log-likelihood for given data and parameters

References

Langrock, R., & Zucchini, W. (2011). Hidden Markov models with arbitrary state dwell-time distributions. Computational Statistics & Data Analysis, 55(1), 715-724.

Koslik, J. O. (2025). Hidden semi-Markov models with inhomogeneous state dwell-time distributions. Computational Statistics & Data Analysis, 209, 108171.

See Also

Other forward algorithms: forward(), forward_g(), forward_ihsmm(), forward_p(), forward_phsmm()

Examples

# currently no examples

Forward algorithm for hidden semi-Markov models with inhomogeneous state durations and/ or conditional transition probabilities

Description

Calculates the (approximate) log-likelihood of a sequence of observations under an inhomogeneous hidden semi-Markov model using a modified forward algorithm.

Usage

forward_ihsmm(
  dm,
  omega,
  allprobs,
  trackID = NULL,
  delta = NULL,
  startInd = NULL,
  eps = 1e-10,
  report = TRUE
)

Arguments

Details

Hidden semi-Markov models (HSMMs) are a flexible extension of HMMs, where the state duration distribution is explicitly modelled by a distribution on the positive integers. This function can be used to fit HSMMs where the state-duration distribution and/ or the conditional transition probabilities vary with covariates. For direct numerical maximum likelhood estimation, HSMMs can be represented as HMMs on an enlarged state space (of size M) and with structured transition probabilities.

This function is designed to be used with automatic differentiation based on the R package RTMB. It will be very slow without it!

Value

log-likelihood for given data and parameters

References

Koslik, J. O. (2025). Hidden semi-Markov models with inhomogeneous state dwell-time distributions. Computational Statistics & Data Analysis, 209, 108171.

See Also

Other forward algorithms: forward(), forward_g(), forward_hsmm(), forward_p(), forward_phsmm()

Examples

# currently no examples

Forward algorithm with for periodically varying transition probability matrices

Description

Calculates the log-likelihood of a sequence of observations under a hidden Markov model with periodically varying transition probabilities using the forward algorithm.

Usage

forward_p(
  delta,
  Gamma,
  allprobs,
  tod,
  trackID = NULL,
  ad = NULL,
  report = TRUE,
  logspace = FALSE
)

Arguments

Details

When the transition probability matrix only varies periodically (e.g. as a function of time of day), there are only L unique matrices if L is the period length (e.g. L=24 for hourly data and time-of-day variation). Thus, it is much more efficient to only calculate these L matrices and index them by a time variable (e.g. time of day or day of year) instead of calculating such a matrix for each index in the data set (which would be redundant). This function allows for that by only expecting a transition probability matrix for each time point in a period and an integer valued (1, \dots, L) time variable that maps the data index to the according time.

Value

log-likelihood for given data and parameters

See Also

Other forward algorithms: forward(), forward_g(), forward_hsmm(), forward_ihsmm(), forward_phsmm()

Examples

## negative log likelihood function
nll = function(par, step, tod) {
 # parameter transformations for unconstrained optimisation
 beta = matrix(par[1:6], nrow = 2)
 Gamma = tpm_p(1:24, beta = beta) # multinomial logit link for each time point
 delta = stationary_p(Gamma, tod[1]) # stationary HMM
 mu = exp(par[7:8])
 sigma = exp(par[9:10])
 # calculate all state-dependent probabilities
 allprobs = matrix(1, length(step), 2)
 ind = which(!is.na(step))
 for(j in 1:2) allprobs[ind,j] = dgamma2(step[ind], mu[j], sigma[j])
 # simple forward algorithm to calculate log-likelihood
 -forward_p(delta, Gamma, allprobs, tod)
}

## fitting an HMM to the nessi data
par = c(-2,-2,            # initial tpm intercepts (logit-scale)
        rep(0, 4),        # initial tpm slopes
        log(c(0.3, 2.5)), # initial means for step length (log-transformed)
        log(c(0.2, 1.5))) # initial sds for step length (log-transformed)
mod = nlm(nll, par, step = trex$step[1:500], tod = trex$tod[1:500])

Forward algorithm for hidden semi-Markov models with periodically inhomogeneous state durations and/ or conditional transition probabilities

Description

Hidden semi-Markov models (HSMMs) are a flexible extension of HMMs, where the state duration distribution is explicitly modelled by a distribution on the positive integers. This function can be used to fit HSMMs where the state-duration distribution and/ or the conditional transition probabilities vary with covariates. For direct numerical maximum likelhood estimation, HSMMs can be represented as HMMs on an enlarged state space (of size M) and with structured transition probabilities.

This function can be used to fit HSMMs where the state-duration distribution and/ or the conditional transition probabilities vary periodically. In the special case of periodic variation (as compared to arbitrary covariate influence), this version is to be preferred over forward_ihsmm because it computes the correct periodically stationary distribution and no observations are lost for the approximation.

This function is designed to be used with automatic differentiation based on the R package RTMB. It will be very slow without it!

Usage

forward_phsmm(
  dm,
  omega,
  allprobs,
  tod,
  trackID = NULL,
  delta = NULL,
  eps = 1e-10,
  report = TRUE
)

Arguments

Details

Calculates the (approximate) log-likelihood of a sequence of observations under a periodically inhomogeneous hidden semi-Markov model using a modified forward algorithm.

Value

log-likelihood for given data and parameters

References

Koslik, J. O. (2025). Hidden semi-Markov models with inhomogeneous state dwell-time distributions. Computational Statistics & Data Analysis, 209, 108171.

See Also

Other forward algorithms: forward(), forward_g(), forward_hsmm(), forward_ihsmm(), forward_p()

Examples

# currently no examples

Forward algorithm for hidden semi-Markov models with homogeneous transition probability matrix

Description

Hidden semi-Markov models (HSMMs) are a flexible extension of HMMs that can be approximated by HMMs on an enlarged state space (of size M) and with structured transition probabilities.

Usage

forward_s(delta, Gamma, allprobs, sizes)

Arguments

Value

log-likelihood for given data and parameters

Examples

## generating data from homogeneous 2-state HSMM
mu = c(0, 6)
lambda = c(6, 12)
omega = matrix(c(0,1,1,0), nrow = 2, byrow = TRUE)
# simulation
# for a 2-state HSMM the embedded chain always alternates between 1 and 2
s = rep(1:2, 100)
C = x = numeric(0)
for(t in 1:100){
  dt = rpois(1, lambda[s[t]])+1 # shifted Poisson
  C = c(C, rep(s[t], dt))
  x = c(x, rnorm(dt, mu[s[t]], 1.5)) # fixed sd 2 for both states
}

## negative log likelihood function
mllk = function(theta.star, x, sizes){
  # parameter transformations for unconstraint optimization
  omega = matrix(c(0,1,1,0), nrow = 2, byrow = TRUE) # omega fixed (2-states)
  lambda = exp(theta.star[1:2]) # dwell time means
  dm = list(dpois(1:sizes[1]-1, lambda[1]), dpois(1:sizes[2]-1, lambda[2]))
  Gamma = tpm_hsmm2(omega, dm)
  delta = stationary(Gamma) # stationary
  mu = theta.star[3:4]
  sigma = exp(theta.star[5:6])
  # calculate all state-dependent probabilities
  allprobs = matrix(1, length(x), 2)
  for(j in 1:2){ allprobs[,j] = dnorm(x, mu[j], sigma[j]) }
  # return negative for minimization
  -forward_s(delta, Gamma, allprobs, sizes)
}

## fitting an HSMM to the data
theta.star = c(log(5), log(10), 1, 4, log(2), log(2))
mod = nlm(mllk, theta.star, x = x, sizes = c(20, 30), stepmax = 5)

Forward algorithm for hidden semi-Markov models with periodically varying transition probability matrices

Description

Hidden semi-Markov models (HSMMs) are a flexible extension of HMMs that can be approximated by HMMs on an enlarged state space (of size M) and with structured transition probabilities. Recently, this inference procedure has been generalised to allow either the dwell-time distributions or the conditional transition probabilities to depend on external covariates such as the time of day. This special case is implemented here. This function allows for that, by expecting a transition probability matrix for each time point in a period, and an integer valued (1, \dots, L) time variable that maps the data index to the according time.

Usage

forward_sp(delta, Gamma, allprobs, sizes, tod)

Arguments

Value

log-likelihood for given data and parameters

Examples

## generating data from homogeneous 2-state HSMM
mu = c(0, 6)
beta = matrix(c(log(4),log(6),-0.2,0.2,-0.1,0.4), nrow=2)
# time varying mean dwell time
Lambda = exp(cbind(1, trigBasisExp(1:24, 24, 1))%*%t(beta))
omega = matrix(c(0,1,1,0), nrow = 2, byrow = TRUE)
# simulation
# for a 2-state HSMM the embedded chain always alternates between 1 and 2
s = rep(1:2, 100)
C = x = numeric(0)
tod = rep(1:24, 50) # time of day variable
time = 1
for(t in 1:100){
  dt = rpois(1, Lambda[tod[time], s[t]])+1 # dwell time depending on time of day
  time = time + dt
  C = c(C, rep(s[t], dt))
  x = c(x, rnorm(dt, mu[s[t]], 1.5)) # fixed sd 2 for both states
}
tod = tod[1:length(x)]

## negative log likelihood function
mllk = function(theta.star, x, sizes, tod){
  # parameter transformations for unconstraint optimization
  omega = matrix(c(0,1,1,0), nrow = 2, byrow = TRUE) # omega fixed (2-states)
  mu = theta.star[1:2]
  sigma = exp(theta.star[3:4])
  beta = matrix(theta.star[5:10], nrow=2)
  # time varying mean dwell time
  Lambda = exp(cbind(1, trigBasisExp(1:24, 24, 1))%*%t(beta))
  dm = list()
  for(j in 1:2){
    dm[[j]] = sapply(1:sizes[j]-1, dpois, lambda = Lambda[,j])
  }
  Gamma = tpm_phsmm2(omega, dm)
  delta = stationary_p(Gamma, tod[1])
  # calculate all state-dependent probabilities
  allprobs = matrix(1, length(x), 2)
  for(j in 1:2){ allprobs[,j] = dnorm(x, mu[j], sigma[j]) }
  # return negative for minimization
  -forward_sp(delta, Gamma, allprobs, sizes, tod)
}

## fitting an HSMM to the data
theta.star = c(1, 4, log(2), log(2), # state-dependent parameters
                 log(4), log(6), rep(0,4)) # state process parameters dm
# mod = nlm(mllk, theta.star, x = x, sizes = c(10, 15), tod = tod, stepmax = 5)

Reparametrised gamma distribution

Description

Density, distribution function, quantile function and random generation for the gamma distribution reparametrised in terms of mean and standard deviation.

Usage

dgamma2(x, mean = 1, sd = 1, log = FALSE)

pgamma2(q, mean = 1, sd = 1, lower.tail = TRUE, log.p = FALSE)

qgamma2(p, mean = 1, sd = 1, lower.tail = TRUE, log.p = FALSE)

rgamma2(n, mean = 1, sd = 1)

Arguments

Details

This implementation allows for automatic differentiation with RTMB.

Value

dgamma2 gives the density, pgamma2 gives the distribution function, qgamma2 gives the quantile function, and rgamma2 generates random deviates.

Examples

x = rgamma2(1)
d = dgamma2(x)
p = pgamma2(x)
q = qgamma2(p)

Computes generalised determinant

Description

Computes generalised determinant

Usage

gdeterminant(x, eps = NULL, log = TRUE)

Arguments

Value

generalised log-determinant of x

Build the generator matrix of a continuous-time Markov chain

Description

This function builds the infinitesimal generator matrix for a continuous-time Markov chain from an unconstrained parameter vector.

Usage

generator(param, byrow = FALSE, report = TRUE)

Arguments

Value

infinitesimal generator matrix of dimension c(N,N)

See Also

Other transition probability matrix functions: tpm(), tpm_cont(), tpm_emb(), tpm_emb_g(), tpm_g(), tpm_p()

Examples

# 2 states: 2 free off-diagonal elements
generator(rep(-1, 2))
# 3 states: 6 free off-diagonal elements
generator(rep(-2, 6))

Extract log-likelihood from qremlModel object

Description

Extract log-likelihood from qremlModel object

Usage

## S3 method for class 'qremlModel'
logLik(object, ...)

Arguments

Value

An object of class "logLik"

Build the design and the penalty matrix for models involving penalised splines based on a formula and a data set

Description

Build the design and the penalty matrix for models involving penalised splines based on a formula and a data set

Usage

make_matrices(formula, data, knots = NULL)

Arguments

Value

a list containing the design matrix Z, a (potentially nested) list of penalty matrices S, the formula, the data, the knots, and the original mod object returned by mgcv. Note that for tensorproduct smooths, the corresponding list entry is itself a list, containing the d marginal penalty matrices if d is the dimension of the tensor product.

Examples

data = data.frame(x = runif(100), 
                  y = runif(100),
                  g = factor(rep(1:10, each = 10)))

# unvariate thin plate regression spline
modmat = make_matrices(~ s(x), data)
# univariate P-spline
modmat = make_matrices(~ s(x, bs = "ps"), data)
# adding random intercept
modmat = make_matrices(~ s(g, bs = "re") + s(x, bs = "ps"), data)
# tensorproduct of x and y
modmat = make_matrices(~ s(x) + s(y) + ti(x,y), data)

Build a standardised P-Spline design matrix and the associated P-Spline penalty matrix

Description

This function builds the B-spline design matrix for a given data vector. Importantly, the B-spline basis functions are normalised such that the integral of each basis function is 1, hence this basis can be used for spline-based density estimation, when the basis functions are weighted by non-negative weights summing to one.

Usage

make_matrices_dens(
  x,
  k,
  type = "real",
  degree = 3,
  knots = NULL,
  diff_order = 2,
  pow = 0.5,
  npoints = 10000
)

Arguments

Value

list containing the design matrix Z, the penalty matrix S, the prediction design matrix Z_predict, the prediction grid xseq, and details for the basis expansion.

Examples

set.seed(1)
# real-valued
x <- rnorm(100)
modmat <- make_matrices_dens(x, k = 20)
# positive-continuouos
x <- rgamma2(100, mean = 5, sd = 2)
modmat <- make_matrices_dens(x, k = 20, type = "positive")
# circular
x <- rvm(100, mu = 0, kappa = 2)
modmat <- make_matrices_dens(x, k = 20, type = "circular")
# bounded in an interval
x <- rbeta(100, 1, 2)
modmat <- make_matrices_dens(x, k = 20)

AD-compatible minimum and maximum functions

Description

These functions compute the parallel minimum/ maximum of two vector-valued inputs and are compatible with automatic differentiation using RTMB.

Usage

min2(x, y)

max2(x, y)

Arguments

Value

min2 returns the parallel minimum and max2 the parallel maximum of x and y

Examples

x <- c(1, 4, 8, 2)
y <- c(2, 5, 3, 7)
min2(x, y)
max2(x, y)

Loch Ness Monster Acceleration Data

Description

A small group of researchers managed to put an accelerometer on the Loch Ness Monster and collected data for a few days. Now we have a data set of the overall dynamic body acceleration (ODBA) of the creature.

Usage

nessi

Format

A data frame with 5.000 rows and 3 variables:

ODBA

overall dynamci body acceleration

logODBA

logarithm of overall dynamic body acceleration

state

hidden state variable

Source

Generated for example purposes.

Computes penalty based on quadratic form

Description

This function computes quadratic penalties of the form

0.5 \sum_{i} \lambda_i b_i^T S_i b_i,

with smoothing parameters \lambda_i, coefficient vectors b_i, and fixed penalty matrices S_i.

It is intended to be used inside the penalised negative log-likelihood function when fitting models with penalised splines or simple random effects via quasi restricted maximum likelihood (qREML) with the qreml function. For qreml to work, the likelihood function needs to be compatible with the RTMB R package to enable automatic differentiation.

Usage

penalty(re_coef, S, lambda)

Arguments

Details

Caution: The formatting of re_coef needs to match the structure of the parameter list in your penalised negative log-likelihood function, i.e. you cannot have two random effect vectors of different names (different list elements in the parameter list), combine them into a matrix inside your likelihood and pass the matrix to penalty. If these are seperate random effects, each with its own name, they need to be passed as a list to penalty. Moreover, the ordering of re_coef needs to match the character vector random specified in qreml.

Value

returns the penalty value and reports to qreml.

References

Koslik, J. O. (2024). Efficient smoothness selection for nonparametric Markov-switching models via quasi restricted maximum likelihood. arXiv preprint arXiv:2411.11498.

See Also

qreml for the qREML algorithm

Examples

# Example with a single random effect
re = rep(0, 5)
S = diag(5)
lambda = 1
penalty(re, S, lambda)

# Example with two random effects, 
# where one element contains two random effects of similar structure
re = list(matrix(0, 2, 5), rep(0, 4))
S = list(diag(5), diag(4))
lambda = c(1,1,2) # length = total number of random effects
penalty(re, S, lambda)

# Full model-fitting example
data = trex[1:1000,] # subset

# initial parameter list
par = list(logmu = log(c(0.3, 1)), # step mean
           logsigma = log(c(0.2, 0.7)), # step sd
           beta0 = c(-2,-2), # state process intercept
           betaspline = matrix(rep(0, 18), nrow = 2)) # state process spline coefs
          
# data object with initial penalty strength lambda
dat = list(step = data$step, # step length
           tod = data$tod, # time of day covariate
           N = 2, # number of states
           lambda = rep(10,2)) # initial penalty strength

# building model matrices
modmat = make_matrices(~ s(tod, bs = "cp"), 
                       data = data.frame(tod = 1:24), 
                       knots = list(tod = c(0,24))) # wrapping points
dat$Z = modmat$Z # spline design matrix
dat$S = modmat$S # penalty matrix

# penalised negative log-likelihood function
pnll = function(par) {
  getAll(par, dat) # makes everything contained available without $
  Gamma = tpm_g(Z, cbind(beta0, betaspline)) # transition probabilities
  delta = stationary_p(Gamma, t = 1) # initial distribution
  mu = exp(logmu) # step mean
  sigma = exp(logsigma) # step sd
  # calculating all state-dependent densities
  allprobs = matrix(1, nrow = length(step), ncol = N)
  ind = which(!is.na(step)) # only for non-NA obs.
  for(j in 1:N) allprobs[ind,j] = dgamma2(step[ind],mu[j],sigma[j])
  -forward_g(delta, Gamma[,,tod], allprobs) +
      penalty(betaspline, S, lambda) # this does all the penalization work
}

# model fitting
mod = qreml(pnll, par, dat, random = "betaspline")

Computes generalised quadratic-form penalties

Description

This function computes a quadratic penalty of the form

0.5 \sum_{i} \lambda_i b^T S_i b,

with smoothing parameters \lambda_i, coefficient vector b, and fixed penalty matrices S_i. This generalises the penalty by allowing subsets of the coefficient vector b to be penalised multiple times with different smoothing parameters, which is necessary for tensor products, functional random effects or adaptive smoothing.

It is intended to be used inside the penalised negative log-likelihood function when fitting models with penalised splines or simple random effects via quasi restricted maximum likelihood (qREML) with the qreml function. For qreml to work, the likelihood function needs to be compatible with the RTMB R package to enable automatic differentiation.

Usage

penalty2(re_coef, S, lambda)

Arguments

Details

Caution: The formatting of re_coef needs to match the structure of the parameter list in your penalised negative log-likelihood function, i.e. you cannot have two random effect vectors of different names (different list elements in the parameter list), combine them into a matrix inside your likelihood and pass the matrix to penalty. If these are seperate random effects, each with its own name, they need to be passed as a list to penalty. Moreover, the ordering of re_coef needs to match the character vector random specified in qreml.

Value

returns the penalty value and reports to qreml.

See Also

qreml for the qREML algorithm

Examples

# Example with a single random effect
re = rep(0, 5)
S = diag(5)
lambda = 1
penalty(re, S, lambda)

# Example with two random effects, 
# where one element contains two random effects of similar structure
re = list(matrix(0, 2, 5), rep(0, 4))
S = list(diag(5), diag(4))
lambda = c(1,1,2) # length = total number of random effects
penalty(re, S, lambda)

# Full model-fitting example
data = trex[1:1000,] # subset

# initial parameter list
par = list(logmu = log(c(0.3, 1)), # step mean
           logsigma = log(c(0.2, 0.7)), # step sd
           beta0 = c(-2,-2), # state process intercept
           betaspline = matrix(rep(0, 18), nrow = 2)) # state process spline coefs
          
# data object with initial penalty strength lambda
dat = list(step = data$step, # step length
           tod = data$tod, # time of day covariate
           N = 2, # number of states
           lambda = rep(10,2)) # initial penalty strength

# building model matrices
modmat = make_matrices(~ s(tod, bs = "cp"), 
                       data = data.frame(tod = 1:24), 
                       knots = list(tod = c(0,24))) # wrapping points
dat$Z = modmat$Z # spline design matrix
dat$S = modmat$S # penalty matrix

# penalised negative log-likelihood function
pnll = function(par) {
  getAll(par, dat) # makes everything contained available without $
  Gamma = tpm_g(Z, cbind(beta0, betaspline)) # transition probabilities
  delta = stationary_p(Gamma, t = 1) # initial distribution
  mu = exp(logmu) # step mean
  sigma = exp(logsigma) # step sd
  # calculating all state-dependent densities
  allprobs = matrix(1, nrow = length(step), ncol = N)
  ind = which(!is.na(step)) # only for non-NA obs.
  for(j in 1:N) allprobs[ind,j] = dgamma2(step[ind],mu[j],sigma[j])
  -forward_g(delta, Gamma[,,tod], allprobs) +
      penalty(betaspline, S, lambda) # this does all the penalization work
}

# model fitting
mod = qreml(pnll, par, dat, random = "betaspline")

Plot pseudo-residuals

Description

Plot pseudo-residuals computed by pseudo_res.

Usage

## S3 method for class 'LaMaResiduals'
plot(x, hist = FALSE, col = "darkblue", lwd = 1.5, main = NULL, ...)

Arguments

Value

NULL, plots the pseudo-residuals in a 2- or 3-panel layout

Examples

## pseudo-residuals for the trex data
step = trex$step[1:200]

nll = function(par){
  getAll(par)
  Gamma = tpm(logitGamma)
  delta = stationary(Gamma)
  mu = exp(logMu); REPORT(mu)
  sigma = exp(logSigma); REPORT(sigma)
  allprobs = matrix(1, length(step), 2)
  ind = which(!is.na(step))
  for(j in 1:2) allprobs[ind,j] = dgamma2(step[ind], mu[j], sigma[j])
  -forward(delta, Gamma, allprobs)
}

par = list(logitGamma = c(-2,-2), 
           logMu = log(c(0.3, 2.5)), 
           logSigma = log(c(0.3, 0.5)))
           
obj = MakeADFun(nll, par)
opt = nlminb(obj$par, obj$fn, obj$gr)

mod = obj$report()

pres = pseudo_res(step, "gamma2", list(mean = mod$mu, sd = mod$sigma),
                  mod = mod)
                  
plot(pres)

Build the prediction design matrix based on new data and model_matrices object created by make_matrices

Description

Build the prediction design matrix based on new data and model_matrices object created by make_matrices

Usage

pred_matrix(model_matrices, newdata, exclude = NULL)

Arguments

Value

prediction design matrix for newdata with the same basis as used for model_matrices

Examples

modmat = make_matrices(~ s(x), data.frame(x = 1:10))
Z_predict = pred_matrix(modmat, data.frame(x = 1:10 - 0.5))

Build the prediction design matrix based on new data and model_matrices object created by make_matrices

Description

Build the prediction design matrix based on new data and model_matrices object created by make_matrices

Usage

## S3 method for class 'LaMa_matrices'
predict(object, newdata, ...)

Arguments

Value

prediction design matrix for newdata with the same basis as used for model_matrices

Examples

modmat = make_matrices(~ s(x), data.frame(x = 1:10))
Z_predict = predict(modmat, data.frame(x = 1:10 - 0.5))

Calculate pseudo-residuals

Description

For HMMs, pseudo-residuals are used to assess the goodness-of-fit of the model. These are based on the cumulative distribution function (CDF)

F_{X_t}(x_t) = F(x_t \mid x_1, \dots, x_{t-1}, x_{t+1}, \dots, x_T)

and can be used to quantify whether an observation is extreme relative to its model-implied distribution.

This function calculates such residuals via probability integral transform, based on the local state probabilities obtained by stateprobs or stateprobs_g and the respective parametric family.

Usage

pseudo_res(
  obs,
  dist,
  par,
  stateprobs = NULL,
  mod = NULL,
  normal = TRUE,
  discrete = NULL,
  randomise = TRUE,
  seed = NULL
)

Arguments

Details

When used for discrete pseudo-residuals, this function is just a wrapper for pseudo_res_discrete.

Value

vector of pseudo residuals

Examples

## continuous-valued observations
obs = rnorm(100)
stateprobs = matrix(0.5, nrow = 100, ncol = 2)
par = list(mean = c(1,2), sd = c(1,1))
pres = pseudo_res(obs, "norm", par, stateprobs)

## discrete-valued observations
obs = rpois(100, lambda = 1)
par = list(lambda = c(1,2))
pres = pseudo_res(obs, "pois", par, stateprobs)

## custom CDF function
obs = rnbinom(100, size = 1, prob = 0.5)
par = list(size = c(0.5, 2), prob = c(0.4, 0.6))
pres = pseudo_res(obs, pnbinom, par, stateprobs, 
                  discrete = TRUE)
# if discrete CDF function is passed, 'discrete' needs to be set to TRUE

## full example with model object
step = trex$step[1:200]

nll = function(par){
  getAll(par)
  Gamma = tpm(logitGamma)
  delta = stationary(Gamma)
  mu = exp(logMu); REPORT(mu)
  sigma = exp(logSigma); REPORT(sigma)
  allprobs = matrix(1, length(step), 2)
  ind = which(!is.na(step))
  for(j in 1:2) allprobs[ind,j] = dgamma2(step[ind], mu[j], sigma[j])
  -forward(delta, Gamma, allprobs)
}

par = list(logitGamma = c(-2,-2), 
           logMu = log(c(0.3, 2.5)), 
           logSigma = log(c(0.3, 0.5)))
           
obj = MakeADFun(nll, par)
opt = nlminb(obj$par, obj$fn, obj$gr)

mod = obj$report()

pres = pseudo_res(step, "gamma2", list(mean = mod$mu, sd = mod$sigma),
                  mod = mod)

qqnorm(pres)
abline(a = 0, b = 1)

Calculate pseudo-residuals for discrete-valued observations

Description

For HMMs, pseudo-residuals are used to assess the goodness-of-fit of the model. These are based on the cumulative distribution function (CDF)

F_{X_t}(x_t) = F(x_t \mid x_1, \dots, x_{t-1}, x_{t+1}, \dots, x_T)

and can be used to quantify whether an observation is extreme relative to its model-implied distribution.

This function calculates such residuals for discrete-valued observations, based on the local state probabilities obtained by stateprobs or stateprobs_g and the respective parametric family.

Usage

pseudo_res_discrete(
  obs,
  dist,
  par,
  stateprobs,
  normal = TRUE,
  randomise = TRUE,
  seed = NULL
)

Arguments

Details

For discrete observations, calculating pseudo residuals is slightly more involved, as the CDF is a step function. Therefore, one can calculate the lower and upper CDF values for each observation. By default, this function does exactly that and then randomly samples the interval in between to give approximately Gaussian psuedo-residuals. If randomise is set to FALSE, the lower, upper and mean pseudo-residuasl are returned.

Value

vector of pseudo residuals

Examples

obs = rpois(100, lambda = 1)
stateprobs = matrix(0.5, nrow = 100, ncol = 2)
par = list(lambda = c(1,2))
pres = pseudo_res_discrete(obs, "pois", par, stateprobs)

Quasi restricted maximum likelihood (qREML) algorithm for models with penalised splines or simple i.i.d. random effects

Description

This algorithm can be used very flexibly to fit statistical models that involve penalised splines or simple i.i.d. random effects, i.e. that have penalties of the form

0.5 \sum_{i} \lambda_i b_i^T S_i b_i,

with smoothing parameters \lambda_i, coefficient vectors b_i, and fixed penalty matrices S_i.

The qREML algorithm is typically much faster than REML or marginal ML using the full Laplace approximation method, but may be slightly less accurate regarding the estimation of the penalty strength parameters.

Under the hood, qreml uses the R package RTMB for automatic differentiation in the inner optimisation. The user has to specify the penalised negative log-likelihood function pnll structured as dictated by RTMB and use the penalty function to compute the quadratic-form penalty inside the likelihood.

Usage

qreml(
  pnll,
  par,
  dat,
  random,
  map = NULL,
  psname = "lambda",
  alpha = 0.3,
  smoothing = 1,
  maxiter = 100,
  tol = 1e-04,
  control = list(reltol = 1e-10, maxit = 1000),
  silent = 1,
  joint_unc = TRUE,
  saveall = FALSE
)

Arguments

Value

model object of class 'qremlModel'. This is a list containing:

References

Koslik, J. O. (2024). Efficient smoothness selection for nonparametric Markov-switching models via quasi restricted maximum likelihood. arXiv preprint arXiv:2411.11498.

See Also

penalty and penalty2 to compute the penalty inside the likelihood function

Examples

data = trex[1:1000,] # subset

# initial parameter list
par = list(logmu = log(c(0.3, 1)), # step mean
           logsigma = log(c(0.2, 0.7)), # step sd
           beta0 = c(-2,-2), # state process intercept
           betaspline = matrix(rep(0, 18), nrow = 2)) # state process spline coefs
          
# data object with initial penalty strength lambda
dat = list(step = data$step, # step length
           tod = data$tod, # time of day covariate
           N = 2, # number of states
           lambda = rep(10,2)) # initial penalty strength

# building model matrices
modmat = make_matrices(~ s(tod, bs = "cp"), 
                       data = data.frame(tod = 1:24), 
                       knots = list(tod = c(0,24))) # wrapping points
dat$Z = modmat$Z # spline design matrix
dat$S = modmat$S # penalty matrix

# penalised negative log-likelihood function
pnll = function(par) {
  getAll(par, dat) # makes everything contained available without $
  Gamma = tpm_g(Z, cbind(beta0, betaspline), ad = TRUE) # transition probabilities
  delta = stationary_p(Gamma, t = 1, ad = TRUE) # initial distribution
  mu = exp(logmu) # step mean
  sigma = exp(logsigma) # step sd
  # calculating all state-dependent densities
  allprobs = matrix(1, nrow = length(step), ncol = N)
  ind = which(!is.na(step)) # only for non-NA obs.
  for(j in 1:N) allprobs[ind,j] = dgamma2(step[ind],mu[j],sigma[j])
  -forward_g(delta, Gamma[,,tod], allprobs) +
      penalty(betaspline, S, lambda) # this does all the penalization work
}

# model fitting
mod = qreml(pnll, par, dat, random = "betaspline")

Quasi restricted maximum likelihood (qREML) algorithm for models with penalised splines or simple i.i.d. random effects

Description

This algorithm can be used very flexibly to fit statistical models that involve penalised splines or simple i.i.d. random effects, i.e. that have penalties of the form

0.5 \sum_{i} \lambda_i b_i^T S_i b_i,

with smoothing parameters \lambda_i, coefficient vectors b_i, and fixed penalty matrices S_i.

The qREML algorithm is typically much faster than REML or marginal ML using the full Laplace approximation method, but may be slightly less accurate regarding the estimation of the penalty strength parameters.

Under the hood, qreml uses the R package RTMB for automatic differentiation in the inner optimisation. The user has to specify the penalised negative log-likelihood function pnll structured as dictated by RTMB and use the penalty function to compute the quadratic-form penalty inside the likelihood.

Usage

qreml_old(
  pnll,
  par,
  dat,
  random,
  map = NULL,
  psname = "lambda",
  alpha = 0.25,
  smoothing = 1,
  maxiter = 100,
  tol = 1e-04,
  control = list(reltol = 1e-10, maxit = 1000),
  silent = 1,
  joint_unc = TRUE,
  saveall = FALSE
)

Arguments

Value

model object of class 'qremlModel'. This is a list containing:

References

Koslik, J. O. (2024). Efficient smoothness selection for nonparametric Markov-switching models via quasi restricted maximum likelihood. arXiv preprint arXiv:2411.11498.

See Also

penalty to compute the penalty inside the likelihood function

Examples

data = trex[1:1000,] # subset

# initial parameter list
par = list(logmu = log(c(0.3, 1)), # step mean
           logsigma = log(c(0.2, 0.7)), # step sd
           beta0 = c(-2,-2), # state process intercept
           betaspline = matrix(rep(0, 18), nrow = 2)) # state process spline coefs
          
# data object with initial penalty strength lambda
dat = list(step = data$step, # step length
           tod = data$tod, # time of day covariate
           N = 2, # number of states
           lambda = rep(100,2)) # initial penalty strength

# building model matrices
modmat = make_matrices(~ s(tod, bs = "cp"), 
                       data = data.frame(tod = 1:24), 
                       knots = list(tod = c(0,24))) # wrapping points
dat$Z = modmat$Z # spline design matrix
dat$S = modmat$S # penalty matrix

# penalised negative log-likelihood function
pnll = function(par) {
  getAll(par, dat) # makes everything contained available without $
  Gamma = tpm_g(Z, cbind(beta0, betaspline), ad = TRUE) # transition probabilities
  delta = stationary_p(Gamma, t = 1, ad = TRUE) # initial distribution
  mu = exp(logmu) # step mean
  sigma = exp(logsigma) # step sd
  # calculating all state-dependent densities
  allprobs = matrix(1, nrow = length(step), ncol = N)
  ind = which(!is.na(step)) # only for non-NA obs.
  for(j in 1:N) allprobs[ind,j] = dgamma2(step[ind],mu[j],sigma[j])
  -forward_g(delta, Gamma[,,tod], allprobs) +
      penalty(betaspline, S, lambda) # this does all the penalization work
}

# model fitting
mod = qreml_old(pnll, par, dat, random = "betaspline")

Monte Carlo version of sdreport

Description

After optimisation of an AD model, sdreportMC can be used to calculate samples of confidence intervals of all model parameters and REPORT()ed quantities including nonlinear functions of random effects and parameters.

Usage

sdreportMC(
  obj,
  what,
  nSamples = 1000,
  Hessian = NULL,
  CI = FALSE,
  probs = c(0.025, 0.975)
)

Arguments

Details

Caution: Currently does not work for models with fixed parameters (i.e. that use the map argument of MakeADFun.)

Value

named list corresponding to the elements of what. Each element has the structure of the corresponding quantity with an additional dimension added for the samples. For example, if a quantity is a vector, the list contains a matrix. If a quantity is a matrix, the list contains an array.

Examples

# fitting an HMM to the trex data and running sdreportMC
## negative log-likelihood function
nll = function(par) {
  getAll(par, dat) # makes everything contained available without $
  Gamma = tpm(eta) # computes transition probability matrix from unconstrained eta
  delta = stationary(Gamma) # computes stationary distribution
  # exponentiating because all parameters strictly positive
  mu = exp(logmu)
  sigma = exp(logsigma)
  kappa = exp(logkappa)
  # reporting statements for sdreportMC
  REPORT(mu)
  REPORT(sigma)
  REPORT(kappa)
  # calculating all state-dependent densities
  allprobs = matrix(1, nrow = length(step), ncol = N)
  ind = which(!is.na(step) & !is.na(angle)) # only for non-NA obs.
  for(j in 1:N){
    allprobs[ind,j] = dgamma2(step[ind],mu[j],sigma[j])*dvm(angle[ind],0,kappa[j])
  }
  -forward(delta, Gamma, allprobs) # simple forward algorithm
}

## initial parameter list
par = list(
 logmu = log(c(0.3, 1)),       # initial means for step length (log-transformed)
  logsigma = log(c(0.2, 0.7)), # initial sds for step length (log-transformed)
  logkappa = log(c(0.2, 0.7)), # initial concentration for turning angle (log-transformed)
  eta = rep(-2, 2)             # initial t.p.m. parameters (on logit scale)
)   
## data and hyperparameters
dat = list(
  step = trex$step[1:500],   # hourly step lengths
  angle = trex$angle[1:500], # hourly turning angles
  N = 2
)

## creating AD function
obj = MakeADFun(nll, par, silent = TRUE) # creating the objective function

## optimising
opt = nlminb(obj$par, obj$fn, obj$gr) # optimization

## running sdreportMC
# `mu` has report statement, `delta` is automatically reported by `forward()`
sdrMC = sdreportMC(obj, 
                   what = c("mu", "delta"), 
                   n = 50)
dim(sdrMC$delta)
# now a matrix with 50 samples (rows)

Report uncertainty of the estimated smoothing parameters or variances

Description

Computes standard deviations for the smoothing parameters of a model object returned by qreml using the delta method.

Usage

sdreport_outer(mod, invert = FALSE)

Arguments

Details

The computations are based on the approximate gradient of the restricted log likelihood. The outer Hessian is computed by finite differencing of this gradient. If the inverse smoothing parameters are requested, the standard deviations are transformed to the variances using the delta method.

Value

list containing report matrix summarising parameters and standard deviations as well as the outer Hessian matrix.

Examples

## no examples

Skew normal distribution

Description

Density, distribution function, quantile function and random generation for the skew normal distribution.

Usage

dskewnorm(x, xi = 0, omega = 1, alpha = 0, log = FALSE)

pskewnorm(q, xi = 0, omega = 1, alpha = 0, ...)

qskewnorm(p, xi = 0, omega = 1, alpha = 0, ...)

rskewnorm(n, xi = 0, omega = 1, alpha = 0)

Arguments

Details

This implementation of dskewnorm allows for automatic differentiation with RTMB while the other functions are imported from the sn package.

Value

dskewnorm gives the density, pskewnorm gives the distribution function, qskewnorm gives the quantile function, and rskewnorm generates random deviates.

Examples

x = rskewnorm(1)
d = dskewnorm(x)
p = pskewnorm(x)
q = qskewnorm(p)

Build the design and penalty matrices for smooth density estimation

Description

This high-level function can be used to prepare objects needed to estimate mixture models of smooth densities using P-Splines.

Usage

smooth_dens_construct(
  data,
  par,
  type = "real",
  k = 25,
  knots = NULL,
  degree = 3,
  diff_order = 2
)

Arguments

Details

Under the hood, make_matrices_dens is used for the actual construction of the design and penalty matrices.

You can provide one or multiple data streams of different types (real, positive, circular) and specify initial means and standard deviations/ concentrations for each data stream. This information is then converted into suitable spline coefficients. smooth_dens_construct then constructs the design and penalty matrices for standardised B-splines basis functions (integrating to one) for each data stream. For types "real" and "circular" the knots are placed equidistant in the range of the data, for type "positive" the knots are placed using polynomial spacing.

Value

a nested list containing the design matrices Z, the penalty matrices S, the initial coefficients coef the prediction design matrices Z_predict, the prediction grids xseq, and details for the basis expansion for each data stream.

Examples

## 3 data streams, each with one distribution
# normal data with mean 0 and sd 1
x1 = rnorm(100, mean = 0, sd = 1)
# gamma data with mean 5 and sd 3
x2 = rgamma2(100, mean = 5, sd = 3)
# circular data
x3 = rvm(100, mu = 0, kappa = 2)

data = data.frame(x1 = x1, x2 = x2, x3 = x3)

par = list(x1 = list(mean = 0, sd = 1),
           x2 = list(mean = 5, sd = 3),
           x3 = list(mean = 0, concentration = 2))

SmoothDens = smooth_dens_construct(data, 
                                   par,
                                   type = c("real", "positive", "circular"))
                             
# extracting objects for x1
Z1 = SmoothDens$Z$x1
S1 = SmoothDens$S$x1
coefs1 = SmoothDens$coef$x1

## one data stream, but mixture of two distributions
# normal data with mean 0 and sd 1
x = rnorm(100, mean = 0, sd = 1)
data = data.frame(x = x)

# now parameters for mixture of two normals
par = list(x = list(mean = c(0, 5), sd = c(1,1)))

SmoothDens = smooth_dens_construct(data, par = par)

# extracting objects 
Z = SmoothDens$Z$x
S = SmoothDens$S$x
coefs = SmoothDens$coef$x

Calculate conditional local state probabilities for homogeneous HMMs

Description

Computes

\Pr(S_t = j \mid X_1, ..., X_T)

for homogeneous HMMs

Usage

stateprobs(delta, Gamma, allprobs, trackID = NULL, mod = NULL)

Arguments

Value

matrix of conditional state probabilities of dimension c(n,N)

See Also

Other decoding functions: stateprobs_g(), stateprobs_p(), viterbi(), viterbi_g(), viterbi_p()

Examples

Gamma = tpm(c(-1,-2))
delta = stationary(Gamma)
allprobs = matrix(runif(10), nrow = 10, ncol = 2)

probs = stateprobs(delta, Gamma, allprobs)

Calculate conditional local state probabilities for inhomogeneous HMMs

Description

Computes

\Pr(S_t = j \mid X_1, ..., X_T)

for inhomogeneous HMMs

Usage

stateprobs_g(delta, Gamma, allprobs, trackID = NULL, mod = NULL)

Arguments

Value

matrix of conditional state probabilities of dimension c(n,N)

See Also

Other decoding functions: stateprobs(), stateprobs_p(), viterbi(), viterbi_g(), viterbi_p()

Examples

Gamma = tpm_g(runif(10), matrix(c(-1,-1,1,-2), nrow = 2, byrow = TRUE))
delta = c(0.5, 0.5)
allprobs = matrix(runif(20), nrow = 10, ncol = 2)

probs = stateprobs_g(delta, Gamma[,,-1], allprobs)

Calculate conditional local state probabilities for periodically inhomogeneous HMMs

Description

Computes

\Pr(S_t = j \mid X_1, ..., X_T)

for periodically inhomogeneous HMMs

Usage

stateprobs_p(delta, Gamma, allprobs, tod, trackID = NULL, mod = NULL)

Arguments

Value

matrix of conditional state probabilities of dimension c(n,N)

See Also

Other decoding functions: stateprobs(), stateprobs_g(), viterbi(), viterbi_g(), viterbi_p()

Examples

L = 24
beta = matrix(c(-1, 2, -1, -2, 1, -1), nrow = 2, byrow = TRUE)
Gamma = tpm_p(1:L, L, beta, degree = 1)
delta = stationary_p(Gamma, 1)
allprobs = matrix(runif(200), nrow = 100, ncol = 2)
tod = rep(1:24, 5)[1:100]

probs = stateprobs_p(delta, Gamma, allprobs, tod)

Compute the stationary distribution of a homogeneous Markov chain

Description

A homogeneous, finite state Markov chain that is irreducible and aperiodic converges to a unique stationary distribution, here called \delta. As it is stationary, this distribution satisfies

\delta \Gamma = \delta,

subject to \sum_{j=1}^N \delta_j = 1, where \Gamma is the transition probability matrix. This function solves the linear system of equations above.

Usage

stationary(Gamma)

Arguments

Value

either a single stationary distribution of the Markov chain (vector of length N) or a matrix of stationary distributions of dimension c(nTracks,N) with one stationary distribution in each row

See Also

tpm to create a transition probabilty matrix using the multinomial logistic link (softmax)

Other stationary distribution functions: stationary_cont(), stationary_p()

Examples

# single matrix
Gamma = tpm(c(rep(-2,3), rep(-3,3)))
delta = stationary(Gamma)
# multiple matrices
Gamma = array(Gamma, dim = c(3,3,10))
Delta = stationary(Gamma)

Compute the stationary distribution of a continuous-time Markov chain

Description

A well-behaved continuous-time Markov chain converges to a unique stationary distribution, here called \pi. This distribution satisfies

\pi Q = 0,

subject to \sum_{j=1}^N \pi_j = 1, where Q is the infinitesimal generator of the Markov chain. This function solves the linear system of equations above for a given generator matrix.

Usage

stationary_cont(Q)

Arguments

Value

either a single stationary distribution of the continuous-time Markov chain (vector of length N) or a matrix of stationary distributions of dimension c(nTracks,N) with one stationary distribution in each row

See Also

generator to create a generator matrix

Other stationary distribution functions: stationary(), stationary_p()

Examples

# single matrix
Q = generator(c(-2,-2))
Pi = stationary_cont(Q)
# multiple matrices
Q = array(Q, dim = c(2,2,10))
Pi = stationary_cont(Q)

Periodically stationary distribution of a periodically inhomogeneous Markov chain

Description

Computes the periodically stationary distribution of a periodically inhomogeneous Markov chain.

Usage

stationary_p(Gamma, t = NULL, ad = NULL)

Arguments

Details

If the transition probability matrix of an inhomogeneous Markov chain varies only periodically (with period length L), it converges to a so-called periodically stationary distribution. This happens, because the thinned Markov chain, which has a full cycle as each time step, has homogeneous transition probability matrix

\Gamma_t = \Gamma^{(t)} \Gamma^{(t+1)} \dots \Gamma^{(t+L-1)}

for all t = 1, \dots, L. The stationary distribution for time t satifies \delta^{(t)} \Gamma_t = \delta^{(t)}.

This function calculates said periodically stationary distribution.

Value

either the periodically stationary distribution at time t or all periodically stationary distributions.

References

Koslik, J. O., Feldmann, C. C., Mews, S., Michels, R., & Langrock, R. (2023). Inference on the state process of periodically inhomogeneous hidden Markov models for animal behavior. arXiv preprint arXiv:2312.14583.

See Also

tpm_p and tpm_g to create multiple transition matrices based on a cyclic variable or design matrix

Other stationary distribution functions: stationary(), stationary_cont()

Examples

# setting parameters for trigonometric link
beta = matrix(c(-1, 2, -1, -2, 1, -1), nrow = 2, byrow = TRUE)
Gamma = tpm_p(beta = beta, degree = 1)
# periodically stationary distribution for specific time point
delta = stationary_p(Gamma, 4)

# all periodically stationary distributions
Delta = stationary_p(Gamma)

Sparse version of stationary_p

Description

This is function computes the periodically stationary distribution of a Markov chain given a list of L sparse transition probability matrices. Compatible with automatic differentiation by RTMB

Usage

stationary_p_sparse(Gamma, t = NULL)

Arguments

Value

either the periodically stationary distribution at time t or all periodically stationary distributions.

Examples

## periodic HSMM example (here the approximating tpm is sparse)
N = 2 # number of states
L = 24 # cycle length
# time-varying mean dwell times
Z = trigBasisExp(1:L) # trigonometric basis functions design matrix
beta = matrix(c(2, 2, 0.1, -0.1, -0.2, 0.2), nrow = 2)
Lambda = exp(cbind(1, Z) %*% t(beta))
sizes = c(20, 20) # approximating chain with 40 states
# state dwell-time distributions
dm = lapply(1:N, function(i) sapply(1:sizes[i]-1, dpois, lambda = Lambda[,i]))
omega = matrix(c(0,1,1,0), nrow = N, byrow = TRUE) # embedded t.p.m.

# calculating extended-state-space t.p.m.s
Gamma = tpm_phsmm(omega, dm)
# Periodically stationary distribution for specific time point
delta = stationary_p_sparse(Gamma, 4)

# All periodically stationary distributions
Delta = stationary_p_sparse(Gamma)

Sparse version of stationary

Description

This is function computes the stationary distribution of a Markov chain with a given sparse transition probability matrix. Compatible with automatic differentiation by RTMB

Usage

stationary_sparse(Gamma)

Arguments

Value

stationary distribution of the Markov chain with the given transition probability matrix

Examples

## HSMM example (here the approximating tpm is sparse)
# building the t.p.m. of the embedded Markov chain
omega = matrix(c(0,1,1,0), nrow = 2, byrow = TRUE)
# defining state aggregate sizes
sizes = c(20, 30)
# defining state dwell-time distributions
lambda = c(5, 11)
dm = list(dpois(1:sizes[1]-1, lambda[1]), dpois(1:sizes[2]-1, lambda[2]))
# calculating extended-state-space t.p.m.
Gamma = tpm_hsmm(omega, dm)
delta = stationary_sparse(Gamma)

Summary method for qremlModel objects

Description

Prints a summary of a model object created by qreml.

Usage

## S3 method for class 'qremlModel'
summary(object, ...)

Arguments

Value

prints a summary of the model object

Examples

# no examples

Build the transition probability matrix from unconstrained parameter vector

Description

Markov chains are parametrised in terms of a transition probability matrix \Gamma, for which each row contains a conditional probability distribution of the next state given the current state. Hence, each row has entries between 0 and 1 that need to sum to one.

For numerical optimisation, we parametrise in terms of unconstrained parameters, thus this function computes said matrix from an unconstrained parameter vector via the inverse multinomial logistic link (also known as softmax) applied to each row.

Usage

tpm(param, byrow = FALSE)

Arguments

Value

Transition probability matrix of dimension c(N,N)

See Also

Other transition probability matrix functions: generator(), tpm_cont(), tpm_emb(), tpm_emb_g(), tpm_g(), tpm_p()

Examples

# 2 states: 2 free off-diagonal elements
par1 = rep(-1, 2)
Gamma1 = tpm(par1)

# 3 states: 6 free off-diagonal elements
par2 = rep(-2, 6)
Gamma2 = tpm(par2)

Calculate continuous time transition probabilities

Description

A continuous-time Markov chain is described by an infinitesimal generator matrix Q. When observing data at time points t_1, \dots, t_n the transition probabilites between t_i and t_{i+1} are caluclated as

\Gamma(\Delta t_i) = \exp(Q \Delta t_i),

where \exp() is the matrix exponential. The mapping \Gamma(\Delta t) is also called the Markov semigroup. This function calculates all transition matrices based on a given generator and time differences.

Usage

tpm_cont(Q, timediff, ad = NULL, report = TRUE)

Arguments

Value

array of continuous-time transition matrices of dimension c(N,N,n-1)

See Also

Other transition probability matrix functions: generator(), tpm(), tpm_emb(), tpm_emb_g(), tpm_g(), tpm_p()

Examples

# building a Q matrix for a 3-state cont.-time Markov chain
Q = generator(rep(-2, 6))

# draw random time differences
timediff = rexp(100, 10)

# compute all transition matrices
Gamma = tpm_cont(Q, timediff)

Build the embedded transition probability matrix of an HSMM from unconstrained parameter vector

Description

Hidden semi-Markov models are defined in terms of state durations and an embedded transition probability matrix that contains the conditional transition probabilities given that the current state is left. This matrix necessarily has diagonal entries all equal to zero as self-transitions are impossible.

This function builds such an embedded/ conditional transition probability matrix from an unconstrained parameter vector. For each row of the matrix, the inverse multinomial logistic link is applied.

For a matrix of dimension c(N,N), the number of free off-diagonal elements is N*(N-2), hence also the length of param. This means, for 2 states, the function needs to be called without any arguments, for 3-states with a vector of length 3, for 4 states with a vector of length 8, etc.

Compatible with automatic differentiation by RTMB

Usage

tpm_emb(param = NULL)

Arguments

Value

embedded/ conditional transition probability matrix of dimension c(N,N)

See Also

Other transition probability matrix functions: generator(), tpm(), tpm_cont(), tpm_emb_g(), tpm_g(), tpm_p()

Examples

# 2 states: no free off-diagonal elements
omega = tpm_emb()

# 3 states: 3 free off-diagonal elements
param = rep(0, 3)
omega = tpm_emb(param)

# 4 states: 8 free off-diagonal elements
param = rep(0, 8)
omega = tpm_emb(param)

Build all embedded transition probability matrices of an inhomogeneous HSMM

Description

Hidden semi-Markov models are defined in terms of state durations and an embedded transition probability matrix that contains the conditional transition probabilities given that the current state is left. This matrix necessarily has diagonal entries all equal to zero as self-transitions are impossible. We can allow this matrix to vary with covariates, which is the purpose of this function.

It builds all embedded/ conditional transition probability matrices based on a design and parameter matrix. For each row of the matrix, the inverse multinomial logistic link is applied.

For a matrix of dimension c(N,N), the number of free off-diagonal elements is N*(N-2) which determines the number of rows of the parameter matrix.

Compatible with automatic differentiation by RTMB

Usage

tpm_emb_g(Z, beta, report = TRUE)

Arguments

Value

array of embedded/ conditional transition probability matrices of dimension c(N,N,n)

See Also

Other transition probability matrix functions: generator(), tpm(), tpm_cont(), tpm_emb(), tpm_g(), tpm_p()

Examples

## parameter matrix for 3-state HSMM
beta = matrix(c(rep(0, 3), -0.2, 0.2, 0.1), nrow = 3)
# no intercept
Z = rnorm(100)
omega = tpm_emb_g(Z, beta)
# intercept
Z = cbind(1, Z)
omega = tpm_emb_g(Z, beta)

Build all transition probability matrices of an inhomogeneous HMM

Description

In an HMM, we often model the influence of covariates on the state process by linking them to the transition probabiltiy matrix. Most commonly, this is done by specifying a linear predictor

\eta_{ij}^{(t)} = \beta^{(ij)}_0 + \beta^{(ij)}_1 z_{t1} + \dots + \beta^{(ij)}_p z_{tp}

for each off-diagonal element (i \neq j) of the transition probability matrix and then applying the inverse multinomial logistic link (also known as softmax) to each row. This function efficiently calculates all transition probabilty matrices for a given design matrix Z and parameter matrix beta.

Usage

tpm_g(Z, beta, byrow = FALSE, ad = NULL, report = TRUE)

Arguments

Value

array of transition probability matrices of dimension c(N,N,n)

See Also

Other transition probability matrix functions: generator(), tpm(), tpm_cont(), tpm_emb(), tpm_emb_g(), tpm_p()

Examples

Z = matrix(runif(200), ncol = 2)
beta = matrix(c(-1, 1, 2, -2, 1, -2), nrow = 2, byrow = TRUE)
Gamma = tpm_g(Z, beta)

Builds the transition probability matrix of an HSMM-approximating HMM

Description

Hidden semi-Markov models (HSMMs) are a flexible extension of HMMs, where the state duration distribution is explicitly modelled. For direct numerical maximum likelhood estimation, HSMMs can be represented as HMMs on an enlarged state space (of size M) and with structured transition probabilities.

This function computes the transition matrix to approximate a given HSMM by an HMM with a larger state space.

Usage

tpm_hsmm(omega, dm, Fm = NULL, sparse = TRUE, eps = 1e-10)

Arguments

Value

extended-state-space transition probability matrix of the approximating HMM

Examples

# building the t.p.m. of the embedded Markov chain
omega = matrix(c(0,1,1,0), nrow = 2, byrow = TRUE)
# defining state aggregate sizes
sizes = c(20, 30)
# defining state dwell-time distributions
lambda = c(5, 11)
dm = list(dpois(1:sizes[1]-1, lambda[1]), dpois(1:sizes[2]-1, lambda[2]))
# calculating extended-state-space t.p.m.
Gamma = tpm_hsmm(omega, dm)

Build the transition probability matrix of an HSMM-approximating HMM

Description

Hidden semi-Markov models (HSMMs) are a flexible extension of HMMs. For direct numerical maximum likelhood estimation, HSMMs can be represented as HMMs on an enlarged state space (of size M) and with structured transition probabilities. This function computes the transition matrix of an HSMM.

Usage

tpm_hsmm2(omega, dm, eps = 1e-10)

Arguments

Value

extended-state-space transition probability matrix of the approximating HMM

Examples

# building the t.p.m. of the embedded Markov chain
omega = matrix(c(0,1,1,0), nrow = 2, byrow = TRUE)
# defining state aggregate sizes
sizes = c(20, 30)
# defining state dwell-time distributions
lambda = c(5, 11)
dm = list(dpois(1:sizes[1]-1, lambda[1]), dpois(1:sizes[2]-1, lambda[2]))
# calculating extended-state-space t.p.m.
Gamma = tpm_hsmm(omega, dm)

Builds all transition probability matrices of an inhomogeneous-HSMM-approximating HMM

Description

Hidden semi-Markov models (HSMMs) are a flexible extension of HMMs. For direct numerical maximum likelhood estimation, HSMMs can be represented as HMMs on an enlarged state space (of size M) and with structured transition probabilities.

This function computes the transition matrices of a periodically inhomogeneos HSMMs.

Usage

tpm_ihsmm(omega, dm, eps = 1e-10)

Arguments

Value

list of dimension length n - max(sapply(dm, ncol)), containing sparse extended-state-space transition probability matrices for each time point (except the first max(sapply(dm, ncol)) - 1).

Examples

N = 2
# time-varying mean dwell times
n = 100
z = runif(n)
beta = matrix(c(2, 2, 0.1, -0.1), nrow = 2)
Lambda = exp(cbind(1, z) %*% t(beta))
sizes = c(15, 15) # approximating chain with 30 states
# state dwell-time distributions
dm = lapply(1:N, function(i) sapply(1:sizes[i]-1, dpois, lambda = Lambda[,i]))

## homogeneous conditional transition probabilites
# diagonal elements are zero, rowsums are one
omega = matrix(c(0,1,1,0), nrow = N, byrow = TRUE)

# calculating extended-state-space t.p.m.s
Gamma = tpm_ihsmm(omega, dm)

## inhomogeneous conditional transition probabilites
# omega can be an array
omega = array(omega, dim = c(N,N,n))

# calculating extended-state-space t.p.m.s
Gamma = tpm_ihsmm(omega, dm)

Build all transition probability matrices of a periodically inhomogeneous HMM

Description

Given a periodically varying variable such as time of day or day of year and the associated cycle length, this function calculates the transition probability matrices by applying the inverse multinomial logistic link (also known as softmax) to linear predictors of the form

\eta^{(t)}_{ij} = \beta_0^{(ij)} + \sum_{k=1}^K \bigl( \beta_{1k}^{(ij)} \sin(\frac{2 \pi k t}{L}) + \beta_{2k}^{(ij)} \cos(\frac{2 \pi k t}{L}) \bigr)

for the off-diagonal elements (i \neq j) of the transition probability matrix. This is relevant for modeling e.g. diurnal variation and the flexibility can be increased by adding smaller frequencies (i.e. increasing K).

Usage

tpm_p(
  tod = 1:24,
  L = 24,
  beta,
  degree = 1,
  Z = NULL,
  byrow = FALSE,
  ad = NULL,
  report = TRUE
)

Arguments

Details

Note that using this function inside the negative log-likelihood function is convenient, but it performs the basis expansion into sine and cosine terms each time it is called. As these do not change during the optimisation, using tpm_g with a pre-calculated (by trigBasisExp) design matrix would be more efficient.

Value

array of transition probability matrices of dimension c(N,N,length(tod))

See Also

Other transition probability matrix functions: generator(), tpm(), tpm_cont(), tpm_emb(), tpm_emb_g(), tpm_g()

Examples

# hourly data 
tod = seq(1, 24, by = 1)
L = 24
beta = matrix(c(-1, 2, -1, -2, 1, -1), nrow = 2, byrow = TRUE)
Gamma = tpm_p(tod, L, beta, degree = 1)

# half-hourly data
## integer tod sequence
tod = seq(1, 48, by = 1)
L = 48
beta = matrix(c(-1, 2, -1, -2, 1, -1), nrow = 2, byrow = TRUE)
Gamma1 = tpm_p(tod, L, beta, degree = 1)

## equivalent specification
tod = seq(0.5, 24, by = 0.5)
L = 24
beta = matrix(c(-1, 2, -1, -2, 1, -1), nrow = 2, byrow = TRUE)
Gamma2 = tpm_p(tod, L, beta, degree = 1)

all(Gamma1 == Gamma2) # same result

Builds all transition probability matrices of an periodic-HSMM-approximating HMM

Description

Hidden semi-Markov models (HSMMs) are a flexible extension of HMMs. For direct numerical maximum likelhood estimation, HSMMs can be represented as HMMs on an enlarged state space (of size M) and with structured transition probabilities.

This function computes the transition matrices of a periodically inhomogeneos HSMMs.

Usage

tpm_phsmm(omega, dm, eps = 1e-10)

Arguments

Value

list of dimension length L, containing sparse extended-state-space transition probability matrices of the approximating HMM for each time point of the cycle.

Examples

N = 2 # number of states
L = 24 # cycle length
# time-varying mean dwell times
Z = trigBasisExp(1:L) # trigonometric basis functions design matrix
beta = matrix(c(2, 2, 0.1, -0.1, -0.2, 0.2), nrow = 2)
Lambda = exp(cbind(1, Z) %*% t(beta))
sizes = c(20, 20) # approximating chain with 40 states
# state dwell-time distributions
dm = lapply(1:N, function(i) sapply(1:sizes[i]-1, dpois, lambda = Lambda[,i]))

## homogeneous conditional transition probabilites
# diagonal elements are zero, rowsums are one
omega = matrix(c(0,1,1,0), nrow = N, byrow = TRUE)

# calculating extended-state-space t.p.m.s
Gamma = tpm_phsmm(omega, dm)

## inhomogeneous conditional transition probabilites
# omega can be an array
omega = array(omega, dim = c(N,N,L))

# calculating extended-state-space t.p.m.s
Gamma = tpm_phsmm(omega, dm)

Build all transition probability matrices of an periodic-HSMM-approximating HMM

Description

Hidden semi-Markov models (HSMMs) are a flexible extension of HMMs. For direct numerical maximum likelhood estimation, HSMMs can be represented as HMMs on an enlarged state space (of size M) and with structured transition probabilities. This function computes the transition matrices of a periodically inhomogeneos HSMMs.

Usage

tpm_phsmm2(omega, dm, eps = 1e-10)

Arguments

Value

array of dimension c(N,N,L), containing the extended-state-space transition probability matrices of the approximating HMM for each time point of the cycle.

Examples

N = 3
L = 24
# time-varying mean dwell times
Lambda = exp(matrix(rnorm(L*N, 2, 0.5), nrow = L))
sizes = c(25, 25, 25) # approximating chain with 75 states
# state dwell-time distributions
dm = list()
for(i in 1:3){
  dmi = matrix(nrow = L, ncol = sizes[i])
  for(t in 1:L){
    dmi[t,] = dpois(1:sizes[i]-1, Lambda[t,i])
  }
  dm[[i]] = dmi
}

## homogeneous conditional transition probabilites
# diagonal elements are zero, rowsums are one
omega = matrix(c(0,0.5,0.5,0.2,0,0.8,0.7,0.3,0), nrow = N, byrow = TRUE)

# calculating extended-state-space t.p.m.s
Gamma = tpm_phsmm(omega, dm)

## inhomogeneous conditional transition probabilites
# omega can be an array
omega = array(rep(omega,L), dim = c(N,N,L))
omega[1,,4] = c(0, 0.2, 0.8) # small change for inhomogeneity

# calculating extended-state-space t.p.m.s
Gamma = tpm_phsmm(omega, dm)

Compute the transition probability matrix of a thinned periodically inhomogeneous Markov chain.

Description

If the transition probability matrix of an inhomogeneous Markov chain varies only periodically (with period length L), it converges to a so-called periodically stationary distribution. This happens, because the thinned Markov chain, which has a full cycle as each time step, has homogeneous transition probability matrix

\Gamma_t = \Gamma^{(t)} \Gamma^{(t+1)} \dots \Gamma^{(t+L-1)}

for all t = 1, \dots, L. This function calculates the matrix above efficiently as a preliminery step to calculating the periodically stationary distribution.

Usage

tpm_thinned(Gamma, t)

Arguments

Value

thinned transition probabilty matrix of dimension c(N,N)

Examples

# setting parameters for trigonometric link
beta = matrix(c(-1, -2, 2, -1, 2, -4), nrow = 2, byrow = TRUE)
# calculating periodically varying t.p.m. array (of length 24 here)
Gamma = tpm_p(beta = beta)
# calculating t.p.m. of thinned Markov chain
tpm_thinned(Gamma, 4)

T-Rex Movement Data

Description

Hourly step lengths and turning angles of a Tyrannosaurus rex, living 66 million years ago.

Usage

trex

Format

A data frame with 10.000 rows and 4 variables:

tod

time of day variable ranging from 1 to 24

step

hourly step lengths in kilometres

angle

hourly turning angles in radians

state

hidden state variable

Source

Generated for example purposes.

Compute the design matrix for a trigonometric basis expansion

Description

Given a periodically varying variable such as time of day or day of year and the associated cycle length, this function performs a basis expansion to efficiently calculate a linear predictor of the form

\eta^{(t)} = \beta_0 + \sum_{k=1}^K \bigl( \beta_{1k} \sin(\frac{2 \pi k t}{L}) + \beta_{2k} \cos(\frac{2 \pi k t}{L}) \bigr).

This is relevant for modeling e.g. diurnal variation and the flexibility can be increased by adding smaller frequencies (i.e. increasing K).

Usage

trigBasisExp(tod, L = 24, degree = 1)

Arguments

Value

design matrix (without intercept column), ordered as sin1, cos1, sin2, cos2, ...

Examples

## hourly data
tod = rep(1:24, 10)
Z = trigBasisExp(tod, L = 24, degree = 2)

## half-hourly data
tod = rep(1:48/2, 10) # in [0,24] -> L = 24
Z1 = trigBasisExp(tod, L = 24, degree = 3)

tod = rep(1:48, 10) # in [1,48] -> L = 48
Z2 = trigBasisExp(tod, L = 48, degree = 3)

all(Z1 == Z2)
# The latter two are equivalent specifications!

Viterbi algorithm for state decoding in homogeneous HMMs

Description

The Viterbi algorithm allows one to decode the most probable state sequence of an HMM.

Usage

viterbi(delta, Gamma, allprobs, trackID = NULL, mod = NULL)

Arguments

Value

vector of decoded states of length n

See Also

Other decoding functions: stateprobs(), stateprobs_g(), stateprobs_p(), viterbi_g(), viterbi_p()

Examples

delta = c(0.5, 0.5)
Gamma = matrix(c(0.9, 0.1, 0.2, 0.8), nrow = 2, byrow = TRUE)
allprobs = matrix(runif(200), nrow = 100, ncol = 2)
states = viterbi(delta, Gamma, allprobs)

Viterbi algorithm for state decoding in inhomogeneous HMMs

Description

The Viterbi algorithm allows one to decode the most probable state sequence of an HMM.

Usage

viterbi_g(delta, Gamma, allprobs, trackID = NULL, mod = NULL)

Arguments

Value

vector of decoded states of length n

See Also

Other decoding functions: stateprobs(), stateprobs_g(), stateprobs_p(), viterbi(), viterbi_p()

Examples

delta = c(0.5, 0.5)
Gamma = tpm_g(runif(10), matrix(c(-2,-2,1,-1), nrow = 2))
allprobs = matrix(runif(20), nrow = 10, ncol = 2)
states = viterbi_g(delta, Gamma[,,-1], allprobs)

Viterbi algorithm for state decoding in periodically inhomogeneous HMMs

Description

The Viterbi algorithm allows one to decode the most probable state sequence of an HMM.

Usage

viterbi_p(delta, Gamma, allprobs, tod, trackID = NULL, mod = NULL)

Arguments

Value

vector of decoded states of length n

See Also

Other decoding functions: stateprobs(), stateprobs_g(), stateprobs_p(), viterbi(), viterbi_g()

Examples

delta = c(0.5, 0.5)
beta = matrix(c(-2, 1, -1,
                -2, -1, 1), nrow = 2, byrow = TRUE)
Gamma = tpm_p(1:24, 24, beta)

tod = rep(1:24, 5)
n = length(tod)

allprobs = matrix(runif(2*n), nrow = n, ncol = 2)
states = viterbi_p(delta, Gamma, allprobs, tod)

von Mises distribution

Description

Density, distribution function and random generation for the von Mises distribution.

Usage

dvm(x, mu = 0, kappa = 1, log = FALSE)

pvm(q, mu = 0, kappa = 1, from = NULL, tol = 1e-20)

rvm(n, mu = 0, kappa = 1, wrap = TRUE)

Arguments

Details

This implementation of dvm allows for automatic differentiation with RTMB. rvm and pvm are simply wrappers of the corresponding functions from circular.

Value

dvm gives the density, pvm gives the distribution function, and rvm generates random deviates.

Examples

set.seed(1)
x = rvm(10, 0, 1)
d = dvm(x, 0, 1)
p = pvm(x, 0, 1)

wrapped Cauchy distribution

Description

Density and random generation for the wrapped Cauchy distribution.

Usage

dwrpcauchy(x, mu = 0, rho, log = FALSE)

rwrpcauchy(n, mu = 0, rho, wrap = TRUE)

Arguments

Details

This implementation of dwrpcauchy allows for automatic differentiation with RTMB. rwrpcauchy is simply a wrapper for rwrappedcauchyimported from circular.

Value

dwrpcauchy gives the density and rwrpcauchy generates random deviates.

Examples

set.seed(1)
x = rwrpcauchy(10, 0, 1)
d = dwrpcauchy(x, 0, 1)

Zero-inflated density constructer

Description

Constructs a zero-inflated density function from a given probability density function

Usage

zero_inflate(dist, discrete = NULL)

Arguments

Details

The definition of zero-inflation is different for discrete and continuous distributions. For discrete distributions with p.m.f. f and zero-inflation probability p, we have

\Pr(X = 0) = p + (1 - p) \cdot f(0),

and

\Pr(X = x) = (1 - p) \cdot f(x), \quad x > 0.

For continuous distributions with p.d.f. f, we have

f_{\text{zinfl}}(x) = p \cdot \delta_0(x) + (1 - p) \cdot f(x),

where \delta_0 is the Dirac delta function at zero.

Value

zero-inflated density function with first argument x, second argument zeroprob, and additional arguments ... that will be passed to dist.

Examples

dzinorm <- zero_inflate(dnorm)
dzinorm(c(NA, 0, 2), 0.5, mean = 1, sd = 1)

zipois <- zero_inflate(dpois)
zipois(c(NA, 0, 1), 0.5, 1)