mnp_gibbs: Gibbs sampler for the Bayesian multinomial probit model

In choicer: Discrete Choice Models for Economic Applications

Gibbs sampler for the Bayesian multinomial probit model

Description

Runs the McCulloch-Rossi (1994) Gibbs sampler with Albert-Chib data augmentation for the multinomial probit model in utility differences against a base alternative. The chain operates on the non-identified parameterization (unrestricted Sigma); identified quantities are obtained by normalizing each draw by sigma_11 (handled by run_mnprobit).

Usage

mnp_gibbs(X, y, p, beta_bar, A, nu, V, R, burn, thin, seed, trace = 0L)

Arguments

X	(N*p) x K stacked design matrix of utility differences. Rows are grouped by choice situation, with the p = J - 1 difference rows of situation i ordered by alternative.
y	N vector of choices: 0 for the base alternative, j in 1..p for the j-th non-base alternative.
p	Number of utility differences (J - 1).
beta_bar	K vector, prior mean of beta.
A	K x K prior precision matrix of beta.
nu	Inverse-Wishart prior degrees of freedom (>= p).
V	p x p inverse-Wishart prior scale matrix.
R	Total number of Gibbs iterations.
burn	Number of initial iterations to discard (0 <= burn < R).
thin	Keep every thin-th post-burn-in draw.
seed	Master RNG seed (non-negative; all streams derive from it).
trace	Print progress every trace iterations (0 = silent).

Details

The latent-utility sweep is parallelized with OpenMP across choice situations (they are conditionally independent given beta and Sigma). Each (iteration, observation) pair uses its own RNG stream, so draws are reproducible given the seed and a fixed thread count; across different thread counts they are invariant only up to floating-point reduction-order round-off (~1e-15), not bitwise (see set_num_threads()).

Value

List with betadraw (R_keep x K), sigmadraw (R_keep x p(p+1)/2, lower triangle of Sigma in row-major order), and R_keep.

Examples

library(data.table)
set.seed(42)
N <- 100; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
d <- prepare_mnp_data(dt, "id", "alt", "choice", c("x1", "x2"))
out <- choicer:::mnp_gibbs(d$X, d$y, d$p,
  beta_bar = rep(0, d$K), A = 0.01 * diag(d$K),
  nu = d$p + 3, V = (d$p + 3) * diag(d$p),
  R = 500, burn = 100, thin = 1, seed = 42)
colMeans(out$betadraw)

choicer documentation built on Sept. 5, 2026, 1:07 a.m.