Package 'LaMa'

Description

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

Usage

buildSmoothDens(data, type = "real", par, k = 20, 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. buildSmoothDens 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 = buildSmoothDens(data, 
                             type = c("real", "positive", "circular"),
                             par)
                             
# 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 = buildSmoothDens(data, par = par)

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

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)

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)

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)

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)

Arguments

Value

log-likelihood for given data and parameters

See Also

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

Examples

## 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(-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], 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

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

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
)

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

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)

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)

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))

Description

Build the design matrix 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, the penalty matrix S, the formula, the data and the knots

Examples

modmat = make_matrices(~ s(x), data.frame(x = 1:10))

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,
  npoints = 10000,
  diff_order = 2,
  pow = 0.5
)

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

modmat = make_matrices_dens(x = (-50):50, k = 20)
modmat = make_matrices_dens(x = 1:100, k = 20, type = "positive")
modmat = make_matrices_dens(x = seq(-pi,pi), k = 20, type = "circular")

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.

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.

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), 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, ad = TRUE) +
      penalty(betaspline, S, lambda) # this does all the penalization work
}

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

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)

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))

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)
stateprobs = matrix(0.5, nrow = 100, ncol = 2)
par = list(lambda = c(1,2))
pres = pseudo_res(obs, "pois", par, stateprobs)

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)

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,
  psname = "lambda",
  alpha = 0,
  smoothing = 1,
  maxiter = 100,
  tol = 1e-04,
  control = list(reltol = 1e-10, maxit = 1000),
  silent = 1,
  joint_unc = TRUE,
  saveall = FALSE
)

Arguments

Value

returns a model list influenced by the users report statements in pnll

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(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, ad = TRUE) +
      penalty(betaspline, S, lambda) # this does all the penalization work
}

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

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,
  Hessian = NULL,
  CI = FALSE,
  n = 1000,
  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)

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(200), nrow = 100, ncol = 2)

probs = stateprobs(delta, Gamma, allprobs)

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(99), matrix(c(-1,-1,1,-2), nrow = 2, byrow = TRUE))
delta = c(0.5, 0.5)
allprobs = matrix(runif(200), nrow = 100, ncol = 2)

probs = stateprobs_g(delta, Gamma, allprobs)

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)

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

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

See Also

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

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

Examples

Gamma = tpm(c(rep(-2,3), rep(-3,3)))
delta = stationary(Gamma)

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

stationary distribution of the continuous-time Markov chain with given generator matrix

See Also

generator to create a generator matrix

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

Examples

Q = generator(c(-2,-2))
Pi = stationary_cont(Q)

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.$ The stationary distribution for time $t$ satifies $\delta^{(t)} \Gamma_t = \delta^{(t)}$.

This function calculates said periodically stationary distribution.

Usage

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

Arguments

Value

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

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)

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)

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)

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)

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)

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)

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)

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)

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)

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

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)

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)

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)

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.

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!