R/RcppExports.R

In RprobitB: Bayesian Probit Choice Modeling

# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

#' Update latent classes
#'
#' @description
#' Low-level sampler kernels for class weights, allocations, sizes, the
#' weight-based class update, and Dirichlet-process class updates.
#'
#' @param prob \[`numeric(C)`\]\cr
#' Class probabilities.
#'
#' @param delta \[`numeric(1)`\]\cr
#' Dirichlet concentration parameter.
#'
#' @param m \[`numeric(C)`\]\cr
#' Class sizes.
#'
#' @param s \[`numeric(C)`\]\cr
#' Class weights.
#'
#' @param beta \[`matrix(P, N)`\]\cr
#' Individual coefficient draws in columns.
#'
#' @param b \[`matrix(P, C)`\]\cr
#' Class means in columns.
#'
#' @param Omega \[`matrix(P * P, C)`\]\cr
#' Vectorized class covariance matrices in columns.
#'
#' @param C \[`integer(1)`\]\cr
#' Number of classes.
#'
#' @param z \[`numeric(N)`\]\cr
#' Class allocations numbered from one to `C`.
#'
#' @param non_zero \[`logical(1)`\]\cr
#' Replace empty class sizes by one?
#'
#' @param identify_classes \[`logical(1)`\]\cr
#' Order the current active classes by size?
#'
#' @param Cmax \[`integer(1)`\]\cr
#' Maximum number of classes.
#'
#' @param epsmin \[`numeric(1)`\]\cr
#' Remove the smallest class when its weight is below this threshold.
#'
#' @param epsmax \[`numeric(1)`\]\cr
#' Split the largest class when its weight exceeds this threshold.
#'
#' @param deltamin \[`numeric(1)`\]\cr
#' Merge the two closest classes when the Euclidean distance between their
#' means is below this threshold.
#'
#' @param deltashift \[`numeric(1)`\]\cr
#' Scale of the mean displacement along the leading covariance eigenvector
#' after splitting a class.
#'
#' @param mu_b_0 \[`numeric(P)`\]\cr
#' Prior mean for class means.
#'
#' @param Sigma_b_0 \[`matrix(P, P)`\]\cr
#' Prior covariance for class means.
#'
#' @param n_Omega_0 \[`integer(1)`\]\cr
#' Prior degrees of freedom for class covariances.
#'
#' @param V_Omega_0 \[`matrix(P, P)`\]\cr
#' Prior scale matrix for class covariances.
#'
#' @return
#' The functions return one sampler update:
#'
#' - `sample_allocation()`: an integer class label.
#' - `update_s()`: a `C` by 1 numeric matrix of class weights.
#' - `update_z()`: an `N` by 1 numeric matrix of allocations.
#' - `update_m()`: a `C` by 1 numeric matrix of class sizes.
#' - `update_classes_wb()`: a list with `s`, `b`, `Omega`, and `update_type`,
#'   where the latter is zero for no update, one for removal, two for
#'   splitting, and three for merging.
#' - `update_classes_dp()`: a list with `z`, `b`, `Omega`, and `C`.
#'
#' @references
#' \insertRef{Neal2000}{RprobitB}
#'
#' \insertRef{Oelschlaeger2021}{RprobitB}
#'
#' @keywords models
#'
#' @examples
#' ### a latent class state of six deciders with one random coefficient
#' set.seed(1)
#' beta <- matrix(c(-1, -1.2, -0.8, 1, 1.3, 0.9), nrow = 1)
#' b <- matrix(c(-1, 1), nrow = 1)
#' Omega <- matrix(c(0.2, 0.2), nrow = 1)
#'
#' ### the weights, the allocations, and the class sizes are updated in turn
#' s <- update_s(delta = 1, m = c(3, 3))
#' z <- update_z(s, beta, b, Omega)
#' m <- update_m(C = 2, z = z)
#' sample_allocation(c(0.5, 0.3, 0.2))
#'
#' ### the weight-based update splits a class that grew too large
#' update_classes_wb(s = c(0.9, 0.1), b = b, Omega = Omega)
#'
#' ### the Dirichlet process update draws the class count from the data
#' update_classes_dp(
#'   beta = beta, z = z, b = b, Omega = Omega, delta = 1,
#'   mu_b_0 = 0, Sigma_b_0 = diag(1), n_Omega_0 = 4, V_Omega_0 = diag(1)
#' )
#'
#' @name class_updates
#' @rdname class_updates
#' @export
sample_allocation <- function(prob) {
    .Call(`_RprobitB_sample_allocation`, prob)
}

