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

#' Sample allocation
#'
#' @param prob \[`numeric(C)`\]\cr
#' The vector of class probabilities.
#'
#' @return
#' An integer \code{1:C}.
#'
#' @examples
#' sample_allocation(c(0.5, 0.3, 0.2))
#'
#' @export
#'
#' @keywords gibbs_sampler
#'
sample_allocation <- function(prob) {
    .Call(`_RprobitB_sample_allocation`, prob)
}

#' Update class weight vector
#'
#' @param m \[`numeric(C)`\]\cr
#' The vector of current class frequencies.
#'
#' @inheritParams check_prior
#'
#' @return
#' An update for \code{s}.
#'
#' @examples
#' update_s(delta = 1, m = 4:1)
#'
#' @export
#'
#' @keywords gibbs_sampler
#'
update_s <- function(delta, m) {
    .Call(`_RprobitB_update_s`, delta, m)
}

#' Update class allocation vector
#'
#' @inheritParams RprobitB_parameter
#'
#' @return
#' An update for \code{z}.
#'
#' @examples
#' update_z(
#'   s = c(0.6, 0.4), beta = matrix(c(-2, 0, 2), ncol = 3),
#'   b = cbind(0, 1), Omega = cbind(1, 1)
#' )
#'
#' @export
#'
#' @keywords gibbs_sampler
#'
update_z <- function(s, beta, b, Omega) {
    .Call(`_RprobitB_update_z`, s, beta, b, Omega)
}

#' Update class sizes
#'
#' @param non_zero \[`logical(1)`\]\cr
#' Enforce strictly positive values in \code{m} (for numerical stability)?
#'
#' @inheritParams RprobitB_parameter
#'
#' @return
#' An update for \code{m}.
#'
#' @examples
#' update_m(C = 4, z = c(1, 1, 1, 2, 2, 3))
#'
#' @export
#'
#' @keywords gibbs_sampler
#'
update_m <- function(C, z, non_zero = FALSE) {
    .Call(`_RprobitB_update_m`, C, z, non_zero)
}

#' Update mean of a single class
#'
#' @param Omega_c \[`matrix(P_r, P_r)`\]\cr
#' The class covariance matrix.
#'
#' @param bar_b_c \[`numeric(P_r)`\]\cr
#' The average observation of this class.
#'
#' @param m_c \[`integer(1)`\]\cr
#' The number of observations in this class.
#'
#' @param Sigma_b_0_inv \[`matrix(P_r, P_r)`\]\cr
#' The prior precision of the class mean.
#'
#' @inheritParams check_prior
#'
#' @return
#' An update for \code{b_c}.
#'
#' @examples
#' update_b_c(
#'   bar_b_c = c(0, 0), Omega_c = diag(2), m_c = 10,
#'   Sigma_b_0_inv = diag(2), mu_b_0 = c(0, 0)
#' )
#'
#' @export
#'
#' @keywords gibbs_sampler
#'
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)
}

#' Update class means
#'
#' @inheritParams RprobitB_parameter
#' @inheritParams update_s
#' @inheritParams update_b_c
#' @inheritParams check_prior
#'
#' @return
#' A matrix of updated means for each class in columns.
#'
#' @examples
#' N <- 100
#' b <- cbind(c(0, 0), c(1, 1))
#' Omega <- matrix(c(1, 0.3, 0.3, 0.5, 1, -0.3, -0.3, 0.8), ncol = 2)
#' z <- c(rep(1, N / 2), rep(2, N / 2))
#' m <- as.numeric(table(z))
#' beta <- sapply(
#'   z, function(z) oeli::rmvnorm(n = 1, b[, z], matrix(Omega[, z], 2, 2))
#' )
#' update_b(
#'   beta = beta, Omega = Omega, z = z, m = m,
#'   Sigma_b_0_inv = diag(2), mu_b_0 = c(0, 0)
#' )
#'
#' @export
#'
#' @keywords gibbs_sampler
#'
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)
}

