Movatterモバイル変換

Type:	Package
Title:	Fast Numerical Maximum Likelihood Estimation for Latent MarkovModels
Version:	2.0.6
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-09-23 13:24:50 UTC; jan-ole
Author:	Jan-Ole Koslik [aut, cre]
Maintainer:	Jan-Ole Koslik <jan-ole.koslik@uni-bielefeld.de>
Repository:	CRAN
Date/Publication:	2025-09-23 13:50:02 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 Koslikjan-ole.koslik@uni-bielefeld.de (ORCID)

See Also

Useful links:

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

Sparsity-retaining matrix multiplication

Description

Standard matrix multiplication destroys automatic sparsity detection byRTMB which is essential for models with high-dimensional random effects.This can be mitigated by changing to "plain" withTapeConfig, but this can make AD tape construction very slow.Here, we provide a different version that retains sparsity. It may be slightly slower than the standard method when constructing the AD tape.

Usage

A %sp% B

Arguments

Value

the matrix product of A and B, which is of dimension n x m

Examples

A <- matrix(1:6, nrow = 2, ncol = 3)B <- matrix(7:12, nrow = 3, ncol = 2)A %sp% B

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 tomake_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 ofperiod).

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 matrixX = 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 temperatureform = ~ 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 byexpressing the transition probability matrix as a periodic function of the time of day usingtpm_p orcosinor, 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) iftime is specified, or overall dwell-time distribution iftime isNULL. If more than onestate 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 linkbeta = 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 stateddwell(1:20, Gamma, time = 1:4, state = 1)# results in 4x20 matrix# or overall distribution for all statesddwell(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 thejoint log-likelihood with penalised splines or i.i.d. random effects that have a multivariate Gaussian distribution with fixed precision/ penalty matrix\lambda S.AsS is fixed and only scaled by\lambda, it is more efficient to precompute the determinant ofS (for the normalisation constant) and only scale the quadratic form by\lambdawhen multiple spline parameters/ random effects with different\lambda's but the same penalty matrixS are evaluated.

Usage

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

Arguments

Details

This implementation allows for automatic differentiation withRTMB.

Value

vector of density values

Examples

x = matrix(runif(30), nrow = 3)# iid random effectsS = diag(10)sigma = c(1, 2, 3) # random effect standard deviationslambda = 1 / sigma^2d = dgmrf2(x, 0, S, lambda)# P-splinesL = diff(diag(10), diff = 2) # second-order difference matrixS = t(L) %*% Llambda = 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 ofddirichlet allows for automatic differentiation withRTMB.

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 theforward 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 functionnll = 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 datapar = 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])

Generalforward 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 theforward 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 usageGamma = 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 functionnll = 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 datapar = 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 modifiedforward 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 sizeM) and with structured transition probabilities.

This function is designed to be used with automatic differentiation based on theR packageRTMB. 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 modifiedforward 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 sizeM) and with structured transition probabilities.

This function is designed to be used with automatic differentiation based on theR packageRTMB. 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 theforward 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 onlyL unique matrices ifL 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 theseL 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 functionnll = 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 datapar = 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 sizeM) 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 overforward_ihsmm because it computes thecorrect periodically stationary distribution and no observations are lost for the approximation.

This function is designed to be used with automatic differentiation based on theR packageRTMB. 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 modifiedforward 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 sizeM) 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 HSMMmu = 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 2s = 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 functionmllk = 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 datatheta.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 sizeM) 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 HSMMmu = c(0, 6)beta = matrix(c(log(4),log(6),-0.2,0.2,-0.1,0.4), nrow=2)# time varying mean dwell timeLambda = 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 2s = rep(1:2, 100)C = x = numeric(0)tod = rep(1:24, 50) # time of day variabletime = 1for(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 functionmllk = 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 datatheta.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 forthe 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 withRTMB.

Value

dgamma2 gives the density,pgamma2 gives the distribution function,qgamma2 gives the quantile function, andrgamma2 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 ofx

Build the generator matrix of a continuous-time Markov chain

Description

This function builds theinfinitesimal generator matrix for acontinuous-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_g2(),tpm_p()

Examples

# 2 states: 2 free off-diagonal elementsgenerator(rep(-1, 2))# 3 states: 6 free off-diagonal elementsgenerator(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 of classLaMa_matrices containing:

See Also

