The math behind choicer: hierarchical Bayesian multinomial probit

This article gives a self-contained description of the hierarchical Bayesian multinomial probit with iid utility shocks (HMNP) implemented by run_hmnprobit() and src/hmnprobit.cpp. Shared hierarchy code is in src/hb_internal.h and sampling primitives are in src/bayes_samplers.h and src/rng.h.

Model boundary. HMNP is not a hierarchical version of the unrestricted-covariance MNP fitted by run_mnprobit(). HMNP imposes iid utility-level normal shocks with a common expanded variance. Its substitution flexibility comes from persistent taste heterogeneity $W$ and pooled alternative effects. The standalone MNP instead estimates a full covariance matrix of utility differences but has no respondent hierarchy.

1. Model and data layout

Respondents are $i=1,\ldots,N$, tasks are $t=1,\ldots,T_i$, and $C_{it}\subseteq\{1,\ldots,J\}$ is the possibly unbalanced set of inside alternatives. The implicit outside option is $o$. Structural covariates are $x_{ijt}\in\mathbb R^K$ and alternative characteristics are $z_j\in\mathbb R^P$. Utilities on the sampler’s raw scale are

$$U_{ijt}=x_{ijt}'\beta_i+\delta_j+\varepsilon_{ijt},\qquad U_{iot}=\varepsilon_{iot},$$

$$\varepsilon_{ijt},\varepsilon_{iot}\stackrel{iid}{\sim}N(0,\sigma^2), \qquad y_{it}=\arg\max_{a\in C_{it}\cup\{o\}}U_{iat}.$$

The outside utility is stochastic, not a deterministic zero threshold. The sampler retains one latent outside utility per task. It operates in undifferenced utility space, so each augmented error is independent and unbalanced choice sets require no special covariance submatrix.

Taste and alternative hierarchies are

$$\beta_i\mid b,W\sim N_K(b,W),$$

$$\delta_j=z_j'\theta+\xi_j,\qquad \xi_j\mid\sigma_d^2\stackrel{ind}{\sim}N(0,\sigma_d^2).$$

All taste coordinates are normal. HMNP does not offer rc_dist: exponentiating a coordinate would make the $\beta_i$ regression conditional non-Gaussian and destroy the fully conjugate update used by the kernel.

2. Exact relationship to an MNP in differences

Difference each inside utility against the outside:

$$w_{ijt}=U_{ijt}-U_{iot}=x_{ijt}'\beta_i+\delta_j+\eta_{ijt},\qquad \eta_{ijt}=\varepsilon_{ijt}-\varepsilon_{iot}.$$

Within a task,

$$\operatorname{Var}(\eta_{ijt})=2\sigma^2,\qquad \operatorname{Cov}(\eta_{ijt},\eta_{ikt})=\sigma^2\quad(j\ne k).$$

Thus the difference-error covariance for a task with $m$ inside alternatives is