#' Update covariance of a single class
#'
#' @param S_c \[`matrix(P_r, P_r)`\]\cr
#' The scatter matrix of this class.
#'
#' @inheritParams update_b_c
#' @inheritParams check_prior
#'
#' @return
#' An update for \code{Omega_c}.
#'
#' @examples
#' update_Omega_c(S_c = diag(2), m_c = 10, n_Omega_0 = 4, V_Omega_0 = diag(2))
#'
#' @export
#'
#' @keywords gibbs_sampler
#'
update_Omega_c <- function(S_c, m_c, n_Omega_0, V_Omega_0) {
    .Call(`_RprobitB_update_Omega_c`, S_c, m_c, n_Omega_0, V_Omega_0)
}

#' Update class covariances
#'
#' @inheritParams RprobitB_parameter
#' @inheritParams update_s
#' @inheritParams check_prior
#'
#' @return
#' A matrix of updated covariance matrices for each class in columns.
#'
#' @examples
#' N <- 100
#' b <- cbind(c(0, 0), c(1, 1))
#' Omega <- matrix(c(1, 0.3, 0.3, 0.5, 1, -0.3, -0.3, 0.8), ncol = 2)
#' z <- c(rep(1, N / 2), rep(2, N / 2))
#' m <- as.numeric(table(z))
#' beta <- sapply(
#'   z, function(z) oeli::rmvnorm(n = 1, b[, z], matrix(Omega[, z], 2, 2))
#' )
#' update_Omega(
#'   beta = beta, b = b, z = z, m = m,
#'   n_Omega_0 = 4, V_Omega_0 = diag(2)
#' )
#'
#' @export
#'
#' @keywords gibbs_sampler
#'
update_Omega <- function(beta, b, z, m, n_Omega_0, V_Omega_0) {
    .Call(`_RprobitB_update_Omega`, beta, b, z, m, n_Omega_0, V_Omega_0)
}