predict.LaMa_matrices for prediction design matrix construction based on the model matrices object created by this function.

Examples

data = data.frame(x = runif(100),                   y = runif(100),                  g = factor(rep(1:10, each = 10)))# unvariate thin plate regression splinemodmat = make_matrices(~ s(x), data)# univariate P-splinemodmat = make_matrices(~ s(x, bs = "ps"), data)# adding random interceptmodmat = make_matrices(~ s(g, bs = "re") + s(x, bs = "ps"), data)# tensorproduct of x and ymodmat = make_matrices(~ s(x) + s(y) + ti(x,y), data)# multiple formulas at oncemodmat = make_matrices(list(mu ~ s(x) + y, sigma ~ s(g, bs = "re")), data = 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 matrixZ, the penalty matrixS, the prediction design matrixZ_predict, the prediction gridxseq, and details for the basis expansion.

Examples

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

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_old(formula, data, knots = NULL)

Arguments

Value

a list containing the design matrixZ, a (potentially nested) list of penalty matricesS, theformula, thedata, theknots, and the originalmod object returned bymgcv.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 splinemodmat = make_matrices(~ s(x), data)# univariate P-splinemodmat = make_matrices(~ s(x, bs = "ps"), data)# adding random interceptmodmat = make_matrices(~ s(g, bs = "re") + s(x, bs = "ps"), data)# tensorproduct of x and ymodmat = make_matrices(~ s(x) + s(y) + ti(x,y), data)

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 usingRTMB.

Usage

min2(x, y)max2(x, y)

Arguments

Value

min2 returns the parallel minimum andmax2 the parallel maximum ofx andy

Examples

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

Smooth approximations to max(x, 0) and min(x, 0)

Description

Smooth approximations to max(x, 0) and min(x, 0)

Usage

max0_smooth(x, rho = 20)min0_smooth(x, rho = 20)

Arguments

Value

the approximate maximum or minimum of x and 0

Examples

x <- seq(-1, 1, by = 0.1)min0_smooth(x)max0_smooth(x)

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 vectorsb_i, and fixed penalty matricesS_i.

It is intended to be used inside thepenalised negative log-likelihood function when fitting models with penalised splines or simple random effects viaquasi restricted maximum likelihood (qREML) with theqreml function.Forqreml to work, the likelihood function needs to be compatible with theRTMB R package to enable automatic differentiation.

Usage

penalty(re_coef, S, lambda)

Arguments

Details

Caution: The formatting ofre_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 topenalty.If these are seperate random effects, each with its own name, they need to be passed as a list topenalty. Moreover, the ordering ofre_coef needs to match the character vectorrandom specified inqreml.

Value

returns the penalty value and reports toqreml.

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 theqREML algorithm

Examples

# Example with a single random effectre = rep(0, 5)S = diag(5)lambda = 1penalty(re, S, lambda)# Example with two random effects, # where one element contains two random effects of similar structurere = list(matrix(0, 2, 5), rep(0, 4))S = list(diag(5), diag(4))lambda = c(1,1,2) # length = total number of random effectspenalty(re, S, lambda)# Full model-fitting exampledata = trex[1:1000,] # subset# initial parameter listpar = 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 lambdadat = 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 matricesmodmat = make_matrices(~ s(tod, bs = "cp"),                        data = data.frame(tod = 1:24),                        knots = list(tod = c(0,24))) # wrapping pointsdat$Z = modmat$Z # spline design matrixdat$S = modmat$S # penalty matrix# penalised negative log-likelihood functionpnll = 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 fittingmod = 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 vectorb, and fixed penalty matricesS_i.This generalises thepenalty by allowing subsets of the coefficient vectorb to be penalised multiple times with different smoothing parameters, which is necessary fortensor products,functional random effects oradaptive smoothing.

It is intended to be used inside thepenalised negative log-likelihood function when fitting models with penalised splines or simple random effects viaquasi restricted maximum likelihood (qREML) with theqreml function.Forqreml to work, the likelihood function needs to be compatible with theRTMB R package to enable automatic differentiation.

Usage

penalty2(re_coef, S, lambda)

Arguments

Details

Caution: The formatting ofre_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 topenalty.If these are seperate random effects, each with its own name, they need to be passed as a list topenalty. Moreover, the ordering ofre_coef needs to match the character vectorrandom specified inqreml.

Value

returns the penalty value and reports toqreml.

See Also