#' @name utility_updates
#' @rdname utility_updates
#' @export
d_to_gamma <- function(d) {
    .Call(`_RprobitB_d_to_gamma`, d)
}

#' @rdname utility_updates
#' @export
log_likelihood_ordered <- function(d, y, sys, Tvec) {
    .Call(`_RprobitB_log_likelihood_ordered`, d, y, sys, Tvec)
}

#' @name coefficient_updates
#' @rdname coefficient_updates
#' @export
update_coefficient <- function(mu_beta_0, Sigma_beta_0_inv, XSigX, XSigU) {
    .Call(`_RprobitB_update_coefficient`, mu_beta_0, Sigma_beta_0_inv, XSigX, XSigU)
}

#' Update coefficient distributions
#'
#' @description
#' Low-level Gibbs sampler kernels for coefficients, their class means, and
#' covariance matrices.
#'
#' @param bar_b_c \[`numeric(P)`\]\cr
#' Average individual coefficient vector in one class.
#'
#' @param Omega_c \[`matrix(P, P)`\]\cr
#' Covariance matrix of one class.
#'
#' @param m_c \[`integer(1)`\]\cr
#' Size of one class.
#'
#' @param Sigma_b_0_inv \[`matrix(P, P)`\]\cr
#' Prior precision matrix for class means.
#'
#' @param mu_b_0 \[`numeric(P)`\]\cr
#' Prior mean for class means.
#'
#' @param beta \[`matrix(P, N)`\]\cr
#' Individual coefficient draws in columns.
#'
#' @param Omega \[`matrix(P * P, C)`\]\cr
#' Vectorized class covariance matrices in columns.
#'
#' @param z \[`numeric(N)`\]\cr
#' Class allocations numbered from one to `C`.
#'
#' @param m \[`numeric(C)`\]\cr
#' Class sizes.
#'
#' @param S_c \[`matrix(P, P)`\]\cr
#' Scatter matrix for one class.
#'
#' @param n_Omega_0 \[`integer(1)`\]\cr
#' Prior degrees of freedom for class covariances.
#'
#' @param V_Omega_0 \[`matrix(P, P)`\]\cr
#' Prior scale matrix for class covariances.
#'
#' @param correlated \[`logical(P)` | `NULL`\]\cr
#' Which random effects are correlated. Covariances between the other
#' effects are zero. In `update_Omega()`, the default (`NULL`) correlates all
#' random effects.
#'
#' @param b \[`matrix(P, C)`\]\cr
#' Class means in columns.
#'
#' @param mu_beta_0 \[`numeric(P)`\]\cr
#' Prior mean for a coefficient vector.
#'
#' @param Sigma_beta_0_inv \[`matrix(P, P)`\]\cr
#' Prior precision matrix for a coefficient vector.
#'
#' @param XSigX \[`matrix(P, P)`\]\cr
#' Sum of design cross-products weighted by inverse error covariance.
#'
#' @param XSigU \[`numeric(P)`\]\cr
#' Sum of design-utility products weighted by inverse error covariance.
#'
#' @return
#' The functions return one sampler update:
#'
#' - `update_b_c()`: a `P` by 1 numeric matrix containing a class mean.
#' - `update_b()`: a `P` by `C` matrix of class means.
#' - `update_Omega_c()`: a `P` by `P` class covariance matrix.
#' - `update_Omega()`: a `P * P` by `C` matrix of vectorized covariances.
#' - `update_coefficient()`: a `P` by 1 numeric coefficient matrix.
#'
#' @keywords models
#'
#' @examples
#' ### four deciders with one random coefficient in two classes
#' set.seed(1)
#' beta <- matrix(c(-1, -1.2, 1, 1.3), nrow = 1)
#' Omega <- matrix(c(0.2, 0.2), nrow = 1)
#' z <- c(1, 1, 2, 2)
#' m <- c(2, 2)
#'
#' ### a coefficient from its conditional posterior
#' update_coefficient(c(0, 0), diag(2), diag(2), c(0, 0))
#'
#' ### the class means, for one class and for all classes at once
#' update_b_c(
#'   bar_b_c = c(0, 0), Omega_c = diag(2), m_c = 4,
#'   Sigma_b_0_inv = diag(2), mu_b_0 = c(0, 0)
#' )
#' update_b(beta, Omega, z, m, Sigma_b_0_inv = diag(1), mu_b_0 = 0)
#'
#' ### the class covariances
#' update_Omega_c(
#'   S_c = diag(2), m_c = 4, n_Omega_0 = 4, V_Omega_0 = diag(2),
#'   correlated = c(TRUE, TRUE)
#' )
#' update_Omega(
#'   beta, b = matrix(c(-1, 1), nrow = 1), z, m,
#'   n_Omega_0 = 4, V_Omega_0 = diag(1)
#' )
#'
#' @rdname coefficient_updates
#' @export
update_b_c <- function(bar_b_c, Omega_c, m_c, Sigma_b_0_inv, mu_b_0) {
    .Call(`_RprobitB_update_b_c`, bar_b_c, Omega_c, m_c, Sigma_b_0_inv, mu_b_0)
}