#' Update coefficient vector
#'
#' @param mu_beta_0 \[`numeric(P)`\]\cr
#' The prior mean for the coefficient vector,
#'
#' @param Sigma_beta_0_inv \[`matrix(P, P)`\]\cr
#' The prior precision for the coefficient vector.
#'
#' @param XSigX \[`matrix(P, P)`\]\cr
#' The matrix \eqn{\sum_{n=1}^N X_n'\Sigma^{-1}X_n}.
#'
#' @param XSigU \[`numeric(P)`\]\cr
#' The vector \eqn{\sum_{n=1}^N X_n'\Sigma^{-1}U_n}.
#'
#' @return
#' An update for the coefficient vector.
#'
#' @examples
#' beta_true <- matrix(c(-1, 1), ncol = 1)
#' Sigma <- matrix(c(1, 0.5, 0.2, 0.5, 1, 0.2, 0.2, 0.2, 2), ncol = 3)
#' N <- 100
#' X <- replicate(N, matrix(rnorm(6), ncol = 2), simplify = FALSE)
#' eps <- replicate(
#'   N, oeli::rmvnorm(n = 1, mean = c(0, 0, 0), Sigma = Sigma),
#'   simplify = FALSE
#' )
#' U <- mapply(
#'   function(X, eps) X %*% beta_true + eps, X, eps, SIMPLIFY = FALSE
#' )
#' mu_beta_0 <- c(0, 0)
#' Sigma_beta_0_inv <- diag(2)
#' XSigX <- Reduce(
#'   `+`, lapply(X, function(X) t(X) %*% solve(Sigma) %*% X)
#' )
#' XSigU <- Reduce(
#'   `+`, mapply(function(X, U) t(X) %*% solve(Sigma) %*% U, X, U,
#'   SIMPLIFY = FALSE)
#' )
#' R <- 10
#' beta_draws <- replicate(
#'   R, update_coefficient(mu_beta_0, Sigma_beta_0_inv, XSigX, XSigU),
#'   simplify = TRUE
#' )
#' rowMeans(beta_draws)
#'
#' @export
#'
#' @keywords gibbs_sampler
#'
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 error covariance matrix
#'
#' @param N \[`integer(1)`\]\cr
#' The sample size.
#'
#' @param S \[`matrix(J - 1, J - 1)`\]\cr
#' The sum over the outer products of the residuals
#' \eqn{(\epsilon_n)_{n = 1, \dots, N}}.
#'
#' @inheritParams check_prior
#'
#' @return
#' An update for \code{Sigma}.
#'
#' @examples
#' (Sigma_true <- matrix(c(1, 0.5, 0.2, 0.5, 1, 0.2, 0.2, 0.2, 2), ncol = 3))
#' beta <- matrix(c(-1, 1), ncol = 1)
#' N <- 100
#' X <- replicate(N, matrix(rnorm(6), ncol = 2), simplify = FALSE)
#' eps <- replicate(
#'   N, oeli::rmvnorm(n = 1, mean = c(0, 0, 0), Sigma = Sigma_true),
#'   simplify = FALSE
#' )
#' U <- mapply(function(X, eps) X %*% beta + eps, X, eps, SIMPLIFY = FALSE)
#' n_Sigma_0 <- 4
#' V_Sigma_0 <- diag(3)
#' outer_prod <- function(X, U) (U - X %*% beta) %*% t(U - X %*% beta)
#' S <- Reduce(`+`, mapply(
#'   function(X, U) (U - X %*% beta) %*% t(U - X %*% beta), X, U,
#'   SIMPLIFY = FALSE
#' ))
#' Sigma_draws <- replicate(100, update_Sigma(n_Sigma_0, V_Sigma_0, N, S))
#' apply(Sigma_draws, 1:2, mean)
#'
#' @export
#'
#' @keywords gibbs_sampler
#'
update_Sigma <- function(n_Sigma_0, V_Sigma_0, N, S) {
    .Call(`_RprobitB_update_Sigma`, n_Sigma_0, V_Sigma_0, N, S)
}

#' Update utility vector
#'
#' @param U \[`numeric(J - 1)`\]\cr
#' The current utility vector.
#'
#' @param y \[`integer(1)`\]\cr
#' The index of the chosen alternative, from \code{1} to \code{J}.
#'
#' @param sys \[`numeric(J - 1)`\]\cr
#' The systematic utility.
#'
#' @param Sigma_inv \[`matrix(J - 1, J - 1)`\]\cr
#' The inverted error covariance matrix.
#'
#' @return
#' An update for (a single) \code{U}.
#'
#' @examples
#' U <- sys <- c(0, 0, 0)
#' Sigma_inv <- diag(3)
#' lapply(1:4, function(y) update_U(U, y, sys, Sigma_inv))
#'
#' @export
#'
#' @keywords gibbs_sampler
#'
update_U <- function(U, y, sys, Sigma_inv) {
    .Call(`_RprobitB_update_U`, U, y, sys, Sigma_inv)
}

#' Update ranked utility vector
#'
#' @inheritParams update_U
#'
#' @return
#' An update for (a single) ranked \code{U}.
#'
#' @examples
#' U <- sys <- c(0, 0)
#' Sigma_inv <- diag(2)
#' update_U_ranked(U, sys, Sigma_inv)
#'
#' @export
#'
#' @keywords gibbs_sampler
#'
update_U_ranked <- function(U, sys, Sigma_inv) {
    .Call(`_RprobitB_update_U_ranked`, U, sys, Sigma_inv)
}

#' Transform increments to thresholds
#'
#' @param d \[`numeric(J - 2)`\]\cr
#' Threshold increments.
#'
#' @return
#' The threshold vector of length \code{J + 1}.
#'
#' @examples
#' d_to_gamma(c(0, 0, 0))
#'
#' @export
#'
#' @keywords gibbs_sampler
#'
d_to_gamma <- function(d) {
    .Call(`_RprobitB_d_to_gamma`, d)
}