qreml for theqREML algorithm

Examples

# Example with a single random effectre = rep(0, 5)S = diag(5)lambda = 1penalty(re, S, lambda)# Example with two random effects, # where one element contains two random effects of similar structurere = list(matrix(0, 2, 5), rep(0, 4))S = list(diag(5), diag(4))lambda = c(1,1,2) # length = total number of random effectspenalty(re, S, lambda)# Full model-fitting exampledata = trex[1:1000,] # subset# initial parameter listpar = 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 lambdadat = 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 matricesmodmat = make_matrices(~ s(tod, bs = "cp"),                        data = data.frame(tod = 1:24),                        knots = list(tod = c(0,24))) # wrapping pointsdat$Z = modmat$Z # spline design matrixdat$S = modmat$S # penalty matrix# penalised negative log-likelihood functionpnll = 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 fittingmod = qreml(pnll, par, dat, random = "betaspline")

Penalty approximation of unimodality constraints for univariates smooths

Description

Penalty approximation of unimodality constraints for univariates smooths

Usage

penalty_uni(coef, m, kappa = 1000, concave = TRUE, rho = 20)

Arguments

Value

a numeric value of the penalty for the given coefficients

Examples

## coefficient vectorcoef <- c(1, 2, 3, 2, 1)# mode at position 3penalty_uni(coef, m = 3) # basically zero#' # mode at position 2penalty_uni(coef, m = 2) # large positive penalty## coefficient matrixcoef <- rbind(coef, coef)m <- c(1, 4)penalty_uni(coef, m)

Plot pseudo-residuals

Description

Plot pseudo-residuals computed bypseudo_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 datastep = 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)plot(pres, hist = TRUE)

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

Description

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

Usage

pred_matrix(model_matrices, newdata, what = NULL, exclude = NULL)

Arguments

Value

prediction design matrix fornewdata with the same basis as used formodel_matrices

Examples

# single formulamodmat = make_matrices(~ s(x), data.frame(x = 1:10))Z_p = pred_matrix(modmat, data.frame(x = 1:10 - 0.5))# with multiple formulasmodmat = make_matrices(list(mu ~ s(x), sigma ~ s(x, bs = "ps")), data = data.frame(x = 1:10))Z_p = pred_matrix(modmat, data.frame(x = 1:10 - 0.5), what = "mu")# nested formula listform = list(stream1 = list(mu ~ s(x), sigma ~ s(x, bs = "ps")))modmat = make_matrices(form, data = data.frame(x = 1:10))Z_p = pred_matrix(modmat, data.frame(x = 1:10 - 0.5), what = c("stream1", "mu"))

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

Description

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

Usage

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

Arguments

Value

prediction design matrix fornewdata with the same basis as used formodel_matrices

See Also

make_matrices for creating objects of classLaMa_matrices which can be used for prediction by this function.

Examples

# single formulamodmat = make_matrices(~ s(x), data.frame(x = 1:10))Z_p = predict(modmat, data.frame(x = 1:10 - 0.5))# with multiple formulasmodmat = make_matrices(list(mu ~ s(x), sigma ~ s(x, bs = "ps")), data = data.frame(x = 1:10))Z_p = predict(modmat, data.frame(x = 1:10 - 0.5), what = "mu")# nested formula listform = list(stream1 = list(mu ~ s(x), sigma ~ s(x, bs = "ps")))modmat = make_matrices(form, data = data.frame(x = 1:10))Z_p = predict(modmat, data.frame(x = 1:10 - 0.5), what = c("stream1", "mu"))

Process and standardise formulas for the state process of hidden Markov models

Description

Process and standardise formulas for the state process of hidden Markov models

Usage

process_hid_formulas(formulas, nStates, ref = NULL)

Arguments

Value

named list of formulas of lengthnStates * (nStates - 1), where each formula corresponds to a transition from statei to statej, excluding transitions fromi toref[i].

Examples