#' @rdname coefficient_updates
#' @export
update_b <- function(beta, Omega, z, m, Sigma_b_0_inv, mu_b_0) {
    .Call(`_RprobitB_update_b`, beta, Omega, z, m, Sigma_b_0_inv, mu_b_0)
}

#' @rdname coefficient_updates
#' @export
update_Omega_c <- function(S_c, m_c, n_Omega_0, V_Omega_0, correlated) {
    .Call(`_RprobitB_update_Omega_c`, S_c, m_c, n_Omega_0, V_Omega_0, correlated)
}

#' @rdname coefficient_updates
#' @export
update_Omega <- function(beta, b, z, m, n_Omega_0, V_Omega_0, correlated = NULL) {
    .Call(`_RprobitB_update_Omega`, beta, b, z, m, n_Omega_0, V_Omega_0, correlated)
}

#' @rdname class_updates
#' @export
update_s <- function(delta, m) {
    .Call(`_RprobitB_update_s`, delta, m)
}

#' @rdname class_updates
#' @export
update_z <- function(s, beta, b, Omega) {
    .Call(`_RprobitB_update_z`, s, beta, b, Omega)
}

#' @rdname class_updates
#' @export
update_m <- function(C, z, non_zero = FALSE) {
    .Call(`_RprobitB_update_m`, C, z, non_zero)
}

#' @rdname class_updates
#' @export
update_classes_wb <- function(s, b, Omega, epsmin = 0.01, epsmax = 0.7, deltamin = 0.1, deltashift = 0.5, identify_classes = FALSE, Cmax = 10L) {
    .Call(`_RprobitB_update_classes_wb`, s, b, Omega, epsmin, epsmax, deltamin, deltashift, identify_classes, Cmax)
}

#' @rdname class_updates
#' @export
update_classes_dp <- function(beta, z, b, Omega, delta, mu_b_0, Sigma_b_0, n_Omega_0, V_Omega_0, identify_classes = FALSE, Cmax = 10L) {
    .Call(`_RprobitB_update_classes_dp`, beta, z, b, Omega, delta, mu_b_0, Sigma_b_0, n_Omega_0, V_Omega_0, identify_classes, Cmax)
}