#' Compute ordered probit log-likelihood
#'
#' @param y \[`matrix(nrow = N, ncol = max(Tvec))`\]\cr
#' Choices \code{1,...,J} for each decider in each choice occasion.
#'
#' @param sys \[`matrix(nrow = N, ncol = max(Tvec))`\]\cr
#' Systematic utilties for each decider in each choice occasion.
#'
#' @param Tvec \[`integer(N)`\]\cr
#' Number of choice occasions per decider.
#'
#' @inheritParams d_to_gamma
#'
#' @return
#' The ordered probit log-likelihood value.
#'
#' @examples
#' d <- c(0, 0, 0)
#' y <- matrix(c(1, 2, 1, NA), ncol = 2)
#' sys <- matrix(c(0, 0, 0, NA), ncol = 2)
#' Tvec <- c(2, 1)
#' ll_ordered(d = d, y = y, sys = sys, Tvec = Tvec)
#'
#' @export
#'
#' @keywords gibbs_sampler
#'
ll_ordered <- function(d, y, sys, Tvec) {
    .Call(`_RprobitB_ll_ordered`, d, y, sys, Tvec)
}

#' Update utility threshold increments
#'
#' @param ll \[`numeric(1)`\]\cr
#' Current log-likelihood value.
#'
#' @param step_scale \[`numeric(1)`\]\cr
#' Scaling the variance for the Gaussian proposal.
#'
#' @inheritParams ll_ordered
#' @inheritParams check_prior
#'
#' @return
#' An update for \code{d}.
#'
#' @examples
#' set.seed(1)
#' N <- 1000
#' d_true <- rnorm(2)
#' gamma <- c(-Inf, 0, cumsum(exp(d_true)), Inf)
#' X <- matrix(rnorm(N * 2L), ncol = 2L)
#' beta <- c(0.8, -0.5)
#' mu <- matrix(as.vector(X %*% beta), ncol = 1L)
#' U <- rnorm(N, mean = mu[, 1], sd = 1)
#' yvec <- as.integer(cut(U, breaks = gamma, labels = FALSE))
#' y <- matrix(yvec, ncol = 1L)
#' Tvec <- rep(1, N)
#' mu_d_0 <- c(0, 0)
#' Sigma_d_0 <- diag(2) * 5
#' d <- rnorm(2)
#' ll <- -Inf
#' R <- 1000
#' for (iter in seq_len(R)) {
#'   ans <- update_d(
#'     d = d, y = y, sys = mu, ll = ll, mu_d_0 = mu_d_0, Sigma_d_0 = Sigma_d_0,
#'     Tvec = Tvec
#'   )
#'   d  <- ans$d
#'   ll <- ans$ll
#' }
#' cbind("true" = d_true, "sample" = d) |> round(2)
#'
#' @export
#'
#' @keywords gibbs_sampler
#'
update_d <- function(d, y, sys, ll, mu_d_0, Sigma_d_0, Tvec, step_scale = 0.1) {
    .Call(`_RprobitB_update_d`, d, y, sys, ll, mu_d_0, Sigma_d_0, Tvec, step_scale)
}