# single formula for all non-reference category elementsformulas = process_hid_formulas(~ s(x), nStates = 3)# now a list of length 6 with names tr.ij, not including reference categories# different reference categoriesformulas = process_hid_formulas(~ s(x), nStates = 3, ref = c(1,1,1))# different formulas for different entries (and only for 2 of 6)formulas = list(tr.12 ~ s(x), tr.23 ~ s(y))formulas = process_hid_formulas(formulas, nStates = 3, ref = c(1,1,1))# also a list of length 6, remaining entries filled with tr.ij ~ 1# matrix input with reference categoriesformulas = matrix(c(".", "~ s(x)", "~ s(y)",                    "~ g", ".", "~ I(x^2)",                    "~ y", "~ 1", "."),                     nrow = 3, byrow = TRUE)# dots define reference categoriesformulas = process_hid_formulas(formulas, nStates = 3)

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 bystateprobs orstateprobs_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 forpseudo_res_discrete.

Value

vector of pseudo residuals

See Also

plot.LaMaResiduals for plotting pseudo-residuals.

Examples

## continuous-valued observationsobs = 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 observationsobs = rpois(100, lambda = 1)par = list(lambda = c(1,2))pres = pseudo_res(obs, "pois", par, stateprobs)## custom CDF functionobs = 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## no CDF available, only density (artificial example)obs = rnorm(100)par = list(mean = c(1,2), sd = c(1,1))cdf = function(x, mean, sd) integrate(dnorm, -Inf, x, mean = mean, sd = sd)$valuepres = pseudo_res(obs, cdf, par, stateprobs)## full example with model objectstep = 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)

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 fordiscrete-valued observations, based on the local state probabilities obtained bystateprobs orstateprobs_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.Ifrandomise is set toFALSE, 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 involvepenalised splines or simplei.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 vectorsb_i, and fixed penalty matricesS_i.

TheqREML 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 packageRTMB for automatic differentiation in the inner optimisation.The user has to specify thepenalised negative log-likelihood functionpnll structured as dictated byRTMB and use thepenalty function to compute the quadratic-form penalty inside the likelihood.

Usage