#' Update utilities and thresholds
#'
#' @description
#' Low-level Gibbs sampler kernels for error covariance matrices, latent
#' utilities, and ordered-response thresholds.
#'
#' @param n_Sigma_0 \[`integer(1)`\]\cr
#' Prior degrees of freedom for the error covariance.
#'
#' @param V_Sigma_0 \[`matrix(J - 1, J - 1)`\]\cr
#' Prior scale matrix for the error covariance.
#'
#' @param N \[`integer(1)`\]\cr
#' Number of residual vectors summed in the scatter matrix `S`.
#'
#' @param S \[`matrix(J - 1, J - 1)`\]\cr
#' Error scatter matrix.
#'
#' @param U \[`numeric(J - 1)`\]\cr
#' Current latent utility differences.
#'
#' @param y \[`integer(1)` | `matrix(N, max(Tvec))`\]\cr
#' Chosen alternative for `update_U()`, where `J` denotes the base
#' alternative, or ordered responses by decider and occasion for the
#' threshold functions.
#'
#' @param sys \[`numeric(J - 1)` | `matrix(N, max(Tvec))`\]\cr
#' Systematic utilities matching `U` or `y`.
#'
#' @param Sigma_inv \[`matrix(J - 1, J - 1)`\]\cr
#' Inverse error covariance matrix.
#'
#' @param available \[`logical(J)` | `NULL`\]\cr
#' Availability of the `J - 1` non-base alternatives followed by the base
#' alternative. Utilities of unavailable alternatives are drawn without
#' truncation. By default (`NULL`), all alternatives are available.
#'
#' @param d \[`numeric(J - 2)`\]\cr
#' Log-increments between finite ordered thresholds.
#'
#' @param Tvec \[`integer(N)`\]\cr
#' Number of observed occasions for each decider.
#'
#' @param mu_d_0 \[`numeric(J - 2)`\]\cr
#' Prior mean for threshold log-increments.
#'
#' @param Sigma_d_0 \[`matrix(J - 2, J - 2)`\]\cr
#' Prior covariance for threshold log-increments.
#'
#' @param step_scale \[`numeric(J - 2)`\]\cr
#' Random-walk proposal standard deviations, one per log-increment.
#'
#' @return
#' The functions return one sampler update or transformation:
#'
#' - `update_Sigma()`: a `J - 1` by `J - 1` covariance matrix.
#' - `update_U()` and `update_U_ranked()`: `J - 1` by 1 numeric latent
#'   utility matrices.
#' - `d_to_gamma()`: a numeric column matrix containing ordered thresholds
#'   and their infinite bounds.
#' - `log_likelihood_ordered()`: one numeric log-likelihood value.
#' - `update_d()`: a `J - 2` by 1 numeric matrix of updated log-increments.
#'
#' @references
#' \insertRef{Robert1995}{RprobitB}
#'
#' @keywords models
#'
#' @examples
#' ### two deciders who choose twice on an ordered scale of four levels
#' set.seed(1)
#' d <- c(0, log(2))
#' y <- matrix(c(1, 2, 3, 2), nrow = 2)
#' sys <- matrix(0, nrow = 2, ncol = 2)
#' Tvec <- c(2, 2)
#'
#' ### the thresholds, their likelihood, and their random-walk update
#' d_to_gamma(d)
#' log_likelihood_ordered(d, y, sys, Tvec)
#' update_d(
#'   d, y, sys, mu_d_0 = c(0, 0), Sigma_d_0 = diag(2), Tvec = Tvec,
#'   step_scale = c(0.1, 0.1)
#' )
#'
#' ### the latent utilities of an unordered and of a ranked choice
#' update_U(
#'   c(0, 0), y = 1, sys = c(0, 0), Sigma_inv = diag(2),
#'   available = c(TRUE, TRUE, TRUE)
#' )
#' update_U_ranked(c(0, 0), sys = c(0, 0), Sigma_inv = diag(2))
#'
#' ### the error covariance
#' update_Sigma(n_Sigma_0 = 4, V_Sigma_0 = diag(2), N = 10, S = diag(2))
#'
#' @rdname utility_updates
#' @export
update_Sigma <- function(n_Sigma_0, V_Sigma_0, N, S) {
    .Call(`_RprobitB_update_Sigma`, n_Sigma_0, V_Sigma_0, N, S)
}

#' @rdname utility_updates
#' @export
update_U <- function(U, y, sys, Sigma_inv, available = NULL) {
    .Call(`_RprobitB_update_U`, U, y, sys, Sigma_inv, available)
}

#' @rdname utility_updates
#' @export
update_U_ranked <- function(U, sys, Sigma_inv) {
    .Call(`_RprobitB_update_U_ranked`, U, sys, Sigma_inv)
}

#' @rdname utility_updates
#' @export
update_d <- function(d, y, sys, mu_d_0, Sigma_d_0, Tvec, step_scale) {
    .Call(`_RprobitB_update_d`, d, y, sys, mu_d_0, Sigma_d_0, Tvec, step_scale)
}

gibbs_sampler <- function(sufficient_statistics, prior, latent_classes, R, B, ordered, ranked, save_beta_draws, progress = NULL) {
    .Call(`_RprobitB_gibbs_sampler`, sufficient_statistics, prior, latent_classes, R, B, ordered, ranked, save_beta_draws, progress)
}