#' Weight-based class updates
#'
#' @param Cmax \[`integer(1)`\]\cr
#' The maximum number of classes, used to allocate space.
#'
#' @param epsmin \[`numeric(1)`\]\cr
#' The threshold weight for removing a class.
#'
#' @param epsmax \[`numeric(1)`\]\cr
#' The threshold weight for splitting a class.
#'
#' @param deltamin \[`numeric(1)`\]\cr
#' The threshold difference in class means for joining two classes.
#'
#' @param deltashift \[`numeric(1)`\]\cr
#' The scale for shifting the class means after a split.
#'
#' @param identify_classes \[`logical(1)`\]\cr
#' Identify classes by decreasing class weights?
#'
#' @inheritParams RprobitB_parameter
#'
#' @details
#' The following updating rules apply:
#'
#' * Class \eqn{c} is removed if \eqn{s_c < \epsilon_{min}}.
#' * Class \eqn{c} is split into two classes, if \eqn{s_c > \epsilon_{max}}.
#' * Two classes \eqn{c_1} and \eqn{c_2} are merged to one class, if
#'   \eqn{||b_{c_1} - b_{c_2}|| < \delta_{min}}.
#'
#' @examples
#' s <- c(0.7, 0.3)
#' b <- matrix(c(1, 1, 1, -1), ncol = 2)
#' Omega <- matrix(c(0.5, 0.3, 0.3, 0.5, 1, -0.1, -0.1, 0.8), ncol = 2)
#'
#' ### no update
#' update_classes_wb(s = s, b = b, Omega = Omega)
#'
#' ### remove class 2
#' update_classes_wb(s = s, b = b, Omega = Omega, epsmin = 0.31)
#'
#' ### split class 1
#' update_classes_wb(s = s, b = b, Omega = Omega, epsmax = 0.69)
#'
#' ### merge classes 1 and 2
#' update_classes_wb(s = s, b = b, Omega = Omega, deltamin = 3)
#'
#' @return
#' A list of updated values for \code{s}, \code{b}, and \code{Omega} and
#' the indicator \code{update_type} which signals the type of class update:
#'
#' - `0`: no update
#' - `1`: removed class
#' - `2`: split class
#' - `3`: merged classes
#'
#' @export
#'
#' @keywords gibbs_sampler
#'
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)
}

#' Dirichlet process class updates
#'
#' @inheritParams RprobitB_parameter
#' @inheritParams check_prior
#' @inheritParams update_classes_wb
#'
#' @return
#' A list of updated values for \code{z}, \code{b}, \code{Omega}, and \code{C}.
#'
#' @examples
#' set.seed(1)
#' z <- c(rep(1, 10),rep(2, 10))
#' b <- matrix(c(5, 5, 5, -5), ncol = 2)
#' Omega <- matrix(c(1, 0.3, 0.3, 0.5, 1, -0.3, -0.3, 0.8), ncol = 2)
#' beta <- sapply(
#'   z, function(z) oeli::rmvnorm(n = 1, b[, z], matrix(Omega[, z], 2, 2))
#' )
#' beta[, 1] <- c(-10, 10)
#' update_classes_dp(
#'   beta = beta, z = z, b = b, Omega = Omega,
#'   delta = 1, mu_b_0 = numeric(2), Sigma_b_0 = diag(2),
#'   n_Omega_0 = 4, V_Omega_0 = diag(2)
#' )
#'
#' @export
#'
#' @keywords gibbs_sampler
#'
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)
}

#' Gibbs sampler for probit models
#'
#' @details
#' This function is not supposed to be called directly, but rather via
#' \code{\link{fit_model}}.
#'
#' @param sufficient_statistics \[`list`\]\cr
#' The output of \code{\link{sufficient_statistics}}.
#'
#' @inheritParams fit_model
#' @inheritParams RprobitB_data
#'
#' @return
#' A list of Gibbs samples for
#' \itemize{
#'   \item \code{Sigma},
#'   \item \code{alpha} (only if \code{P_f > 0}),
#'   \item \code{s}, \code{z}, \code{b}, \code{Omega} (only if \code{P_r > 0}),
#'   \item \code{d} (only if \code{ordered = TRUE}),
#' }
#' and a vector \code{class_sequence} of length \code{R}, where the \code{r}-th
#' entry is the number of classes after iteration \code{r}.
#'
#' @keywords gibbs_sampler
#'
gibbs_sampler <- function(sufficient_statistics, prior, latent_classes, fixed_parameter, R, B, print_progress, ordered, ranked, save_beta_draws = FALSE) {
    .Call(`_RprobitB_gibbs_sampler`, sufficient_statistics, prior, latent_classes, fixed_parameter, R, B, print_progress, ordered, ranked, save_beta_draws)
}