$$\Sigma_m=\sigma^2(I_m+\mathbf1_m\mathbf1_m'),$$

with pairwise correlation $1/2$. This is a fixed equicorrelation structure, not an estimated unrestricted covariance. On HMNP’s reported normalization $\sigma=1$, the diagonal is 2 and the off-diagonal is 1. On run_mnprobit()’s convention of dividing by the first diagonal element, the corresponding matrix would have diagonal 1 and off-diagonal $1/2$. This distinction is essential when comparing the two classes.

3. Identification and per-draw normalization

The observational choice rule is unchanged by the positive scaling map

$$(U,\beta_i,b,\delta,\theta,\sigma, W,\sigma_d^2) \mapsto (cU,c\beta_i,cb,c\delta,c\theta,c\sigma,c^2W,c^2\sigma_d^2).$$

The likelihood therefore has no information about one common utility scale. The chain samples a proper raw-scale posterior with free $\sigma^2$ and raw-scale priors. Every retained draw is then mapped to the identified convention $\tilde\sigma=1$:

$$\tilde b=\frac b\sigma,\quad \tilde\beta_i=\frac{\beta_i}\sigma,\quad \tilde W=\frac W{\sigma^2},\quad \tilde\delta=\frac\delta\sigma,\quad \tilde\theta=\frac\theta\sigma,\quad \tilde\sigma_d^2=\frac{\sigma_d^2}{\sigma^2}.$$

Normalization is applied draw by draw. For example, $E(b/\sigma)\ne E(b)/E(\sigma)$ in general. Top-level coefficient, covariance, alternative-effect, respondent-effect, and diagnostic summaries all use these normalized draws. The raw $\sigma^2$ chain is retained as an expansion diagnostic but is deliberately excluded from the automatic convergence warning because it is not an identified estimand.

The proper raw-scale priors are not invariant under this map: they induce a joint prior over all normalized quantities. Consequently the expanded variance prior is not ignorable. Prior-predictive analysis and sensitivity to $(a_0,s_0)$, as well as to the $W$ and $\sigma_d$ priors, are part of serious use.

The outside option fixes utility location: its systematic utility is zero, so there is no base inside alternative and no sum-to-zero constraint on $\delta$. v0.2.0 requires the outside option.

For $J=1$, $U_1-U_o$ has variance $2\sigma^2$. With reported $\sigma=1$, choicer’s coefficient is therefore $\sqrt2$ times the coefficient from glm(..., binomial(link="probit")), whose latent error variance is normalized to one. The test suite uses this binary identity as a validation check.

4. Priors and augmented posterior

The priors shared with HMNL are

$$b\sim N_K(\bar b,A^{-1}),\qquad W\sim IW_K(\nu,V),\qquad \theta\sim N_P(\bar\theta,A_\theta^{-1}).$$

The inverse-Wishart convention is

$$p(W)\propto |W|^{-(\nu+K+1)/2} \exp\{-\tfrac12\operatorname{tr}(VW^{-1})\}.$$

Defaults are $\bar b=0$, $A=.01I_K$, $\nu=K+3$, $V=\nu I_K$, $\bar\theta=0$, and $A_\theta=.01I_P$. The raw prior mean of $W$, when it exists, is $V/(\nu-K-1)=(K+3)I_K/2$.

The default alternative-effect prior is $\sigma_d\sim C^+(0,s_d)$ through

$$\sigma_d^2\mid a_d\sim IG(1/2,1/a_d),\qquad a_d\sim IG(1/2,1/s_d^2),$$

or optionally $\sigma_d^2\sim IG(c_0,d_0)$. HMNP adds

$$\sigma^2\sim IG(a_0,s_0),$$

with defaults $a_0=s_0=3$. Throughout, $IG(a,d)$ is shape–scale: $1/X\sim\operatorname{Gamma}(a,\text{rate}=d)$.

Let $\mathcal B(y)$ be the set of latent utilities in which each observed choice exceeds every other latent in its task. The augmented posterior is proportional to

$$\begin{aligned} &\mathbf1\{U\in\mathcal B(y)\} (\sigma^2)^{-n_L/2} \exp\left[-\frac1{2\sigma^2}\left\{ \sum_{it}\sum_{j\in C_{it}}(U_{ijt}-x_{ijt}'\beta_i-\delta_j)^2 +\sum_{it}U_{iot}^2\right\}\right]\\ &\quad\times\prod_i\phi_K(\beta_i;b,W) \prod_j\phi(\delta_j;z_j'\theta,\sigma_d^2) p(b)p(W)p(\theta)p(\sigma_d^2,a_d)p(\sigma^2), \end{aligned}$$

where $n_L=\sum_{it}(|C_{it}|+1)$ counts every inside latent and one outside latent per task.

5. One fully conjugate Gibbs iteration

5.1 Latent utilities

Within each task the kernel scans inside rows in stored order and then the outside latent. Because shocks are independent, an untruncated inside conditional is

$$U_{ijt}\mid\cdot\sim N(x_{ijt}'\beta_i+\delta_j,\sigma^2),$$

and the outside conditional is $N(0,\sigma^2)$. If alternative $a$ is chosen, its current update is truncated below by the maximum current latent among all $b\ne a$; each nonchosen update is truncated above by the current latent of $a$. The bounds use already-updated values for preceding coordinates and old values for later coordinates, making this an ordinary scalar Gibbs scan of the truncated joint normal. Ties have probability zero.

5.2 Respondent coefficients

Let $X_i$ stack respondent $i$’s inside rows, $U_i$ their latent utilities, and $d_i$ the vector repeating the relevant $\delta_j$. Define the precomputed Gram matrix $G_i=X_i'X_i$. Then

$$Q_i=W^{-1}+G_i/\sigma^2,\qquad r_i=W^{-1}b+X_i'(U_i-d_i)/\sigma^2,$$

$$\beta_i\mid\cdot\sim N_K(Q_i^{-1}r_i,Q_i^{-1}).$$

After each draw the row cache $x_{ijt}'\beta_i$ is refreshed. The outside latents do not enter this regression because their design vectors are zero.

5.3 Alternative effects

Given augmented utilities, different $\delta_j$ appear in disjoint normal regressions and are conditionally independent. Let $\mathcal R_j$ be all inside rows for alternative $j$, $n_j=|\mathcal R_j|$, and $m_j=z_j'\theta$. Then

$$q_j=n_j/\sigma^2+1/\sigma_d^2,\qquad r_j=\frac1{\sigma^2}\sum_{r\in\mathcal R_j}(U_r-x_r'\beta_{i(r)}) +\frac{m_j}{\sigma_d^2},$$

$$\delta_j\mid\cdot\sim N(r_j/q_j,q_j^{-1}).$$

These $J$ draws can be work-shared safely. This is the key computational contrast with HMNL’s serial, coupled softmax sweep.

5.4 Hierarchy and shock variance

The shared hierarchy conditionals are

$$b\mid\cdot\sim N(P_b^{-1}r_b,P_b^{-1}),\quad P_b=A+NW^{-1},\quad r_b=A\bar b+W^{-1}\sum_i\beta_i,$$

$$W\mid\cdot\sim IW_K\left(\nu+N,V+\sum_i(\beta_i-b)(\beta_i-b)'\right),$$

$$\theta\mid\cdot\sim N(P_\theta^{-1}r_\theta,P_\theta^{-1}),\quad P_\theta=A_\theta+Z'Z/\sigma_d^2,\quad r_\theta=A_\theta\bar\theta+Z'\delta/\sigma_d^2.$$

With $\xi=\delta-Z\theta$, the half-Cauchy blocks are

$$\sigma_d^2\mid\cdot\sim IG((J+1)/2,1/a_d+\xi'\xi/2), \qquad a_d\mid\cdot\sim IG(1,1/s_d^2+1/\sigma_d^2).$$

After the new $\beta$ and $\delta$ draws, a separate fixed-order pass computes

$$RSS=\sum_{it}\sum_{j\in C_{it}}(U_{ijt}-x_{ijt}'\beta_i-\delta_j)^2 +\sum_{it}U_{iot}^2.$$

Then

$$\sigma^2\mid\cdot\sim IG(a_0+n_L/2,s_0+RSS/2).$$

Using a separate RSS pass matters: reusing residuals from the coefficient phase would condition on the stale $\delta$.

6. Initialization, storage, and diagnostics

The kernel starts $\beta_i=0$, $b=0$, $\delta=0$, $\theta=0$, $W=I_K$, and $\sigma_d^2=a_d=\sigma^2=1$, regardless of nonzero prior means. Chosen latents start at $+1$ and nonchosen latents at $-1$, including the outside latent. Burn-in, rather than an optimizer, moves the chain away from this feasible state.

The C++ kernel records raw draws of $b$, vech$(W)$, $\delta$, $\theta$, $\sigma_d^2$, and $\sigma^2$. At recording time it divides $\beta_i$, $\delta$, and $\xi$ by the current $\sigma$ before Welford updates or cube storage. The R wrapper applies the same per-draw normalization to all hierarchical chains and constructs object$draws; raw chain-1 versions are retained with _raw names. object$chains contains normalized $b,W,\delta,\theta,\sigma_d^2$ plus raw $\sigma^2$ for every requested chain. Chains run sequentially; their seeds are offset by one.

The same rank-normalized folded split-$\widehat R$, bulk/tail ESS, MCSE, trace, and worst-$\delta$ conventions described for HMNL apply. summary() may display the raw $\sigma^2$ diagnostic for transparency, but the automatic convergence warning excludes it because movement along the expansion scale is not convergence of an identified estimand. vech(W) is retained but also omitted from the automatic table; it should be checked explicitly across object$chains, particularly in cross-sectional data. Normalized parameters, every $\delta_j$, prior sensitivity, and posterior predictive behavior are the relevant checks.

keep_beta_i="draws" is guarded using approximately $1.9\times8KNR_{keep}$ bytes per chain times the number of sequential chains, because all returned C++ cubes coexist before the final object is assembled. The fit object exposes chain-1 respondent summaries/draws; hierarchical draws from every chain remain in object$chains.

7. Choice probabilities and Gauss–Hermite quadrature

On the identified scale $\sigma=1$, condition on the shock $u$ of candidate $j$. With systematic utilities $V_a$ and $V_o=0$,

$$P_j=\int_{-\infty}^{\infty}\phi(u) \prod_{k\ne j}\Phi(V_j-V_k+u)\,du.$$

The product includes the outside option when $j$ is inside and every inside alternative when $j=o$. predict.choicer_hb() evaluates this one-dimensional integral separately for every candidate using 20-node physicists’ Gauss–Hermite quadrature. If $(t_m,w_m)$ integrate $e^{-t^2}$, then

$$P_j\approx\sum_{m=1}^{20}\frac{w_m}{\sqrt\pi} \prod_{k\ne j}\Phi(V_j-V_k+\sqrt2t_m).$$

The resulting candidate probabilities are divided by their numerical sum to remove quadrature drift. This final renormalization is numerical, not part of the probit model.

As with HMNL, population prediction pairs posterior draw $r$ with one $\beta^{(r)}\sim N(b^{(r)},W^{(r)})$; it is a Monte Carlo average over heterogeneity. Individual prediction uses stored $\beta_i^{(r)}$ when available and otherwise plugs in the posterior mean. New alternatives receive $\delta_a^{(r)}\sim N(z_a'\theta^{(r)},\sigma_d^{2(r)})$. Post-estimation simulations use R’s RNG.

Finite-difference elasticities() and diversion_ratios() use these predicted aggregate shares and common random-number paths, exactly as for HMNL. wtp() uses drawwise population-mean coefficient ratios; all HMNP coefficients are normal, so those means are the components of $b^{(r)}$.

The EV1 logsum identity does not apply to normal shocks. Although $E[\max_a U_a]$ is mathematically defined and can be evaluated numerically or by simulation, v0.2.0 does not implement that calculation. logsum() and consumer_surplus() therefore error for HMNP rather than applying an incorrect logit formula.

8. Computational contract

One persistent OpenMP region contains the whole chain. A work-shared respondent pass performs latent updates, $\beta_i$ draws, and cache refreshes; a work-shared alternative pass draws $\delta_j$; another respondent pass computes RSS. Hierarchy draws, fixed-order reductions, recording, and interrupt checks are master-only. Worker paths use hand-written dense Cholesky and triangular solvers and never call R’s API or BLAS/LAPACK.

At iteration $r$, RNG tags $i$ drive respondent latent sweeps, $N+i$ drive respondent coefficient draws, $2N+j$ drive alternative effects, and $2N+J+\{0,1,2,3\}$ drive $b,W,\theta$, and then $(\sigma_d^2,\sigma^2)$ in a fixed order from the last shared stream. This partition is collision-free even when $J>N$. Results are deterministic for a fixed seed and execution environment; byte identity across thread counts or toolchains is not an API guarantee.

References

• Albert, J. H., & Chib, S. (1993). Bayesian analysis of binary and polychotomous response data. Journal of the American Statistical Association, 88(422), 669–679.
• Gelman, A. (2006). Prior distributions for variance parameters in hierarchical models. Bayesian Analysis, 1(3), 515–534.
• Liu, J. S., & Wu, Y. N. (1999). Parameter expansion for data augmentation. Journal of the American Statistical Association, 94(448), 1264–1274.
• Makalic, E., & Schmidt, D. F. (2016). A simple sampler for the horseshoe estimator. IEEE Signal Processing Letters, 23(1), 179–182.
• McCulloch, R. E., & Rossi, P. E. (1994). An exact likelihood analysis of the multinomial probit model. Journal of Econometrics, 64(1–2), 207–240.
• Rossi, P. E., Allenby, G. M., & McCulloch, R. (2005). Bayesian Statistics and Marketing. Wiley.
• Train, K. E. (2009). Discrete Choice Methods with Simulation (2nd ed.). Cambridge University Press.
• Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., & Bürkner, P.-C. (2021). Rank-normalization, folding, and localization: an improved $\widehat R$ for assessing convergence of MCMC. Bayesian Analysis, 16(2), 667–718.