qreml(  pnll,  par,  dat,  random,  map = NULL,  silent = 1,  psname = "lambda",  alpha = 0.3,  smoothing = 1,  maxiter = 100,  tol = 1e-04,  method = "BFGS",  control = list(),  conv_crit = "relchange",  joint_unc = FALSE,  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 andpenalty2 to compute the penalty inside the likelihood function

Examples

data = trex[1:1000,] # subset# initial parameter listpar = 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 lambdadat = 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 matricesmodmat = make_matrices(~ s(tod, bs = "cp"),                        data = data.frame(tod = 1:24),                        knots = list(tod = c(0,24))) # wrapping pointsdat$Z = modmat$Z # spline design matrixdat$S = modmat$S # penalty matrix# penalised negative log-likelihood functionpnll = 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 fittingmod = 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 involvepenalised splines or simplei.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 vectorsb_i, and fixed penalty matricesS_i.

TheqREML 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 packageRTMB for automatic differentiation in the inner optimisation.The user has to specify thepenalised negative log-likelihood functionpnll structured as dictated byRTMB and use thepenalty 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 listpar = 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 lambdadat = 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 matricesmodmat = make_matrices(~ s(tod, bs = "cp"),                        data = data.frame(tod = 1:24),                        knots = list(tod = c(0,24))) # wrapping pointsdat$Z = modmat$Z # spline design matrixdat$S = modmat$S # penalty matrix# penalised negative log-likelihood functionpnll = 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 fittingmod = qreml_old(pnll, par, dat, random = "betaspline")

Monte Carlo version ofsdreport

Description

After optimisation of an AD model,sdreportMC can be used to calculate samples of confidence intervals of all model parameters andREPORT()ed quantitiesincluding 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

This function simply samples from the approximate multivariate normal distribution of the maximum likelihood estimate (MLE) of the parameters

\hat{\theta} \sim N(\theta, H^{-1}),

whereH is the Hessian matrix of the negative log-likelihood function at the MLE. It then returns either the sampled parameters orREPORT()ed transformations of them. IfCI = TRUE, it does not return the samples but directly computes confidence intervals.

If you are interested in several quantities, callingsdreportMC once with a vectorwhat will generally be faster than calling it several times with single elements ofwhat.

Value

named list corresponding to the elements ofwhat. 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. If quantity is an array, the list contains an array with one extra dimension.

Examples

# fitting an HMM to the trex data and running sdreportMC## negative log-likelihood functionnll = 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 listpar = 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 hyperparametersdat = list(  step = trex$step[1:500],   # hourly step lengths  angle = trex$angle[1:500], # hourly turning angles  N = 2)## creating AD functionobj = MakeADFun(nll, par, silent = TRUE) # creating the objective function## optimisingopt = 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"),                    nSamples = 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 byqreml 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 containingreport matrix summarising parameters and standard deviations as well as the outerHessian matrix.

Examples

## no examples

Skew normal distribution

Description

Density, distribution function, quantile function and random generation forthe 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 ofdskewnorm allows for automatic differentiation withRTMB while the other functions are imported from thesn package.

Value

dskewnorm gives the density,pskewnorm gives the distribution function,qskewnorm gives the quantile function, andrskewnorm 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 matricesZ, the penalty matricesS, the initial coefficientscoef the prediction design matricesZ_predict, the prediction gridsxseq, 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 1x1 = rnorm(100, mean = 0, sd = 1)# gamma data with mean 5 and sd 3x2 = rgamma2(100, mean = 5, sd = 3)# circular datax3 = 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 x1Z1 = SmoothDens$Z$x1S1 = SmoothDens$S$x1coefs1 = SmoothDens$coef$x1## one data stream, but mixture of two distributions# normal data with mean 0 and sd 1x = rnorm(100, mean = 0, sd = 1)data = data.frame(x = x)# now parameters for mixture of two normalspar = list(x = list(mean = c(0, 5), sd = c(1,1)))SmoothDens = smooth_dens_construct(data, par = par)# extracting objects Z = SmoothDens$Z$xS = SmoothDens$S$xcoefs = 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 = 24beta = 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 lengthN) or a matrix of stationary distributions of dimensionc(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 matrixGamma = tpm(c(rep(-2,3), rep(-3,3)))delta = stationary(Gamma)# multiple matricesGamma = 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,whereQ 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 lengthN) or a matrix of stationary distributions of dimensionc(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 matrixQ = generator(c(-2,-2))Pi = stationary_cont(Q)# multiple matricesQ = 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 lengthL), 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 allt = 1, \dots, L.The stationary distribution for timet 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 andtpm_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 linkbeta = 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 pointdelta = stationary_p(Gamma, 4)# all periodically stationary distributionsDelta = stationary_p(Gamma)

Sparse version ofstationary_p

Description

This is function computes the periodically stationary distribution of a Markov chain given a list of Lsparse transition probability matrices.Compatible with automatic differentiation byRTMB

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 statesL = 24 # cycle length# time-varying mean dwell timesZ = trigBasisExp(1:L) # trigonometric basis functions design matrixbeta = 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 distributionsdm = 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.sGamma = tpm_phsmm(omega, dm)# Periodically stationary distribution for specific time pointdelta = stationary_p_sparse(Gamma, 4)# All periodically stationary distributionsDelta = stationary_p_sparse(Gamma)

Sparse version ofstationary

Description

This is function computes the stationary distribution of a Markov chain with a givensparse transition probability matrix.Compatible with automatic differentiation byRTMB

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 chainomega = matrix(c(0,1,1,0), nrow = 2, byrow = TRUE)# defining state aggregate sizessizes = c(20, 30)# defining state dwell-time distributionslambda = 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 forqremlModel objects

Description

Prints a summary of a model object created byqreml.

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_g2(),tpm_p()

Examples

# 2 states: 2 free off-diagonal elementspar1 = rep(-1, 2)Gamma1 = tpm(par1)# 3 states: 6 free off-diagonal elementspar2 = rep(-2, 6)Gamma2 = tpm(par2)

Calculate continuous time transition probabilities

Description

A continuous-time Markov chain is described by an infinitesimal generator matrixQ. When observing data at time pointst_1, \dots, t_n the transition probabilites betweent_i andt_{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 theMarkov 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_g2(),tpm_p()

Examples

# building a Q matrix for a 3-state cont.-time Markov chainQ = generator(rep(-2, 6))# draw random time differencestimediff = rexp(100, 10)# compute all transition matricesGamma = 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 anembedded transition probability matrix that contains the conditional transition probabilities given that thecurrent 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 ofparam.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 byRTMB

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_g2(),tpm_p()

Examples

# 2 states: no free off-diagonal elementsomega = tpm_emb()# 3 states: 3 free off-diagonal elementsparam = rep(0, 3)omega = tpm_emb(param)# 4 states: 8 free off-diagonal elementsparam = 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 anembedded transition probability matrix that contains the conditional transition probabilities given that thecurrent 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 byRTMB

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_g2(),tpm_p()

Examples

## parameter matrix for 3-state HSMMbeta = matrix(c(rep(0, 3), -0.2, 0.2, 0.1), nrow = 3)# no interceptZ = rnorm(100)omega = tpm_emb_g(Z, beta)# interceptZ = 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 matrixZ and parameter matrixbeta.

Usage

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

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_g2(),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)

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 linear predictors

\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 matrixZ and parameter matrixbeta.

Usage

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

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_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 sizeM) 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 chainomega = matrix(c(0,1,1,0), nrow = 2, byrow = TRUE)# defining state aggregate sizessizes = c(20, 30)# defining state dwell-time distributionslambda = 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 sizeM) 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 chainomega = matrix(c(0,1,1,0), nrow = 2, byrow = TRUE)# defining state aggregate sizessizes = c(20, 30)# defining state dwell-time distributionslambda = 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 sizeM) 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 lengthn - max(sapply(dm, ncol)), containing sparse extended-state-space transition probability matrices for each time point (except the firstmax(sapply(dm, ncol)) - 1).

Examples

N = 2# time-varying mean dwell timesn = 100z = 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 distributionsdm = lapply(1:N, function(i) sapply(1:sizes[i]-1, dpois, lambda = Lambda[,i]))## homogeneous conditional transition probabilites# diagonal elements are zero, rowsums are oneomega = matrix(c(0,1,1,0), nrow = N, byrow = TRUE)# calculating extended-state-space t.p.m.sGamma = tpm_ihsmm(omega, dm)## inhomogeneous conditional transition probabilites# omega can be an arrayomega = array(omega, dim = c(N,N,n))# calculating extended-state-space t.p.m.sGamma = 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. increasingK).

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, usingtpm_g with a pre-calculated (bytrigBasisExp) 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(),tpm_g2()

Examples

# hourly data tod = seq(1, 24, by = 1)L = 24beta = 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 sequencetod = seq(1, 48, by = 1)L = 48beta = matrix(c(-1, 2, -1, -2, 1, -1), nrow = 2, byrow = TRUE)Gamma1 = tpm_p(tod, L, beta, degree = 1)## equivalent specificationtod = seq(0.5, 24, by = 0.5)L = 24beta = 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 sizeM) 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 statesL = 24 # cycle length# time-varying mean dwell timesZ = trigBasisExp(1:L) # trigonometric basis functions design matrixbeta = 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 distributionsdm = lapply(1:N, function(i) sapply(1:sizes[i]-1, dpois, lambda = Lambda[,i]))## homogeneous conditional transition probabilites# diagonal elements are zero, rowsums are oneomega = matrix(c(0,1,1,0), nrow = N, byrow = TRUE)# calculating extended-state-space t.p.m.sGamma = tpm_phsmm(omega, dm)## inhomogeneous conditional transition probabilites# omega can be an arrayomega = array(omega, dim = c(N,N,L))# calculating extended-state-space t.p.m.sGamma = 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 sizeM) 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 = 3L = 24# time-varying mean dwell timesLambda = exp(matrix(rnorm(L*N, 2, 0.5), nrow = L))sizes = c(25, 25, 25) # approximating chain with 75 states# state dwell-time distributionsdm = 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 oneomega = 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.sGamma = tpm_phsmm(omega, dm)## inhomogeneous conditional transition probabilites# omega can be an arrayomega = 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.sGamma = 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 lengthL), 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 allt = 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 linkbeta = 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 chaintpm_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. increasingK).

Usage

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

Arguments

Value

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

Examples

## hourly datatod = rep(1:24, 10)Z = trigBasisExp(tod, L = 24, degree = 2)## half-hourly datatod = rep(1:48/2, 10) # in [0,24] -> L = 24Z1 = trigBasisExp(tod, L = 24, degree = 3)tod = rep(1:48, 10) # in [1,48] -> L = 48Z2 = 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 ofdvm allows for automatic differentiation withRTMB.rvm andpvm are simply wrappers of the corresponding functions fromcircular.

Value

dvm gives the density,pvm gives the distribution function, andrvm 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 ofdwrpcauchy allows for automatic differentiation withRTMB.rwrpcauchy is simply a wrapper forrwrappedcauchyimported fromcircular.

Value

dwrpcauchy gives the density andrwrpcauchy 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 probabilityp, 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 argumentx, second argumentzeroprob, and additional arguments... that will be passed todist.

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)