The math behind choicer: mixed logit

This document provides a detailed mathematical description of the mixed logit (random coefficients logit) model as implemented in src/mxlogit.cpp.

Table of Contents

0. Notation
1. Model Definition
2. Log-Likelihood Function
3. Gradient Computation
4. Hessian Computation
5. Implementation Details
6. Elasticity Computation
7. Diversion Ratios
8. BLP Contraction Mapping
9. Choice-Based Sampling and WESML Weighting
10. Willingness to Pay
11. Consumer Surplus and the Simulated Logsum
12. Goodness of Fit
13. Conditional (Individual-Level) Tastes

---

Notation

Symbol	Description
$i = 1, \ldots, N$	Index for choice situations
$j = 1, \ldots, J_i$	Index for alternatives available to individual $i$
$s = 1, \ldots, S$	Index for simulation draws
$j_i$	The alternative chosen by individual $i$
$w_i$	Weight for individual $i$
$X_{ij}$	Row vector of covariates with fixed coefficients ($1 \times K_x$)
$W_{ij}$	Row vector of covariates with random coefficients ($1 \times K_w$)
$\beta$	Fixed coefficient vector ($K_x \times 1$)
$\mu$	Mean parameters for random coefficients ($K_w \times 1$)
$L$	Lower-triangular Cholesky factor ($K_w \times K_w$)
$\Sigma = LL'$	Covariance matrix of the latent Gaussian deviations (and of realized coefficients only in the all-normal case)
$\delta_j$	Alternative-specific constant (ASC) for alternative $j$
$\eta_i^s$	Standard-normal integration draw for choice situation $i$ and draw $s$ ($K_w \times 1$)
$\gamma_i^s = L\eta_i^s$	Latent Gaussian random-coefficient deviation
$n = 1, \ldots, U$	Index for decision makers (respondents) when person_col groups situations into a panel; $U$ is the number of decision makers
$t = 1, \ldots, T_n$	Index for the choice situations belonging to decision maker $n$; $T_n$ is their count
$\omega_{ns}$	Posterior weight of draw $s$ given decision maker $n$’s observed choices ($$2.6)

Sections 2–12 are written per choice situation $i$, the default likelihood unit. §1.6 introduces the decision maker $n$ as an optional coarser unit; §2.6, §3.5 and §4.5 give the corresponding panel likelihood, score and Hessian. Without person_col, $U = N$ and $T_n \equiv 1$ for every $n$, so $i$ and $(n, t{=}1)$ index the same choice situation and every panel formula collapses to its §2–4 counterpart.

---

1. Model Definition

1.1 Utility Specification

The utility that individual $i$ derives from alternative $j$ under simulation draw $s$ is:

$$V_{ij}^s = X_{ij}\beta + W_{ij}\left(\mu^* + \gamma_i^{s*}\right) + \delta_j$$

where: - $X_{ij}\beta$ captures the systematic utility from fixed coefficients - $W_{ij}\mu^*$ captures the mean contribution of random coefficients - $W_{ij}\gamma_i^{s*}$ captures the individual-specific random deviation - $\delta_j$ is the alternative-specific constant (the first inside ASC is normalized without an outside option; with an outside option all inside ASCs are free)

1.2 Random Coefficient Structure

Random coefficients are generated through a Cholesky decomposition structure:

$$\gamma_i^s = L \eta_i^s, \quad \text{where } \eta_i^s \sim N(0, I_{K_w})$$

This implies:

$$\gamma_i^s \sim N(0, \Sigma), \quad \text{where } \Sigma = LL'$$

For finite parameter values the exponentiated diagonal makes $L$ nonsingular, so $\Sigma$ is positive definite (positive semidefinite in the limiting sense as a diagonal element tends to zero). Here $\Sigma$ is the covariance of the latent Gaussian deviation $\gamma$. It is the covariance of the realized coefficient vector only when all entries are normal. After any componentwise log-normal transformation, the realized coefficient covariance is a nonlinear function of $L$ and is not $LL'$.

1.3 Distribution Transformations

The implementation supports two distributions for each random coefficient $k$:

Normal Distribution (rc_dist[k] = 0): $$\mu_k^* = \mu_k, \quad \gamma_{ik}^{s*} = \gamma_{ik}^s$$

Log-Normal Distribution (rc_dist[k] = 1): $$\mu_k^* = \exp(\mu_k), \quad \gamma_{ik}^{s*} = \exp(\gamma_{ik}^s)$$

Note: This is a shifted log-normal parameterization: $\beta_k = \exp(\mu_k) + \exp(L_k \eta)$, which differs from the textbook parameterization $\beta_k = \exp(\mu_k + L_k \eta)$.

When rc_mean = FALSE, $\mu^*=0$ for every distribution. Thus a normal random coefficient is centered at zero, while a log-normal random coefficient is $\exp((L\eta)_k)$ and has median one. Correlation is introduced on the latent-normal scale by the rows of $L$; componentwise transformations then induce a generally non-Gaussian joint coefficient distribution.

1.4 Cholesky Parameterization

For a full covariance matrix (rc_correlation = TRUE), the lower-triangular matrix $L$ is parameterized as:

$$L = \begin{pmatrix} \exp(\ell_{11}) & 0 & 0 & \cdots \\ \ell_{21} & \exp(\ell_{22}) & 0 & \cdots \\ \ell_{31} & \ell_{32} & \exp(\ell_{33}) & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{pmatrix}$$

where: - Diagonal elements: $L_{pp} = \exp(\ell_{pp})$ (ensures positivity) - Off-diagonal elements: $L_{pq} = \ell_{pq}$ for $p > q$ (unconstrained)

For diagonal-only covariance (rc_correlation = FALSE):

$$L = \text{diag}\left(\exp(\ell_{11}), \exp(\ell_{22}), \ldots, \exp(\ell_{K_w K_w})\right)$$

The parameter vector for $L$ has length: - Full: $K_w(K_w + 1)/2$ - Diagonal: $K_w$

1.5 Scaling the Optimizer’s Coordinates

Motivation. scale_vars changes how the optimizer moves, never what it maximizes. L-BFGS starts each iteration from a single scalar multiple of the identity as its curvature, so it is slow, and can fail, when the log-likelihood’s curvature differs by orders of magnitude across parameters. Two things spread it. Units: a price in thousands next to a 0/1 dummy gives their coefficients curvatures far apart. The model: the constant of a small alternative (a school with few applicants, a rarely chosen hospital) has far less curvature than a large one’s, and the Cholesky entries of $L$ vary as well. Column scales reach the first; only the curvature itself reaches the second.

A change of coordinates. The kernels evaluate the model at $\theta$, on the data as supplied. The optimizer works on $\tilde\theta$, with

$$\theta = m \circ \tilde\theta + c, \qquad m > 0 \text{ elementwise}.$$

The objective in its coordinates is $\tilde f(\tilde\theta) = f(m \circ \tilde\theta + c)$, with gradient $\nabla\tilde f = m \circ \nabla f$: the map is affine and diagonal, $J = \mathrm{diag}(m)$, so the chain rule has no second-order term and the Hessian in $\tilde\theta$ is $J H J$. The maximizer is the same: $\hat\theta = m \circ \hat{\tilde\theta} + c$. Start values and bounds are given in natural units and mapped by $(\theta - c) / m$, which keeps their order. scale_vars = "bhhh" is the default of run_mxlogit() and run_mnlogit(); "none" bypasses the map entirely.

Column scales ("sd", "mad", "iqr"). Let $s_k^X$ and $s_k^W$ be the columns’ scales. The map

Parameter	$m$	$c$
$\beta_k$	$1 / s_k^X$	0
$\mu_k$ (normal RC)	$1 / s_k^W$	0
$\ell_{pp}$ (Cholesky diagonal, stored as $\log L_{pp}$)	1	$-\log s_p^W$
$\ell_{pq}$ (Cholesky off-diagonal, $p > q$)	$1 / s_p^W$	0
$\text{ASC}_j$	1	0

makes $\tilde f$ the log-likelihood of the design with every column divided by its scale: $(x / s)(s \beta) = x \beta$, and $\tilde W_p = W_p / s_p^W$ implies $L^{\mathrm{scaled}}_{p\cdot} = s_p^W L_{p\cdot}$, row $p$ of $L$ (on the log diagonal, an additive shift; the scale is keyed on the row $p$, not the column $q$). The data are never divided: the optimizer sees the scaled problem through the map. "sd" uses the sample standard deviation; "mad" (stats::mad, $1.4826 \times \text{median}|X_{\cdot k} - \text{median}(X_{\cdot k})|$) and "iqr" (stats::IQR(x) / 1.349) are SD-equivalent under normality and preferable when outliers dominate a column’s SD.

Why centering is omitted. Without an outside option, subtracting a column’s mean from every alternative shifts every utility in a situation by the same amount, which cancels from the softmax; the same argument applies draw by draw to $W$. With an implicit outside option fixed at utility zero, an inside-only shift does not cancel against the outside: a fixed-$X$ shift could be offset only by shifting every inside ASC, and a random-$W$ shift cannot generally be absorbed by fixed ASCs. So the column-scale map corresponds to dividing the data by their scales, never to centering them.

Log-normal coefficients. For a log-normal coefficient ($\beta_k = \exp(\mu_k) + \exp((L\eta)_k)$, rc_dist[k] == 1), dividing the column by $s_k^W$ has no closed form in $(\mu_k, L_{k\cdot})$, because the two exponentials cannot share a scale factor; the column-scale map leaves those parameters at $m = 1$, $c = 0$.

BHHH scales ("bhhh"). With $B_0 = \sum_u w_u s_u s_u^\top$ the BHHH (outer product of gradients) matrix at the start values $\theta_0$, from the unit scores $s_u$ and the fit’s weights,

$$m_k = 2^{\mathrm{round}(\log_2 (1 / \sqrt{B_{0,kk}}))}, \qquad c = 0,$$

for every parameter, the constants, the Cholesky entries and log-normal coefficients included. The diagonal of $J B_0 J$ is then within a factor of two of one, for entries neither floored nor clamped: insofar as $B_0$ tracks the curvature, the optimizer starts from curvatures near one in those coordinates. One extra gradient pass returns $\mathrm{diag}(B_0)$, accumulated beside the gradient (the full matrix is never formed). A $B_{0,kk}$ that is not finite, or below $10^{-12}$ times the median of the finite positive ones, carries no information at $\theta_0$ and takes the median’s scale; scales are clamped to $[2^{-400}, 2^{400}]$. Powers of two make the map exact: $\theta \to \tilde\theta \to \theta$ round-trips without rounding, bounds included. The scales are fixed for the whole run, because L-BFGS’s line search and stored curvature pairs assume one fixed function of $\tilde\theta$. A start near the estimate (theta_init) gives $B_0$ closer to the information at the optimum. The scales fix the optimizer’s metric, not a start that is wrong in a column’s units; the default start is set in the columns’ units (below).

The default start. With theta_init = NULL the start is zero for $\beta$, $\mu$, the off-diagonal $\ell_{pq}$ and the constants, and on the diagonal

$$\ell_{pp} = -\log h_p, \qquad L_{pp} = 1 / h_p,$$

for a normal random coefficient, where $h_p$ is the typical step of column $p$: the distance between its values for a column with two (a dummy, a two-level attribute), the mean of $|W_p - v|$ over the other rows when one value $v$ fills more than half of the column (a mostly-zero column), and its standard deviation otherwise. At the start, a step of $h_p$ in $W_p$ moves the utility by a random amount with standard deviation one (the logit error’s is $\pi / \sqrt{6} \approx 1.28$). For a dummy with a share $\phi$ of ones the standard deviation, $\sqrt{\phi(1-\phi)}$, would understate the step its coefficient multiplies, which is 1; with one of $J$ alternatives marked per situation, $\phi \approx 1/J$, it would start $L_{pp}$ about $\sqrt{J}$ times higher. For a column far from zero with an implicit outside option, the contrast with the outside good is the level $W_p$ itself, so a root mean square might seem the step; on a census-style design with a price of mean 100 and standard deviation 5 it put the start at the true heterogeneity and took 169 evaluations, against 104 from the standard deviation (the level errs high, below). A column with no variation (a step below $10^{-12}$ of its level) or a step that is not finite starts at $L_{pp} = 1$, and steps are clamped to $[2^{-500}, 2^{500}]$. A default start outside lower/upper is moved to the nearest bound, with a message.

Rescaling the column, $W_p \to a W_p$, gives the same model with row $p$ of $L$ (and $\mu_p$, with rc_mean = TRUE) divided by $a$, its $\ell_{pp}$ shifted by $-\log a$; $h_p$ is multiplied by $a$, so this start transforms with it. So does the map of "sd" (and of "mad" and "iqr"), whose shift $-\log s_p^W$ puts $\tilde\ell_{pp}$ at $\log(s_p^W / h_p)$, a constant that does not depend on the units ($0$ under "sd" for a column whose step is its standard deviation): in exact arithmetic the optimizer’s problem is then the same in any units. Under "bhhh", a parameter divided by $a$ has its score multiplied by $a$, so $B_{0,kk}$ is multiplied by $a^2$ and $m_k$ (before rounding, for entries neither floored nor clamped) divided by $a$: $\tilde\theta_k = \theta_k / m_k$ is unchanged up to the power-of-two rounding. $\ell_{pp}$, shifted by $-\log a$, keeps its score and its scale, so $\tilde\ell_{pp}$ is translated by $-\log a / m_{pp}$, which L-BFGS’s steps do not see; only nloptr’s relative step test (xtol_rel, $10^{-8}$ by default) does.

From $L_{pp} = 0.5$ in each column’s own units (the 0.2.x default), no choice of scale_vars reached the maximum of a census-style design with columns rescaled by $10^3$ or $10^{-3}$ in 2,000 likelihood evaluations; from this start "bhhh" reaches it in 54, against 56 on the original units. The level, one, errs high on purpose. $B_0$ is fixed at the start, and while $L_{pp}$ is small the simulated likelihood depends on it only at second order (the draws have mean zero, or nearly): the score in $\ell_{pp}$ is $O(L_{pp}^2)$, so $B_{0,pp}$ is $O(L_{pp}^4)$ and $m_{pp} \propto L_{pp}^{-2}$ at the start, while the curvature in $\ell_{pp}$ grows like $L_{pp}^2$ as the optimizer raises it. A start with too little heterogeneity therefore yields a scale too generous for where the optimizer goes: from $0.25 / \mathrm{sd}(W_p)$ one design reported convergence 554 log-likelihood units short.

A log-normal coefficient starts at $\ell_{pp} = \log 0.5$: there $L_{pp}$ is the standard deviation of the logarithm of the coefficient’s random part, which has no units. That random part, $\exp((L\eta)_p)$, has median one in any units, so no start makes a log-normal column’s units irrelevant: its covariate belongs in units where the coefficient is of order one.

Variance. The variance is computed at $\hat\theta$, in natural units, as for an unscaled fit: $V = A^{-1}$ with $A = -H(\hat\theta)$ the information (the negated Hessian of the log-likelihood), the BHHH variance $\big(\sum_u w_u s_u s_u^\top\big)^{-1}$ at $\hat\theta$, or the sandwich $A^{-1} B A^{-1}$ of §9.3 (with $B$ the meat, its weights squared). A scaled fit’s natural information can be ill-conditioned for the reason it was scaled, so the matrix $M$ being inverted ($A$, or the BHHH matrix) is equilibrated:

$$V = D \, (D M D)^{-1} D, \qquad D = \mathrm{diag}\big(2^{\mathrm{round}(-\log_2 |M_{kk}| / 2)}\big),$$

and for the sandwich $V = D \, [(DAD)^{-1} (DBD) (DAD)^{-1}] \, D$, with $D$ from $A$. The products with $D$ are exact (powers of two), so only the inversion rounds (for the sandwich, also the products with the meat and their symmetrization), on a matrix whose diagonal is within a factor of two of one in magnitude. On one test design the natural bread had condition number $2.5 \times 10^{13}$ and the equilibrated one $7.4$; inverted raw, the robust variance was off by $1.6 \times 10^{-8}$ relative, equilibrated by $8 \times 10^{-16}$. The magnitude of the diagonal is used so that an information that is not positive definite away from the optimum is equilibrated too, and its negative variances reported. Post-hoc vcov(fit, type = ) uses the same route; an unscaled fit keeps the raw inversion.

The other models. run_mnlogit() and run_nestlogit() take the same map: column scales reach the coefficients only ($m = 1/s_k^X$ on $\beta_k$; the constants and the nest parameters keep $m = 1$), and "bhhh" reaches every parameter, the nest parameters $\lambda$ included, whose lower bound is mapped with the start values. run_nestlogit() keeps "none" as its default: the BHHH scales are fixed at $\theta_0$, and nest parameters that end far from their start can leave them misleading.

1.6 Likelihood Unit: Cross-Section or Panel

By default $i$ indexes a choice situation, and the likelihood integrates a separate coefficient distribution for each situation — the right model when tastes are redrawn every situation, or each decision maker is observed once. Passing person_col to run_mxlogit() / prepare_mxl_data() switches the likelihood unit from the choice situation to the decision maker $n$ (Revelt and Train 1998): the $T_n$ choice situations belonging to $n$ then share one taste draw $\beta_n \sim f(\beta \mid \theta)$, and the likelihood integrates their joint probability, $$L_n(\theta) = \int \prod_{t=1}^{T_n} P_{nt}(j_{nt} \mid \beta)\, f(\beta \mid \theta)\, d\beta,$$ in place of the per-situation $L_i(\theta) = \int P_{ij_i}(\beta)\, f(\beta \mid \theta)\, d\beta$ used elsewhere in this document (§2.2–2.3). id_col must still identify choice situations uniquely across the data set — situation ids do not restart within person. person_col = NULL (the default) makes every situation its own decision maker, $T_n \equiv 1$, and recovers the cross-sectional likelihood exactly, so every formula in §2–§9 written per situation $i$ applies unchanged, read with $i \equiv (n, t{=}1)$.

On panel data, imposing $T_n \equiv 1$ on genuinely repeated tasks (i.e., fitting cross-sectionally) maximizes a composite (pseudo-) likelihood: each situation’s marginal probability is still correctly specified under $f(\beta \mid \theta)$, so the estimator stays consistent for $\theta$, provided the per-situation marginals identify $\theta$, but it is statistically inefficient — it never uses the fact that the same $\beta_n$ generated a person’s repeated choices, which is exactly the variation that identifies $\Sigma$ most sharply (the marginals identify it too, but only through the shape of the mixture). Cluster-robust standard errors (cluster_col=, or vcov(fit, type = "cluster"); §9.5) repair the inference for the resulting score dependence but leave the likelihood unchanged; §2.6, §3.5 and §4.5 give the panel likelihood, score and Hessian that use the repetition directly. run_hmnlogit() is the Bayesian counterpart of the panel likelihood, adding partially-pooled alternative effects and posterior rather than asymptotic inference.

---

2. Log-Likelihood Function

2.1 Choice Probability

For a given draw $s$, the probability that individual $i$ chooses alternative $j$ follows the logit formula:

$$P_{ij}^s = \frac{\exp(V_{ij}^s)}{\sum_{k=1}^{J_i} \exp(V_{ik}^s)}$$

If an outside option is included, it is normalized to $V_{i0}^s = 0$.

2.2 Simulated Probability

Since the random coefficients $\gamma_i$ are unobserved, we integrate over their distribution using simulation:

$$\bar{P}_{ij_i} = \frac{1}{S} \sum_{s=1}^{S} P_{ij_i}^s$$

where $j_i$ is the alternative chosen by individual $i$.

2.3 Log-Likelihood

The simulated log-likelihood is:

$$\ell(\theta) = \sum_{i=1}^{N} w_i \log\left(\bar{P}_{ij_i}\right) = \sum_{i=1}^{N} w_i \log\left(\frac{1}{S} \sum_{s=1}^{S} P_{ij_i}^s\right)$$

where $\theta = (\beta, \mu, \text{vec}(L), \delta)$ is the full parameter vector.

2.4 Log-Sum-Exp Trick for Numerical Stability

To avoid numerical overflow/underflow, the implementation uses:

1. Probability computation: Before computing $P_{ij}^s$, subtract $\max_k V_{ik}^s$: $$P_{ij}^s = \frac{\exp(V_{ij}^s - \max_k V_{ik}^s)}{\sum_k \exp(V_{ik}^s - \max_k V_{ik}^s)}$$

2. Log-probability averaging: Use log-sum-exp to compute $\log(\sum_s P_{ij_i}^s)$: $$\text{logSumExp}(a, b) = \max(a,b) + \log\left(\exp(a - \max(a,b)) + \exp(b - \max(a,b))\right)$$

   Applied iteratively across draws $s$. With person_col (§2.6), each draw’s log-probability is first summed across a decision maker’s situations, $\lambda_{ns} = \sum_{t=1}^{T_n} \log P_{ntj_{nt}}^s$, and only then does log-sum-exp run over draws $s$ to obtain $\ell_n(\theta) = \mathrm{LSE}_s(\lambda_{n\cdot}) - \log S$.

3. Log-probability of the chosen alternative: each term of $\lambda_{ns}$ is formed directly in log space, as the chosen utility minus the log-sum-exp of the choice set, $$\log P_{ntj}^s = V_{ntj}^s - \max_k V_{ntk}^s - \log \sum_k \exp\left(V_{ntk}^s - \max_k V_{ntk}^s\right),$$ never as the logarithm of $P_{ntj}^s$, which turns subnormal beyond a utility gap of about 708 and underflows to zero beyond about 745 (the multinomial logit kernels use the same form). Neither a single choice probability nor a product of many is ever formed outside log space, so $\lambda_{ns}$ and $\ell_n(\theta)$ are finite, and exact, whenever the utilities are finite, however poorly the parameter vector fits the observed choices.

2.5 Statistical Properties of the Finite-Draw Criterion

For independent pseudo-random draws, $\bar P_{ij_i}$ is an unbiased simulator of the exact mixed-logit probability. Its logarithm is not: Jensen’s inequality gives $\mathbb E[\log\bar P]<\log\mathbb E[\bar P]$ when the simulator is nondegenerate, and a second-order expansion gives the familiar $O(S^{-1})$ criterion bias. Under the regularity conditions in Hajivassiliou and Ruud (1994), Lee (1995), Gouriéroux and Monfort (1996), and Train (2009, Ch. 10), the classical pseudo-random MSL benchmarks are:

• Fixed $S$: MSL is inconsistent — the simulation bias does not vanish as $N \to \infty$.
• $S \to \infty$ with $N$ (at any rate): MSL is consistent.
• $S$ rising faster than $\sqrt{N}$: MSL is asymptotically equivalent to exact maximum likelihood — simulation bias is $o(N^{-1/2})$ and simulation variance is negligible under the maintained draw scheme.

choicer uses a Halton sequence rather than independent pseudo-random draws (§5.4), so the first sentence above is not a literal property of its deterministic store-mode simulator, and the pseudo-random rate conditions should not be quoted as an exact theorem for this implementation. Quasi-Monte Carlo often reduces integration error markedly at a given $S$ (Bhat 2003), but its error depends on integrand smoothness, dimension, sequence construction, and randomization. The operational implication is to treat $S$ and, for generated scrambled draws, the seed as part of the numerical specification; increase $S$ until estimates, the Cholesky block, the objective, and standard errors stabilize. mc_asymptotics() supplies Monte Carlo diagnostics rather than a proof that a particular finite-$S$ fit satisfies exact-ML asymptotics.

2.6 Panel Likelihood

With person_col, §1.6 replaces the per-situation likelihood by the per-decision-maker one. Writing $\lambda_{ns} = \sum_{t=1}^{T_n} \log P_{nt j_{nt}}^s$ for the log-probability of decision maker $n$’s choices at draw $s$ — a sum, over the situations sharing that draw, of the situation-level log-probabilities of §2.1 — the simulated log-likelihood is

$$\ell_n(\theta) = \log\left(\frac1S\sum_{s=1}^S \exp(\lambda_{ns})\right), \qquad \ell(\theta) = \sum_{n=1}^U w_n\, \ell_n(\theta),$$

with $w_n$ a decision-maker weight, constant across $n$’s situations — a requirement of this panel objective, not of §9.1–§9.2’s WESML weights, which vary with the alternative a person chose and can therefore differ across their situations; WESML is accordingly not supported with person_col. Define the posterior draw weight

$$\omega_{ns} = \frac{\exp(\lambda_{ns})}{\sum_{s'=1}^S \exp(\lambda_{ns'})} = \exp\bigl(\lambda_{ns} - \mathrm{LSE}_s(\lambda_{n\cdot})\bigr), \qquad \sum_{s=1}^S \omega_{ns} = 1,$$

the simulated posterior probability that draw $s$ is decision maker $n$’s taste, given the choices they made — computed via the max-shifted log-sum-exp of §2.4, so it never forms the product of many situation probabilities outside log space. §3.5 and §4.5 write the panel score and Hessian in terms of $\omega_{ns}$. Setting $T_n \equiv 1$ collapses $\lambda_{ns}$ to $\log P_{ij_i}^s$ and $\omega_{ns}$ to $P_{ij_i}^s / \sum_{s'} P_{ij_i}^{s'}$, the per-draw share already implicit in §2.3’s $\bar P_{ij_i} = \frac1S\sum_s P_{ij_i}^s$, so the panel and cross-sectional log-likelihoods agree termwise.

---

3. Gradient Computation

3.1 General Formula

Define for draw $s$ the log-probability of the chosen alternative:

$$\log P_{ij_i}^s = V_{ij_i}^s - \log\left(\sum_k \exp(V_{ik}^s)\right)$$

The gradient of this with respect to utility is:

$$\frac{\partial \log P_{ij_i}^s}{\partial V_{ia}^s} = \mathbf{1}_{a = j_i} - P_{ia}^s$$

3.2 Per-Draw Gradient

Define $z_{ij}^s = \frac{\partial V_{ij}^s}{\partial \theta}$ as the gradient of utility with respect to parameters. The per-draw gradient is:

$$g_i^s = \frac{\partial \log P_{ij_i}^s}{\partial \theta} = \sum_{j=1}^{J_i} \left(\mathbf{1}_{j = j_i} - P_{ij}^s\right) z_{ij}^s$$

Vectorized implementation. Define the logit residual vector $\mathbf{d}_i^s \in \mathbb{R}^{J_i}$ with elements $d_{ij}^s = \mathbf{1}_{j = j_i} - P_{ij}^s$. The per-block gradient components collapse to matrix-vector products — computed once per draw rather than looped over alternatives:

Block	Scalar sum	Vectorized form
$\beta$	$\sum_j d_{ij}^s X_{ij}^T$	$X_i^T \mathbf{d}_i^s$
$\mu_p$	$\sum_j d_{ij}^s W_{ijp} \cdot \partial\mu_p^*/\partial\mu_p$	$(\tilde{w}_i^s)_p \cdot \partial\mu_p^*/\partial\mu_p$
$\ell_{pq}$	$\sum_j d_{ij}^s W_{ijp} \cdot \partial\gamma_p^{s*}/\partial\gamma_p^s \cdot \partial L_{pq}/\partial\ell_{pq} \cdot \eta_q^s$	$(\tilde{w}_i^s)_p \cdot \partial\gamma_p^{s*}/\partial\gamma_p^s \cdot \partial L_{pq}/\partial\ell_{pq} \cdot \eta_q^s$

where $\tilde{w}_i^s = W_i^T \mathbf{d}_i^s \in \mathbb{R}^{K_w}$ (Wt_diff in the code) is computed once and shared across the $\mu$ and $L$ blocks. The $\delta$ block uses an irregular alt-index scatter and is kept as a loop.

Code reference: src/mxlogit.cpp

3.3 Gradient of Simulated Log-Likelihood

Using the identity for the derivative of a logarithm of a sum:

$$\frac{\partial}{\partial \theta} \log\left(\frac{1}{S}\sum_s P_{ij_i}^s\right) = \frac{\sum_s P_{ij_i}^s \cdot g_i^s}{\sum_s P_{ij_i}^s}$$

The full gradient is:

$$\nabla_\theta \ell = \sum_{i=1}^{N} w_i \cdot \frac{\sum_s P_{ij_i}^s \cdot g_i^s}{\sum_s P_{ij_i}^s}$$

Code reference: src/mxlogit.cpp

3.4 Utility Derivatives by Parameter Block

3.4.1 Fixed Coefficients ($\beta$)

$$\frac{\partial V_{ij}^s}{\partial \beta} = X_{ij}^T$$

Summing over alternatives: $\sum_j d_{ij}^s X_{ij}^T = X_i^T \mathbf{d}_i^s$.

Code reference: src/mxlogit.cpp

3.4.2 Random Coefficient Means ($\mu$)

For Normal distribution: $$\frac{\partial V_{ij}^s}{\partial \mu_p} = W_{ijp}$$

For Log-Normal distribution: $$\frac{\partial V_{ij}^s}{\partial \mu_p} = W_{ijp} \cdot \exp(\mu_p)$$

Both cases are captured by dmu_final_dmu(p) ($= 1$ for normal, $= \exp(\mu_p)$ for log-normal). Summing over alternatives: $g_\mu = \tilde{w}_i^s \odot \text{dmu\_final\_dmu}$.

Code reference: src/mxlogit.cpp

3.4.3 Cholesky Parameters ($L$)

Let $\ell_{pq}$ denote the parameter for $L_{pq}$ (where $p \geq q$).

Step 1: Derivative of $L_{pq}$ with respect to $\ell_{pq}$: $$\frac{\partial L_{pq}}{\partial \ell_{pq}} = \begin{cases} L_{pp} = \exp(\ell_{pp}) & \text{if } p = q \text{ (diagonal)} \\ 1 & \text{if } p > q \text{ (off-diagonal)} \end{cases}$$

Step 2: Derivative of $\gamma_p^s$ with respect to $\ell_{pq}$: $$\frac{\partial \gamma_p^s}{\partial \ell_{pq}} = \frac{\partial L_{pq}}{\partial \ell_{pq}} \cdot \eta_q^s$$

Note: $\gamma_p^s = \sum_{r \leq p} L_{pr} \eta_r^s$, so only the term with $r = q$ contributes.

Step 3: Derivative of transformed $\gamma_p^{s*}$:

For Normal distribution: $$\frac{\partial \gamma_p^{s*}}{\partial \gamma_p^s} = 1$$

For Log-Normal distribution: $$\frac{\partial \gamma_p^{s*}}{\partial \gamma_p^s} = \exp(\gamma_p^s) = \gamma_p^{s*}$$

Step 4: Chain rule for utility derivative: $$\frac{\partial V_{ij}^s}{\partial \ell_{pq}} = W_{ijp} \cdot \frac{\partial \gamma_p^{s*}}{\partial \gamma_p^s} \cdot \frac{\partial \gamma_p^s}{\partial \ell_{pq}}$$

Step 5: Summing over alternatives using $\tilde{w}_i^s = W_i^T \mathbf{d}_i^s$: $$g_{\ell_{pq}}^s = (\tilde{w}_i^s)_p \cdot \frac{\partial \gamma_p^{s*}}{\partial \gamma_p^s} \cdot \frac{\partial L_{pq}}{\partial \ell_{pq}} \cdot \eta_q^s$$

For the diagonal case, all $K_w$ parameters are computed simultaneously as Wt_diff % dgamma_final_dgamma % L.diag() % eta_i_s. For the full correlated case, Wt_p = Wt_diff(p) * dgamma_final_dgamma(p) is factored out per row $p$.

Code reference: src/mxlogit.cpp

3.4.4 Alternative-Specific Constants ($\delta$)

$$\frac{\partial V_{ij}^s}{\partial \delta_a} = \mathbf{1}_{j = a}$$

where $\delta_1 = 0$ (normalized) if no outside option, or all $\delta_j$ for $j \geq 1$ are free if outside option is included. The sum over alternatives is a simple scatter (one non-zero contribution per alternative), implemented as a loop.

Code reference: src/mxlogit.cpp

3.5 Panel Score and Posterior Draw Weights

With person_col, the per-draw score of §3.2 is first summed across a decision maker’s situations, $g_{ns} = \sum_{t=1}^{T_n} g_{nts}$ (each $g_{nts}$ the situation-level score of §3.2–§3.4), and then averaged over draws with the posterior weight of §2.6 in place of the per-draw probability share of §3.3:

$$\nabla_\theta \ell_n = \sum_{s=1}^S \omega_{ns}\, g_{ns}, \qquad \nabla_\theta \ell = \sum_{n=1}^U w_n \nabla_\theta \ell_n.$$

This is §3.3’s $\sum_s P_{ij_i}^s g_i^s / \sum_s P_{ij_i}^s$ with $P_{ij_i}^s / \sum_{s'} P_{ij_i}^{s'}$ replaced by $\omega_{ns}$ and $g_i^s$ replaced by the situation-summed $g_{ns}$; $T_n \equiv 1$ recovers §3.3 termwise. The block decomposition of §3.4 (the $\beta$/$\mu$/$L$ BLAS collapse via $\tilde w^s$, and the $\delta$ scatter) is unchanged within a situation.

Implementation. The double sum over situations and draws needs no explicit loop over $t$: mxl_unit_simulate() already stacks the logit residual $d_{nt}^s$ of every row of every situation belonging to unit $n$ into one $R_n \times S$ matrix DiffW ($R_n$ = stacked alternative rows across all of $n$’s $T_n$ situations), because $\omega_{ns}$ is a single per-draw weight shared by every situation in the unit. $\bar d_n = \texttt{DiffW}\,\omega_n$ then gives the $\beta$, $\mu$ and $\delta$ blocks directly, exactly as §3.2’s vectorized table builds them from $\mathbf d_i^s$ — the situation sum is absorbed into the row-stacking, and the draw sum into the matrix-vector product, in a single pass.

DiffW is formed $B$ draws at a time, $R_n \times B$, with $B = S$ unless the two $R_n \times S$ matrices (DiffW and $W_n \Gamma_n$) would exceed a small per-thread budget (4 MiB, sized to stay in cache): a unit is split into batches (of sizes differing by at most one draw) only when $R_n S > 2^{18}$, i.e. for a long panel, a choice set of hundreds of alternatives or a very large $S$, and a thread’s memory then grows with $R_n$ rather than $R_n S$. Each batch is folded in before the next one overwrites it. $\bar d_n$ accumulates by a streaming max shift: the batch’s weights $e^{\lambda_{ns} - r}$ are taken relative to the largest $\lambda_{ns}$ seen so far. The running sum is rescaled by $e^{r_\text{old} - r}$ whenever $r$ rises, then divided by $\sum_s e^{\lambda_{ns} - r}$ at the end, the same denominator used to normalize $\omega_n$. Keeping the maximum separate from this denominator avoids losing normalization accuracy through subtraction of a large-magnitude log-sum-exp. Batching changes $\bar d_n$ only by rounding (about one rounding per batch), and a unit in one batch computes exactly the unbatched values. The $L$ block cannot collapse through $\bar d_n$ alone: its per-draw derivative couples the residuals with that draw’s $\eta_{ns}$ (§3.4.3), so the draw sum has to happen after that coupling; each batch writes its columns of $W_n'\,\texttt{DiffW}$ ($K_w \times S$, small whatever $R_n$ is), and mxl_unit_score_cont() then forms $A = \bigl((W_n' \,\texttt{DiffW} \odot D\gamma_1)\,\mathrm{diag}(\omega_n)\bigr)\,\eta_n'$ ($K_w \times K_w$; a situation-stacked dgemm for the inner product, then a draw-block dgemm against $\eta_n'$), and reads $\ell_{pq}$’s score off $\partial L_{pq}/\partial \ell_{pq} \cdot A(p, q)$. (The Hessian, §4.5, needs an explicit per-situation loop because its outer-product terms cannot be pre-summed across situations before being squared.)

Code reference: mxl_unit_simulate(), mxl_unit_score_cont() (the continuous block, with the $A$ collapse), mxl_unit_score() and mxl_unit_score_local() (the delta block, over every parameter or the unit’s own, §4.4.6) in src/mxlogit.cpp, shared by mxl_loglik_gradient_parallel, mxl_scores_parallel, mxl_bhhh_parallel and mxl_cluster_meat_parallel.

---

4. Hessian Computation

4.1 General Formula

Define: - $g_i = \frac{\partial \log \bar{P}_i}{\partial \theta}$ — the gradient for individual $i$ - $g_i^s = \frac{\partial \log P_{ij_i}^s}{\partial \theta}$ — the per-draw gradient - $H_i^s = \frac{\partial^2 \log P_{ij_i}^s}{\partial \theta \partial \theta'}$ — the per-draw Hessian

The Hessian of the simulated log-probability for individual $i$ is:

$$H_i = \frac{\partial^2 \log \bar{P}_i}{\partial \theta \partial \theta'} = \frac{\sum_s P_{ij_i}^s \left(g_i^s (g_i^s)' + H_i^s\right)}{\sum_s P_{ij_i}^s} - g_i g_i'$$

The full Hessian is:

$$\nabla^2_\theta \ell = \sum_{i=1}^{N} w_i \cdot H_i$$

Code reference: src/mxlogit.cpp

4.2 Per-Draw Hessian

The per-draw Hessian $H_i^s$ has the form:

$$H_i^s = -\sum_j P_{ij}^s z_{ij}^s (z_{ij}^s)' + \left(\sum_j P_{ij}^s z_{ij}^s\right)\left(\sum_j P_{ij}^s z_{ij}^s\right)' + \sum_j \left(\mathbf{1}_{j=j_i} - P_{ij}^s\right) H_{V,ij}^s$$

where $H_{V,ij}^s = \frac{\partial^2 V_{ij}^s}{\partial \theta \partial \theta'}$ is the Hessian of utility.

In the code, this is computed as: - sum_Pzz $= \sum_j P_{ij}^s z_{ij}^s (z_{ij}^s)'$ - sum_Pz $= \sum_j P_{ij}^s z_{ij}^s$ - sum_diff_H_V $= \sum_j (\mathbf{1}_{j=j_i} - P_{ij}^s) H_{V,ij}^s$

Then: $H_i^s = -\text{sum\_Pzz} + \text{sum\_Pz} \cdot \text{sum\_Pz}' + \text{sum\_diff\_H\_V}$

Code reference: src/mxlogit.cpp

4.3 Second Derivatives of Utility

Most second derivatives are zero since utility is linear in most parameters. Non-zero second derivatives arise from:

4.3.1 Mean Parameters ($\mu$) — Log-Normal Only

For log-normal distribution: $$\frac{\partial^2 V_{ij}^s}{\partial \mu_p^2} = W_{ijp} \cdot \exp(\mu_p)$$

Code reference: src/mxlogit.cpp

4.3.2 Cholesky Diagonal Parameters ($\ell_{pp}$)

Since $L_{pp} = \exp(\ell_{pp})$: $$\frac{\partial^2 L_{pp}}{\partial \ell_{pp}^2} = \exp(\ell_{pp}) = L_{pp}$$

For normal distribution: $$\frac{\partial^2 V_{ij}^s}{\partial \ell_{pp}^2} = W_{ijp} \cdot L_{pp} \cdot \eta_p^s$$

For log-normal distribution (combining chain rule): $$\frac{\partial^2 V_{ij}^s}{\partial \ell_{pp}^2} = W_{ijp} \left[ \gamma_p^{s*} \left(\frac{\partial \gamma_p^s}{\partial \ell_{pp}}\right)^2 + \frac{\partial \gamma_p^{s*}}{\partial \gamma_p^s} \cdot \frac{\partial^2 \gamma_p^s}{\partial \ell_{pp}^2} \right]$$

where $\frac{\partial^2 \gamma_p^s}{\partial \ell_{pp}^2} = L_{pp} \cdot \eta_p^s$.

Code reference: src/mxlogit.cpp

4.3.3 Cholesky Cross-Derivatives Within a Row

Let $q\neq r$ with $q,r\leq p$ and define

$$a_{pq}^s=\frac{\partial\gamma_p^s}{\partial\ell_{pq}} =d_{pq}\eta_q^s, \qquad d_{pq}=\begin{cases}L_{pp}&q=p,\\1&q<p.\end{cases}$$

Distinct Cholesky parameters affect the same transformed coefficient only when they belong to the same row $p$. Their cross derivative is

$$\frac{\partial^2 V_{ij}^s}{\partial \ell_{pq} \partial \ell_{pr}} =W_{ijp}\,f_p''(\gamma_p^s)\,a_{pq}^s a_{pr}^s,$$

where $f_p(\gamma)=\gamma$ for a normal coefficient and $f_p(\gamma)=\exp(\gamma)$ for a log-normal coefficient. Hence the cross derivative is zero for normal coefficients and equals $W_{ijp}\exp(\gamma_p^s)a_{pq}^sa_{pr}^s$ for log-normal coefficients. This expression also covers a diagonal/off-diagonal pair (say $q=p,r<p$); restricting both indices to be off-diagonal would omit terms that the C++ Hessian includes. Cross derivatives across different rows of $L$ are zero. The separate shifted term $\exp(\mu_p)$ implies no $\mu$–$L$ cross derivative.

Code reference: src/mxlogit.cpp

---

4.4 Implementation Tricks in mxl_hessian_parallel

The subsections below describe the computational restructuring used by mxl_hessian_parallel. They cover (i) a block decomposition of the score vector and Hessian terms, (ii) buffer accumulation that eliminates the per-draw dense matrix allocation, (iii) BLAS-3 batching of two outer-product sums, (iv) the symmetry strategy used for final assembly, and (v) the compact block of the parameters a unit’s likelihood involves, in which (i)–(iv) run, and how the blocks are accumulated into the result.

A reader who has never seen the C++ can understand what the function does and why the result is numerically equivalent to the per-draw naive form (§4.1–4.2) from these sections alone.

4.4.1 Parameter Block Layout and Score Vector Decomposition

The parameter vector $\theta$ is partitioned into a continuous block of length $K_c$ and a delta block (ASCs) of length $J_d$:

$$\theta = \bigl(\underbrace{\beta,\; \mu,\; \mathrm{vec}(L)}_{\text{continuous, length } K_c},\; \underbrace{\delta}_{\text{delta, length } J_d}\bigr), \qquad n_\text{params} = K_c + J_d.$$

The boundary index $K_c = K_x + (\mathbf{1}_{\mathrm{rc\_mean}} \cdot K_w) + L\_\text{size}$ is stored as idx_delta_start in the code.

The per-draw score vector $g_{is} \in \mathbb{R}^{n_\text{params}}$ inherits this partition: $g_{is} = [g_{is,c};\; g_{is,d}]$. The continuous sub-vector is

$$g_{is,c} = \sum_{a} d_{ia}^s\, z_{c,ia}^s, \qquad d_{ia}^s = \mathbf{1}_{a = j_i} - P_{ia}^s,$$

where $z_{c,ia}^s \in \mathbb{R}^{K_c}$ is the continuous-block utility gradient for alternative $a$ at draw $s$. The delta sub-vector scatters logit residuals $d_{ia}^s$ into the ASC index for each alternative.

Define also: - $\mathrm{sum\_Pz}_{is} = \sum_a P_{ia}^s z_{c,ia}^s \in \mathbb{R}^{K_c}$ (the cc weighted-score sum, plus a $J_d$-length delta counterpart $\mathrm{sum\_Pz\_d}_{is}$), - $\mathrm{sum\_Pzz}_{cc,is} = \sum_a P_{ia}^s z_{c,ia}^s (z_{c,ia}^s)^\top$ ($K_c \times K_c$, symmetric), and its cd ($K_c \times J_d$) and dd (diagonal $J_d$-vector) sub-blocks.

4.4.2 Block Buffer Accumulation — Identity C

Two of the three terms in $H_{is}$ (§4.2) are linear in the draw weight $P_{is}$:

$$-\mathrm{sum\_Pzz}_{is} \quad \text{and} \quad +\mathrm{sum\_diff\_H\_V}_{is}.$$

Their draw-weighted sum can therefore be accumulated directly into per-individual buffer matrices without constructing the full $n_\text{params} \times n_\text{params}$ dense matrix $H_{is}$ at each draw:

$$\mathrm{buf\_Pzz} = \sum_{s=1}^{S} P_{is} \cdot \mathrm{sum\_Pzz}_{is}, \qquad \mathrm{buf\_diff\_HV} = \sum_{s=1}^{S} P_{is} \cdot \mathrm{sum\_diff\_HV}_{is}.$$

Each buffer is updated at the end of draw $s$ with a single scalar-times-matrix addition (buf += P_s * block_s), replacing the original path that assembled the full $H_{is}$ and then computed hess_term1_numerator += P_s * (g_is * g_isᵀ + H_is).

The block structure of buf_Pzz minimizes memory by storing only the three structurally distinct sub-blocks:

Buffer	Shape	Content	Notes
buf_Pzz_cc	$K_c \times K_c$	$\sum_s P_{is}\,\mathrm{sum\_Pzz}_{cc,s}$	Symmetric
buf_Pzz_cd	$K_c \times J_d$	$\sum_s P_{is}\,\mathrm{sum\_Pzz}_{cd,s}$	Rectangular; not symmetric
buf_Pzz_dd	$J_d$-vector	$\sum_s P_{is}\,\mathrm{sum\_Pzz}_{dd,s}$	Diagonal only
buf_diff_HV_cc	$K_c \times K_c$	$\sum_s P_{is}\,\mathrm{sum\_diff\_HV}_{s}$	Symmetric; cc block only ($H_V$ has no delta rows)

All four buffers are thread-private scratch, zeroed at the start of each individual $i$. In the code the cd and dd buffers hold only the constants that $i$’s likelihood involves (§4.4.6), so they grow only to the largest such block their thread meets.

4.4.3 BLAS-3 Batching of Outer-Product Sums — Identities A and B

The third term in $H_{is}$ — $+\mathrm{sum\_Pz}_{is}\,\mathrm{sum\_Pz}_{is}^\top$ — is not linear in $P_{is}$, so it cannot go into a linear buffer. Instead, together with the OPG term $\sum_s P_{is} g_{is} g_{is}^\top$, it is handled by reformulating both weighted outer-product sums as a single BLAS level-3 matrix multiply.

Identity A (OPG batching). With the normalized draw weights $\omega_{is} = P_{is} / \hat P_i$, $\hat P_i = \sum_s P_{is}$, and the weighted mean score $g_i = \sum_s \omega_{is} g_{is}$, define the $n_\text{params} \times S$ matrix $G$ with columns

$$G_{[:,s]} = \sqrt{\omega_{is}}\; (g_{is} - g_i).$$

Then

$$\sum_{s=1}^{S} \omega_{is}\, g_{is} g_{is}^\top - g_i g_i^\top = G G^\top.$$

Proof. $G G^\top = \sum_s \omega_{is} (g_{is} - g_i)(g_{is} - g_i)^\top = \sum_s \omega_{is} g_{is} g_{is}^\top - g_i g_i^\top$, because $\sum_s \omega_{is} = 1$. The identity requires only $\omega_{is} \geq 0$, which holds because $P_{is}$ is a softmax probability (guaranteed by stable_softmax). It folds the OPG term and the final $- g_i g_i^\top$ of the Hessian (§4.4.5) into one product, and the centered form is also the numerically stable one. The implementation normalizes max-shifted exponentials directly, then subtracts the highest-weight draw’s score before computing the weighted mean of the deviations. This prevents a large common score from amplifying rounding in the weights into spurious covariance.

Identity B (sum-Pz outer-product batching). Define the $n_\text{params} \times S$ matrix $F$ with columns

$$F_{[:,s]} = \sqrt{\omega_{is}}\;\bigl[\mathrm{sum\_Pz\_c}_{is};\;\mathrm{sum\_Pz\_d}_{is}\bigr].$$

Then

$$\sum_{s=1}^{S} \omega_{is}\;\mathrm{sum\_Pz}_{is}\,\mathrm{sum\_Pz}_{is}^\top = F F^\top,$$

by the same expansion without the centering. The linear buffers of §4.4.2 are likewise accumulated with $\omega_{is}$ in place of $P_{is}$, which divides them by $\hat P_i$ once and for all.

Combined BLAS-3 form. Concatenating $G$ and $F$ into $[G \mid F] \in \mathbb{R}^{n_\text{params} \times 2S}$ and forming one matrix product gives

$$[G \mid F]\;[G \mid F]^\top = G G^\top + F F^\top,$$

because $[G \mid F]^\top = [G^\top;\; F^\top]$ and the block multiply sums over $2S$ columns with no $G$-$F$ cross term. This is the form implemented in mxl_hessian_parallel: G_stash and F_stash are filled column-by-column during the $S$-loop, and opg_pz accumulates $F F^\top$ and $G G^\top$ as two separate rank-$S$ BLAS level-3 updates rather than one join_rows product (§4.5 explains why: with person_col, $F$ is squared and added once per situation while $G$ accumulates across situations and is squared once per unit; the two coincide here because §4.4 is the single-situation, $T_n \equiv 1$ case). The result is $\sum_s P_{is}(g_{is} g_{is}^\top + \mathrm{sum\_Pz}_{is}\,\mathrm{sum\_Pz}_{is}^\top)$ in cache-friendly level-3 kernels.

Floating-point precision. Changing the summation order (per-draw increments vs. one level-3 multiply) is a reassociation of additions that is exact in real arithmetic. In floating-point the results differ only by rounding. Across 10 equivalence scenarios sweeping all combinations of rc_dist, rc_correlation, and rc_mean, the worst-case element-wise deviation from the per-draw form was $4.588 \times 10^{-13}$, several orders of magnitude below the $10^{-6}$ tolerance enforced by the package’s equivalence tests.

4.4.4 Optional Per-Draw Batching of sum_Pzz_cc — Identity D

Within each draw $s$, the cc sub-block of sum_Pzz accumulates per-alternative rank-1 updates:

$$\mathrm{sum\_Pzz\_cc}_{s} = \sum_{a} P_{ia}^s\, z_{c,ia}^s (z_{c,ia}^s)^\top.$$

This can equivalently be expressed as a single matrix product:

$$\mathrm{sum\_Pzz\_cc}_{s} = \tilde{Z}_{c}^s\, (\tilde{Z}_{c}^s)^\top, \qquad \tilde{Z}_{c,[:,a]}^s = \sqrt{P_{ia}^s}\; z_{c,ia}^s,$$

where $\tilde{Z}_c^s \in \mathbb{R}^{K_c \times J_{\text{eff}}^*}$ and $J_{\text{eff}}^*$ counts the active (non-outside-option) alternatives. One dsyrk call per draw over $J_{\text{eff}}^*$ columns would replace $J_{\text{eff}}^*$ individual dsyr rank-1 updates. The identity holds by the same $P_{ia}^s \geq 0$ argument as Identities A and B.

Implementation choice. The implementation retains the per-alternative rank-1 accumulation (sum_Pzz_cc += P_a * (zc_a * zc_a.t())). Both forms are correct; the result flows into buf_Pzz_cc via the Identity C scalar-times-matrix step at the end of each draw. Identity D is recorded here to establish that the batched form is available without changing any result; the per-alternative form is kept for auditability.

4.4.5 Final Assembly and Symmetry

After the $S$-loop for individual $i$, let $\mathrm{opg\_pz} = [G \mid F][G \mid F]^\top$, with the weights of §4.4.3. The Hessian $\sum_s \omega_{is}(H_{is} + (g_{is} - g_i)(g_{is} - g_i)^\top)$ is assembled block-by-block into hess_t1 (the code builds it in place in opg_pz, reading each block before overwriting it):

$$\mathrm{hess\_t1}\big|_{\mathrm{cc}} = \mathrm{opg\_pz}_{\mathrm{cc}} - \mathrm{buf\_Pzz\_cc} + \mathrm{buf\_diff\_HV\_cc},$$

$$\mathrm{hess\_t1}\big|_{\mathrm{cd}} = \mathrm{opg\_pz}_{\mathrm{cd}} - \mathrm{buf\_Pzz\_cd}, \qquad \mathrm{hess\_t1}\big|_{\mathrm{dc}} = \bigl(\mathrm{hess\_t1}|_{\mathrm{cd}}\bigr)^\top,$$

$$\mathrm{hess\_t1}\big|_{\mathrm{dd}} = \mathrm{opg\_pz}_{\mathrm{dd}} - \mathrm{diag}(\mathrm{buf\_Pzz\_dd}).$$

The cd block is rectangular ($K_c \times J_d$) and not symmetric, so both the cd and dc sub-matrices are written explicitly into hess_t1. The normalized draw weight $\omega_{is} \equiv P_{is}/\hat P_i$ replaces a separate division by $\hat P_i$, and the centered $G$ of Identity A already carries $- g_i g_i^\top$, matching §4.5’s Louis-identity assembly:

$$H_i = \mathrm{hess\_t1}, \qquad \text{result} \mathrel{+}= w_i H_i \quad \text{(upper triangle, §4.4.6)}.$$

Symmetry strategy. $H_i$ is exactly symmetric in floating point, not only in exact arithmetic: $G G^\top$ and $F_t F_t^\top$ are symmetric products, buf_Pzz_cc and buf_diff_HV_cc are accumulated as full symmetric matrices (the per-alternative loop writes both $(r, p)$ and $(p, r)$ for off-diagonal $L$-$L$ entries in sum_diff_H_V_cc), and the dc block is written as the transpose of the cd block. So only the upper triangle of $w_i H_i$ is added into the result, which is mirrored once after the parallel region (§4.4.6); the mirror is the sum the lower triangle would have accumulated.

4.4.6 Compact Blocks and Accumulation

A likelihood unit’s score and Hessian are zero outside the continuous block and the constants of the alternatives the unit faced (unit $u$: a decision maker, or a choice situation in the cross-section, §1.6). Let $J_{d,u}$ be the number of free ASCs among the alternatives in $u$’s choice situations (the normalized reference has none without an outside option, and the outside option has none), and $J_{d,t} \le J_{d,u}$ the number in situation $t$. Everything in §4.4.1–§4.4.5 runs in $u$’s compact block of $K_c + J_{d,u}$ local parameters, the continuous parameters first and then $u$’s free ASCs in increasing global order, so that each $J_d$-sized piece there is restricted to the constants $u$’s likelihood involves. Each situation has its own slots as well: the per-situation, per-draw pieces (sum_Pzz_cd, sum_Pzz_dd, $F_t$) are formed on its $K_c + J_{d,t}$ parameters and added into $u$’s slots, and an alternative listed twice in a situation adds into one slot, as the dense scatter did. The product $G G^\top$ then costs $O(S\,(K_c+J_{d,u})^2)$ per unit and each $F_t F_t^\top$ costs $O(S\,(K_c+J_{d,t})^2)$ per situation, instead of $O(S\,n_\text{params}^2)$: on a school-census design with 2,620 parameters, a student’s likelihood involves 59 of them on average.

$w_u H_u$ is added into the $n_\text{params} \times n_\text{params}$ result at its global indices, upper triangle only (§4.4.5). With one thread the blocks are added into the result itself. With several, each thread adds into its own packed upper triangle while the threads’ triangles, $n_\text{params}(n_\text{params}+1)/2$ doubles each, fit in 2 GiB together (up to about 7,000 parameters at eleven threads), and the triangles are added up in thread order. Above that the threads share the result: each keeps private only the rows of the $K_c$ continuous parameters, which every unit touches, and adds the ASC block atomically. The two differ only in the order of the additions. Per-thread triangles are kept where they fit because they need no atomic additions: at census scale with eleven threads, the clustered meat of 200,000 students in 2,000 clusters (§9.5) took 1.37 s with them and 1.82 s with a shared result, and the BHHH matrix of 200,000 students 1.23 s and 1.32 s; the Hessian took the same time either way, and only with two draws, where accumulation dominates, was the shared result faster (the BHHH matrix in 0.47 s against 0.52 s). The BHHH kernel’s rank-one updates $w_u g_u g_u^\top$ (§4.1’s score) use the same blocks and accumulators, and the cluster meat’s update with each cluster’s summed weighted scores (§9.5) the union of its units’ blocks.

Non-finite terms spread as the dense products spread them: a non-finite row of $F_t$ or of the centered $G$ makes the rows and columns of the parameters concerned NaN in the result ($0 \cdot \mathrm{NaN}$), and a non-finite weight makes every entry outside the unit’s block NaN ($w \cdot 0$; NA for a missing weight), or, with several such units, every entry outside the intersection of their blocks.

4.5 Panel Hessian: the Louis (1982) Identity and Its Two-Pass Assembly

With person_col, §4.1’s per-situation Hessian generalizes to decision maker $n$ by the Louis (1982) missing-information identity, treating the draw $s$ as missing data and $n$’s choices as the observed data:

$$H_n = \nabla^2_\theta \ell_n = \sum_{s=1}^S \omega_{ns}\bigl(H_{ns} + g_{ns} g_{ns}^\top\bigr) - \bar g_n \bar g_n^\top, \qquad \bar g_n = \sum_{s=1}^S \omega_{ns}\, g_{ns}, \qquad \nabla^2_\theta \ell = \sum_{n=1}^U w_n H_n,$$

with $g_{ns} = \sum_t g_{nts}$ as in §3.5 and $H_{ns} = \sum_t H_{nts}$ the situation-summed per-draw Hessian of §4.2. Because $\omega_{ns}$ is already normalized ($\sum_s \omega_{ns} = 1$, §2.6), this is $H_n$ directly, with no separate division step — the same normalized form as §4.4.5’s $\omega_{is} \equiv P_{is}/\hat P_i$ (defined there). Setting $T_n \equiv 1$ makes §4.4.5’s formula the special case of the one above, term by term.

Two passes per decision maker. Unlike the score (§3.5), the Hessian cannot be read off the situation-stacked residual matrix alone: its outer-product terms must be squared within a situation before being summed across situations, since summing residuals across situations first and squaring once would cross-multiply different situations’ $\mathrm{sum\_Pz}$ vectors — an identity the algebra does not license. mxl_hessian_parallel therefore runs two passes per unit, sharing the same draw loop that computes $\omega_{ns}$ (§2.6):

• Pass 1 is mxl_unit_simulate() without residuals: it yields $\lambda_{ns}$ (the per-draw log-likelihoods of the unit’s choices), their log-sum-exp over draws, and $\omega_{ns}$ — exactly as in §3.5.
• Pass 2 loops situations $t$ outer, draws $s$ inner, repeating §4.4’s per-alternative block accumulation with $\omega_{ns}$ everywhere §4.4 used $P_{is}$. A unit whose draws pass 1 had to split into batches ($R_n S > 2^{19}$) does not keep the $R_n \times S$ matrix $W_n \Gamma_n$: each situation’s $m_t \times S$ rows of it are formed in turn, so a thread’s memory grows with the largest choice set rather than with the unit:
  • Identity C (linear pieces): buf_Pzz_cc/cd/dd and buf_diff_HV_cc accumulate $\omega_{ns}\,\mathrm{sum\_Pzz}_{nts}$ and $\omega_{ns}\,\mathrm{sum\_diff\_HV}_{nts}$ at every $(t,s)$ pair in the unit — linear accumulation does not care whether it is summed within or across situations.
  • Identity B (F_stash, the $\mathrm{sum\_Pz}\,\mathrm{sum\_Pz}^\top$ term): column $s$ is overwritten each situation with $\sqrt{\omega_{ns}}\,\mathrm{sum\_Pz}_{nts}$, and $F_t F_t^\top$, formed on the situation’s slots (§4.4.6), is added into opg_pz once per situation, immediately after that situation’s $S$-loop — situations’ outer products are added, never combined before squaring: $\sum_t F_t F_t^\top = \sum_t \sum_s \omega_{ns}\,\mathrm{sum\_Pz}_{nts}\, \mathrm{sum\_Pz}_{nts}^\top$.
  • Identity A (G_stash, the score outer product): column $s$ instead accumulates $g_{nts}$ across every situation in the unit, and only after the situation loop ends is it centered at $\bar g_n$, scaled by $\sqrt{\omega_{ns}}$ and squared once: $G G^\top = \sum_s \omega_{ns} (g_{ns} - \bar g_n)(g_{ns} - \bar g_n)^\top$ with $g_{ns} = \sum_t g_{nts}$ — the situation sum must happen before squaring here, the opposite order from Identity B, because it is the person’s total score that enters the outer product, not each situation’s on its own. Since the $\omega_{ns}$ sum to one, this is the Louis identity’s $\sum_s \omega_{ns} g_{ns} g_{ns}^\top - \bar g_n \bar g_n^\top$, computed without its cancellation. In code, g_bar first holds a highest-weight draw’s score, which is subtracted from every column; it then holds the weighted mean of those deviations, subtracted in turn. This is equivalent to centering at $\bar g_n$ while avoiding amplification of weight-rounding error by a large common score.

Final assembly mirrors §4.4.5’s block-by-block combination (hess_t1 = opg_pz - buf_Pzz + buf_diff_HV, cc/cd/dc/dd), whose centered Identity A already carries the Louis identity’s $-\bar g_n \bar g_n^\top$, and closes by adding the upper triangle of $w_n\,$hess_t1 into the result (§4.4.6) — no $/\hat P_n$ step, because $\omega_{ns}$ is already a normalized weight (§2.6), as in §4.4.5. Setting $T_n \equiv 1$ collapses pass 2 to a single situation per unit: F_stash (holding $\sqrt{\omega_{ns}}\,\mathrm{sum\_Pz}_{nts}$) and G_stash (holding $\sqrt{\omega_{ns}}\, g_{nts}$) remain different vectors, but with only one situation to loop over, F’s per-situation overwrite and G’s cross-situation accumulate produce the same result, and the two-pass structure reduces exactly to §4.4’s single-situation accumulation.

The non-finite-likelihood guard. Pass 1 skips a unit only when its log-sum-exp is non-finite, !std::isfinite(lse). Because every $\log P_{ntj}^s$ is formed in log space (§2.4), lse is finite whenever the unit’s utilities are, however poorly the unit is fit; the guard fires only for utilities that overflow, where the unit’s log-likelihood is undefined. It replaces the former $\hat P_i^{\mathrm{sum}} \le 10^{-12}$ threshold of §5.4, which could also discard a unit that was merely very poorly fit under the current parameter vector. See §5.4 and NEWS.md.

Code reference: mxl_hessian_parallel() in src/mxlogit.cpp.

---

5. Implementation Details

5.1 Parameter Vector Structure

The full parameter vector $\theta$ is organized as:

Block	Indices	Length	Description
$\beta$	$[0, K_x)$	$K_x$	Fixed coefficients
$\mu$	$[K_x, K_x + K_w)$	$K_w$ (if rc_mean=TRUE)	Random coefficient means
$L$	next block	$K_w(K_w+1)/2$ or $K_w$	Cholesky parameters
$\delta$	remaining	$0$, $J-1$, or $J$	Free ASCs

The ASC length is zero when use_asc = FALSE, $J-1$ without an outside option, and $J$ with an outside option. The low-level parser requires at least one fixed-coefficient column ($K_x>0$), even when the substantive focus is entirely on random coefficients.

5.2 Delta Method for Variance Parameters

When reporting results, the Cholesky parameters $\ell$ are transformed to variance matrix elements $\Sigma = LL'$.

The Jacobian $J = \frac{\partial \text{vech}(\Sigma)}{\partial \ell}$ is computed analytically in jacobian_vech_Sigma().

For element $\Sigma_{ij} = \sum_k L_{ik} L_{jk}$: $$\frac{\partial \Sigma_{ij}}{\partial \ell_{pq}} = \frac{\partial L_{pq}}{\partial \ell_{pq}} \left( L_{jq} \mathbf{1}_{i=p} + L_{iq} \mathbf{1}_{j=p} \right)$$

This follows from the product rule applied to $\Sigma = LL^T$: $\partial \Sigma / \partial \ell_{pq} = E \cdot L^T + L \cdot E^T$, where $E$ is the matrix with a single non-zero entry $E_{pq} = \partial L_{pq} / \partial \ell_{pq}$.

The variance-covariance matrix of $\text{vech}(\Sigma)$ is then: $$V_\Sigma = J \cdot V_\ell \cdot J'$$

Code reference: src/mxlogit.cpp

5.3 OpenMP Parallelization

The implementation parallelizes over likelihood units (§1.6) using OpenMP, with dynamic scheduling for load balancing (#pragma omp for schedule(dynamic)): - The log-likelihood and gradient keep thread-local accumulators, combined in a #pragma omp critical section. - The Hessian and BHHH kernels add each unit’s compact block, and the cluster meat each cluster’s, into per-thread packed upper triangles, added up in thread order after the parallel region, or into a shared result when those would exceed 2 GiB (§4.4.6). - Before allocating, each estimation kernel compares an upper estimate of its working memory, per thread and shared, with the machine’s physical memory (or options(choicer.max_memory)), and stops with an error naming the number of threads that would fit instead of running out of memory inside the parallel region.

5.4 Practical Notes

1. Negated Objectives: mxl_loglik_gradient_parallel returns $-\ell(\theta)$ and $-\nabla\ell(\theta)$; mxl_hessian_parallel returns $-\nabla^2\ell(\theta)$. The BHHH kernel instead returns the positive-semidefinite information approximation $\sum_iw_is_is_i^T$.

2. ASC Identification: When include_outside_option = FALSE, the first inside alternative’s ASC ($\delta_1$) is fixed to zero for identification. When include_outside_option = TRUE, all inside ASCs are free parameters (the outside option with $V=0$ serves as the reference).

3. Gradient vectorization: The per-draw gradient is computed without a per-alternative loop. After forming $\mathbf{d}_i^s = \mathbf{e}_{j_i} - P_i^s$, the $\beta$ and $\mu$/$L$ blocks are evaluated with two BLAS dgemv calls ($X_i^T \mathbf{d}_i^s$ and $W_i^T \mathbf{d}_i^s$). The delta block uses an irregular scatter and remains a loop.

4. Exact draw layout and source: In both modes, likelihood unit $u$ (§1.6: a decision maker under person_col, a choice situation otherwise) receives the contiguous global indices $g=(u-1)S+s$ for $s=1,\ldots,S$; dimension $k$ uses the $k$th prime base and the resulting uniform is mapped through the standard-normal inverse CDF. In draws = "store", randtoolbox::halton(n = U*S, dim = K_w, normal = TRUE) is materialized as a $K_w\times S\times U$ cube during estimation, one $K_w \times S$ block per unit ($U$ decision makers under person_col, or $N$ choice situations otherwise). It is filled a block of units at a time with halton(start = ), which returns exactly those points. Store mode needs $US \le 2^{31} - 1$, the largest starting index halton() accepts. The prediction family (predict(), elasticities(), diversion_ratios(), blp(), logsum(), consumer_surplus(), and gof() through predict()) integrates over one block per choice situation instead ($N$ blocks for a panel fit, not $U$), which it forms on the fly from the same points of the sequence (§5.5), with no limit on $NS$ for up to 128 random coefficients. The fitted object retains only draws_info metadata, not the cube; the estimation-type functions (vcov(type = ), wesml_vcov(), conditional_tastes()) rebuild it. This mode is deterministic and has no seed or scramble. In draws = "generate", HaltonGen computes the radical inverse on demand. scramble = "none" uses identity permutations and is tested against the randtoolbox sequence. scramble = "permuted" (the default) constructs one independently seeded Fisher–Yates permutation for every (dimension, digit-position) pair and shares that table across sequence indices. Evaluation stops after the last nonzero base-$b$ digit, so trailing zero digits remain unchanged. This is not Owen nested-uniform scrambling: the digit permutation does not depend on the preceding digit prefix, and the construction does not make each point marginally uniform. Consequently it carries no standard RQMC unbiasedness, variance-rate, or replicate-based error-estimation guarantee. The historical API value "owen" is a deprecated alias for "permuted"; old fitted objects carrying that serialized label remain readable. The result is deterministic conditional on the seed and independent of OpenMP scheduling. Numerical sensitivity should be assessed by increasing $S$ and, for permuted draws, varying the seed. Generate mode is available only through the convenience data workflow; the advanced input_data workflow requires an explicit eta_draws cube.

5. Finite-value guards: The likelihood/gradient kernel replaces a non-finite global negated objective by $10^{10}$ and a zero gradient, allowing the optimizer’s line search to backtrack; isolated non-finite gradient entries are zeroed. The Hessian contributes a likelihood unit — decision maker $n$ under person_col, choice situation $i$ otherwise (§1.6, §4.5) — only when its log-sum-exp $\mathrm{lse}$ is finite. Since every $\log P_{ntj}^s$ is the chosen utility minus a log-sum-exp (§2.4), $\lambda_{ns}$, $\mathrm{lse}$ and the objective are finite whenever the utilities are, however poorly the parameter vector fits the observed choices: the $10^{10}$ sentinel, the Hessian’s skip and the NA conditional tastes (§13) all fire only when the utilities themselves overflow (for instance an exploding Cholesky diagonal during a line search), where the log-likelihood is undefined. (Earlier releases skipped a situation whenever $\widehat P_i^{\mathrm{sum}} = \sum_s P_{ij_i}^s \le 10^{-12}$, which could also discard a situation that was merely very poorly fit rather than genuinely zero-probability, biasing the reported Hessian and its standard errors; and until the log-space form of §2.4, a single choice probability that underflowed to zero, at a utility gap beyond about 745, made the whole objective the sentinel with a zero gradient, while subnormal probabilities, at gaps beyond about 708, lost precision. See NEWS.md.) Relatedly, the posterior draw weights $\omega_{ns}$ (§2.6) that both the score (§3.5) and the Hessian (§4.5) are built from are computed in log space via log-sum-exp rather than as a ratio of raw probabilities, so they stay finite and accurate even when the product over situations falls far below the smallest double — the earlier raw-probability construction could return an identically zero gradient at a genuinely poorly-fit point, a false stationary point for the optimizer. These are numerical failure guards, not changes to the mathematical criterion in its regular finite region.

5.5 Data Layout

• $X$ has $\sum_iM_i$ rows; situation $i$ occupies rows $\sum_{i'<i}M_{i'}+1,\ldots,\sum_{i'\le i}M_{i'}$ (64-bit offsets in the kernels).
• $W$ may be row-aligned with $X$ or have one row per global inside alternative; the kernels either read the situation’s block of rows or gather rows by the global 1-based alt_idx. The convenience R workflow constructs the row-aligned form.
• choice_idx is a 1-based local row position within each situation block, or 0 for an implicit outside choice. alt_idx is instead the global alternative ID used to gather ASCs and alternative-level $W$ rows.
• Stored draws have dimensions $K_w\times S\times U$, one block per likelihood unit ($U$ decision makers under person_col, choice situations otherwise). Generate-mode kernels receive an empty placeholder cube plus gen_seed, gen_scramble, and gen_S.
• Prediction kernels read one block per choice situation: situation $i$ gets points $(i-1)S+1,\ldots,iS$ of the Halton sequence. With stored draws the prediction-type functions (predict(), logsum(), elasticities(), diversion_ratios(), blp(), and through them consumer_surplus() and gof()) never build that $K_w\times S\times N$ cube (though blp() may keep the draws across its iterations, §8.4): they pass an empty placeholder and gen_scramble = 2, and each thread forms its situation’s block on the fly from HaltonGen’s identity-permutation uniforms, which are randtoolbox::halton()’s points, mapped to normals by R’s qnorm() (R::qnorm, libR’s own qnorm5), as get_halton_normals() maps them. The draws, and so the results, are the cube’s bit for bit (per situation at any thread count, and sums over situations at one thread; at more threads their rounding depends on the order in which threads take the situations, as it always did) where randtoolbox and choicer are compiled with the same floating-point contraction: both add each digit’s term $d_j b^{-(j+1)}$ to the radical inverse in the same order, and agree bit for bit when both fuse those multiply-adds or neither does, as on CRAN’s platforms. Otherwise under 1% of the uniforms, all of bases 5 and up, differ in their last bit. The on-the-fly generator covers up to 128 random coefficients; beyond them the functions build the cube.

---

6. Elasticity Computation

6.1 Elasticity Definitions

The elasticity measures the percentage change in choice probability with respect to a percentage change in an attribute. For the mixed logit model, we compute elasticities using simulated probabilities.

Simulated Probability:

$$\bar{P}_{ij} = \frac{1}{S} \sum_{s=1}^{S} P_{ij}^s$$

6.2 Elasticity for Fixed Coefficients (X Variables)

For a variable $k$ in the design matrix $X$ with fixed coefficient $\beta_k$:

Own-Elasticity (elasticity of $P_{ij}$ with respect to $x_{ijk}$):

$$E_{jj}^k = \frac{\partial \bar{P}_{ij}}{\partial x_{ijk}} \cdot \frac{x_{ijk}}{\bar{P}_{ij}}$$

Using the chain rule and averaging over draws:

$$E_{jj}^k = \frac{1}{\bar{P}_{ij}} \cdot \frac{1}{S} \sum_{s=1}^{S} \beta_k \cdot x_{ijk} \cdot P_{ij}^s \cdot (1 - P_{ij}^s)$$

Cross-Elasticity (elasticity of $P_{ij}$ with respect to $x_{imk}$ where $m \neq j$):

$$E_{jm}^k = \frac{1}{\bar{P}_{ij}} \cdot \frac{1}{S} \sum_{s=1}^{S} \left( -\beta_k \cdot x_{imk} \cdot P_{ij}^s \cdot P_{im}^s \right)$$

6.3 Elasticity for Random Coefficients (W Variables)

For a variable $k$ in the design matrix $W$ with random coefficient, the effective coefficient for individual $i$ and draw $s$ is:

$$\beta_{ik}^s = \mu_k^* + \gamma_{ik}^{s*}$$

where: - For normal distribution: $\mu_k^* = \mu_k$, $\gamma_{ik}^{s*} = \gamma_{ik}^s$ - For log-normal distribution: $\mu_k^* = \exp(\mu_k)$, $\gamma_{ik}^{s*} = \exp(\gamma_{ik}^s)$

Own-Elasticity:

$$E_{jj}^k = \frac{1}{\bar{P}_{ij}} \cdot \frac{1}{S} \sum_{s=1}^{S} \beta_{ik}^s \cdot w_{ijk} \cdot P_{ij}^s \cdot (1 - P_{ij}^s)$$

Cross-Elasticity:

$$E_{jm}^k = \frac{1}{\bar{P}_{ij}} \cdot \frac{1}{S} \sum_{s=1}^{S} \left( -\beta_{ik}^s \cdot w_{imk} \cdot P_{ij}^s \cdot P_{im}^s \right)$$

6.4 Weighted Average of Individual Elasticities

The implementation computes weighted average individual probability elasticities:

$$\bar{E}_{jm}^k = \frac{\sum_{i=1}^{N} w_i \cdot E_{ijm}^k}{\sum_{i=1}^{N} w_i}$$

This is not generally the elasticity of an aggregate share under a common market-level perturbation. Rows index responding alternatives and columns perturbed alternatives. The kernel sets a contribution to zero when $\bar P_{ij}\leq10^{-12}$ and divides the accumulated matrix by total weight over all situations; with ordinary finite utilities, the threshold is only a numerical guard.

Code reference: mxlogit.cpp:mxl_elasticities_parallel

---

7. Diversion Ratios

7.1 Definition

The attribute-based diversion ratio from alternative $j$ to alternative $m$ with respect to attribute $k$ is

$$\mathrm{DR}_{j \to m}(x_k) \;=\; -\,\frac{\partial s_m / \partial x_{j,k}}{\partial s_j / \partial x_{j,k}},$$

where $s_j = \frac{1}{\sum_i w_i}\sum_i w_i \bar{P}_{ij}$ is the aggregate share. It is the fraction of demand lost by $j$ that is captured by $m$ when a marginal change in $j$’s attribute $k$ reduces $s_j$.

7.2 MNL Special Case (β cancels)

For multinomial logit, the per-individual derivatives of the choice probability are

$$\frac{\partial P_{ij}}{\partial x_{j,k}} = \beta_k\, P_{ij}(1 - P_{ij}), \qquad \frac{\partial P_{im}}{\partial x_{j,k}} = -\beta_k\, P_{ij}P_{im} \quad (m \neq j).$$

Letting $\bar{w} = \sum_i w_i$, the aggregate share derivatives are

$$\frac{\partial s_j}{\partial x_{j,k}} \;=\; \frac{1}{\bar{w}}\sum_i w_i \frac{\partial P_{ij}}{\partial x_{j,k}} \;=\; \beta_k \cdot \frac{1}{\bar{w}}\sum_i w_i\, P_{ij}(1 - P_{ij}),$$

$$\frac{\partial s_m}{\partial x_{j,k}} \;=\; -\,\beta_k \cdot \frac{1}{\bar{w}}\sum_i w_i\, P_{ij} P_{im}.$$

Because $\beta_k$ is a scalar constant, it pulls out of both weighted sums. Forming the ratio,

$$\mathrm{DR}_{j \to m}^{\text{MNL}} \;=\; -\,\frac{\partial s_m / \partial x_{j,k}}{\partial s_j / \partial x_{j,k}} \;=\; \frac{\beta_k \cdot \tfrac{1}{\bar{w}}\sum_i w_i\, P_{ij} P_{im}}{\beta_k \cdot \tfrac{1}{\bar{w}}\sum_i w_i\, P_{ij}(1 - P_{ij})} \;=\; \frac{\sum_i w_i\, P_{ij} P_{im}}{\sum_i w_i\, P_{ij}(1 - P_{ij})},$$

with $\beta_k$ and $1/\bar{w}$ cancelling. The result is independent of which attribute is perturbed — this is a special property of MNL.

7.3 MXL: β Does Not Cancel

For mixed logit the realized coefficient on a random variable for individual $i$ at draw $s$ is

$$\beta_{ik}^s = \mu_k^* + \gamma_{ik}^{s*},$$

with the same distribution transforms as in §6.3. The derivatives evaluated at draw $s$ are

$$\frac{\partial P_{ij}^s}{\partial x_{j,k}} = \beta_{ik}^s\, P_{ij}^s\,(1 - P_{ij}^s), \qquad \frac{\partial P_{im}^s}{\partial x_{j,k}} = -\beta_{ik}^s\, P_{ij}^s P_{im}^s.$$

Integrating over draws gives $\partial \bar{P}_{ij}/\partial x_{j,k} = \frac{1}{S}\sum_s \beta_{ik}^s P_{ij}^s(1-P_{ij}^s)$ and similarly for the cross term. Aggregating over individuals, the diversion ratio is

$$\mathrm{DR}_{j \to m}^{\text{MXL}}(x_k) \;=\; \frac{\sum_i w_i \cdot \tfrac{1}{S}\sum_s \beta_{ik}^s\, P_{ij}^s P_{im}^s}{\sum_i w_i \cdot \tfrac{1}{S}\sum_s \beta_{ik}^s\, P_{ij}^s (1 - P_{ij}^s)}.$$

Because $\beta_{ik}^s$ varies across $(i, s)$, it cannot be pulled out of the sums and does not cancel. The MXL diversion ratio therefore depends on the perturbed variable. For a variable with a fixed coefficient the dependence again vanishes ($\beta_{ik}^s$ is a constant scalar); for a random-coefficient variable it does not.

7.4 Properties

Column-sum identity. For each $j$, $\sum_{m \neq j} \mathrm{DR}_{j \to m} = 1$. This follows from $\sum_{m \neq j} P_{ij}^s P_{im}^s = P_{ij}^s(1 - P_{ij}^s)$, which holds inside the draw sum and is preserved through the same $\beta_{ik}^s$ weight on both sides.

Sign of $\beta_{ik}^s$. If $\beta_{ik}^s$ is constant negative (e.g. a fixed price coefficient), both numerator and denominator change sign together and the ratio is unchanged — DR is non-negative. With a random coefficient of mixed sign across $(i, s)$, terms can partially cancel; the aggregate ratio is still well-defined but can be numerically sensitive. The implementation guards against this with an abs(denominator) > 1e-15 check.

Implementation. Cross-products $P_{ij}^s P_{im}^s$ and $P_{ij}^s(1 - P_{ij}^s)$ must be accumulated inside the per-draw loop, multiplied by $\beta_{ik}^s$ at that point, and only then averaged across draws and weighted by $w_i$. Computing $(\overline{P}_{ij})(\overline{P}_{im})$ or $\beta_{ik}^s$-weighted averages computed outside the draw loop is biased.

Code reference: mxlogit.cpp:mxl_diversion_ratios_parallel

---

8. BLP Contraction Mapping

8.1 Problem Statement

Given observed market shares $s_j$, find the ASC parameters $\delta_j$ such that the model-predicted shares match the observed shares:

$$\hat{s}_j(\delta) = s_j \quad \forall j$$

where $\hat{s}_j(\delta)$ is the predicted market share for alternative $j$ computed using simulated probabilities.

8.2 Simulated Market Shares

The predicted market share for alternative $j$ is:

$$\hat{s}_j(\delta) = \frac{1}{\sum_i w_i} \sum_{i=1}^{N} w_i \cdot \bar{P}_{ij}$$

where:

$$\bar{P}_{ij} = \frac{1}{S} \sum_{s=1}^{S} P_{ij}^s(\delta)$$

and $P_{ij}^s(\delta)$ is the probability for individual $i$, alternative $j$, and draw $s$, which depends on $\delta$ through the utility specification.

8.3 Contraction Mapping Algorithm

The Berry–Levinsohn–Pakes fixed-point iteration is:

$$\delta_j^{(t+1)} = \delta_j^{(t)} + \log(s_j) - \log(\hat{s}_j^{(t)})$$

This can be written in vector form as:

$$\delta^{(t+1)} = \delta^{(t)} + \log(s) - \log(\hat{s}^{(t)})$$

The target must be feasible and strictly positive. The kernel checks positivity of the inside target and predicted shares but does not verify that the full target sums to one.

8.4 Implementation Details

1. No outside option. delta may be the full length-$J$ vector or the free length-$(J-1)$ block; the latter is padded with a leading zero. The first element is subtracted initially and after every update, and the returned vector has length $J$ with its baseline entry equal to zero.
2. Outside option. delta has length $J$ and contains inside ASCs only. The outside utility remains zero, the update uses only the inside components of target and predicted shares, and no inside ASC is re-centered. A coherent target outside share is matched only indirectly because all shares sum to one.
3. Stopping. The residual is the largest absolute inside-delta change. On reaching max_iter, the kernel prints a warning and returns the last completed iterate.

The public blp.choicer_mxl() method resolves the fitted draw metadata through .mxl_pred_draws(): store-mode fits have the kernel form the plain-Halton points on the fly (§5.5), while generate-mode fits pass the original seed, permutation mode, and $S$ to the contraction kernel. Share inversion therefore draws from the same deterministic Halton sequence as the fit, one block per choice situation (for a cross-sectional fit, the estimation draws themselves; a panel fit’s estimation draws are per decision maker instead, §5.5). Neither the draws nor $L$ change across iterations, so blp() (through the internal mxl_blp_contraction_cached()) keeps each situation’s $\Gamma_i = L\eta_i$, with its log-normal rows exponentiated, from the first evaluation of the shares on, when the draws are formed on the fly (stored draws with up to 128 random coefficients, or generated ones), the $N$ blocks take at most keep_draws_bytes ($8K_wSN$ bytes; 4 GiB by default) and the allocation succeeds; otherwise every evaluation forms them again. The kept blocks are the values each evaluation would compute, so the shares, and the returned $\delta$, are the same either way. The exported mxl_blp_contraction() keeps nothing.

Code reference: mxlogit.cpp:mxl_blp_contraction

---

9. Choice-Based Sampling and WESML Weighting

9.1 Endogenous Stratified (Choice-Based) Sampling

In random (exogenous) sampling each choice situation is drawn independently of its outcome. Under choice-based (endogenous stratified) sampling the strata are defined by the chosen alternative $j_i$, and situations are drawn with stratum-specific frequencies – for example, deliberately oversampling the few individuals who chose a rare alternative. Let

$$Q(j) = \text{population share of alternative } j, \qquad H(j) = \text{sample share of choosers of } j .$$

Maximum simulated likelihood applied naively to such a sample is generally inconsistent for the population parameters. (A known special case: in a pure multinomial logit with a full set of alternative-specific constants, choice-based sampling biases only the constants; in general – and for the mixed logit – both the slopes and the constants are affected.)

9.2 The WESML Weighted Log-Likelihood

Manski and Lerman (1977) restore consistency by weighting each situation’s contribution by the ratio of population to sample share of its chosen alternative,

$$w_i = \frac{Q(j_i)}{H(j_i)},$$

which yields the weighted simulated log-likelihood already defined in $$2.3,

$$\ell^{W}(\theta) = \sum_{i=1}^{N} w_i \log \bar{P}_{i j_i}(\theta).$$

The maximizer $\hat{\theta}$ is invariant to multiplying every $w_i$ by a common positive constant (the objective is merely rescaled), so the weights may be normalized to mean 1 without changing the estimates.

9.3 Robust (Sandwich) Variance: the $w^2$ Meat

WESML is a weighted M-estimator solving $\sum_i w_i\, s_i(\hat{\theta}) = 0$, where $s_i = \partial \log \bar{P}_{i j_i} / \partial\theta$ is the per-situation score ($$3) and $H_i = \partial^2 \log\bar{P}_{i j_i}/\partial\theta\,\partial\theta^\top$ its Hessian ($$4). Its robust (Huber–White) asymptotic variance is the sandwich

$$V = A^{-1} B A^{-1}, \qquad A = \sum_{i} w_i\,(-H_i), \qquad B = \sum_{i} w_i^{2}\, s_i s_i^{\top}.$$

The weight enters the bread $A$ linearly, but the meat $B$ carries the weight squared: the contribution of situation $i$ to the estimating equation is $w_i s_i$, so the variance of that contribution involves $(w_i s_i)(w_i s_i)^{\top} = w_i^{2}\, s_i s_i^{\top}$.

This is exactly why the two naive variances are wrong under weighting:

• the inverse-Hessian $A^{-1}$ assumes the information-matrix equality $A = B$, which fails once the $w_i$ are non-degenerate;
• the ordinary BHHH/OPG $\big(\sum_i w_i\, s_i s_i^{\top}\big)^{-1}$ uses the weight to the first power, not the second.

A consistency check: rescaling all weights by a constant $c$ sends $A \to cA$ and $B \to c^{2} B$, so $V \to (cA)^{-1}(c^{2}B)(cA)^{-1} = A^{-1}BA^{-1}$ – the sandwich is invariant to the weight scale, whereas $A^{-1} \to c^{-1} A^{-1}$ is not.

9.4 Scope: Robust vs. Design-Based Variance

The estimator implemented here is the robust weighted-M-estimator (Huber–White) variance, whose meat $B = \sum_i w_i^{2}\, s_i s_i^{\top}$ is uncentered. This is the asymptotic variance under variable-probability sampling (each population unit sampled independently with a stratum-dependent probability). Under a fixed-quota design – a prescribed number drawn from each choice stratum, i.e. standard stratified sampling – the design-based variance centers the scores within strata and may be smaller; that stratum-centered (design-based) estimator is out of scope here. See Manski & McFadden (1981) for the design-based treatment and Cosslett (1981) for the asymptotically efficient conditional-maximum-likelihood alternative.

9.5 Implementation and Cluster-Robust Extension

The variance routes never form the $U\times n_\text{params}$ matrix $S$ of the units’ weight-free scores: two million decision makers by 2,620 parameters would take 42 GB. With the bread $A=\texttt{mxl\_hessian\_parallel}(\widehat\theta;w)$, the robust meat is the BHHH kernel with squared weights, summed over each unit’s compact block (§4.4.6),

$$B=\sum_i w_i^2 s_i s_i^T=\texttt{mxl\_bhhh\_parallel}(\widehat\theta;w^2),$$

and vcov(type = "bhhh") is the same kernel with weights $w$. The matrices are combined as $V = A^{-1} B A^{-1}$, and one set of draws, stored or generated from the fit’s draw metadata, serves both. For clusters $g(i)$, the meat

$$B_{\mathrm{cluster}}=\sum_g \left(\sum_{i:g(i)=g}w_is_i\right) \left(\sum_{i:g(i)=g}w_is_i\right)^T$$

is mxl_cluster_meat_parallel: the units are sorted by cluster and cut into chunks of similar numbers of rows, about four per thread; each cluster’s weighted scores are summed over the parameters they touch, and a cluster that a chunk boundary cuts is summed in parts, added after the parallel region, so a few large clusters still keep every thread busy. There is no finite-cluster correction. mxl_scores_parallel (internal; no variance route of the mixed logit calls it) still returns $S$ for small problems and refuses a matrix above 2 GiB. Clustering changes inference, not the estimand of whichever likelihood was fit (§1.6): cross-sectional by default, or the panel likelihood of §2.6 under person_col, where cluster labels must nest decision makers. It is only §9.1–§9.2’s WESML weighting that is defined solely for the cross-sectional likelihood: $w_i$ is a weight that varies with the alternative chosen in situation $i$, whereas the panel objective of §2.6 requires one weight $w_n$ per decision maker, constant across their situations, so a choice-based weight cannot enter it — run_mxlogit(..., person_col = ) therefore rejects choice-based-sampling provenance. The sandwich variance of §9.3 and its cluster-robust extension here in §9.5 carry no such restriction: both are built from whichever likelihood unit’s score was actually fit (a decision maker under person_col, a choice situation otherwise), so they apply to panel fits too, with cluster labels then nesting decision makers exactly as above. A scaled fit’s variance is computed in natural units and inverted equilibrated ($$1.5); the Cholesky-to-$$ delta method ($$5.2) applies on top, exactly as for any stored vcov. User-facing entry points: wesml_weights() and sample_by_choice() build the weights; run_mxlogit(..., se_method = "sandwich") estimates with the robust variance; wesml_vcov() returns it post hoc.

Code reference: mxl_bhhh_parallel(), mxl_cluster_meat_parallel() and mxl_scores_parallel() in src/mxlogit.cpp; compute_sandwich_vcov(), .assemble_score_vcov(), .mxl_score_meat() and (MNL, NL) .score_meat() in R/classes.R; wesml_weights() / sample_by_choice() in R/sampling.R.

---

10. Willingness to Pay

10.1 Fixed Coefficients

When the price variable has a fixed coefficient $\alpha = \beta_p$ (a column of $X$), the WTP for any fixed attribute $k$ is the usual marginal rate of substitution, with the same analytic delta-method SE as in the MNL case (§8 of the MNL document):

$$\mathrm{WTP}_k = -\frac{\beta_k}{\alpha}, \qquad \nabla g = \begin{pmatrix} -1/\alpha \\ \beta_k/\alpha^2 \end{pmatrix}, \qquad \widehat{\mathrm{SE}} = \sqrt{\nabla g^T \hat{V}_{(k,p)} \nabla g}$$

10.2 Random Coefficients

For a random attribute coefficient $\beta_k(\eta)$ the WTP is itself a random variable $-\beta_k(\eta)/\alpha$, and a scalar summary must be chosen. The summaries reported follow the package’s parameterization (§1.3):

Normal RC with estimated mean (rc_mean = TRUE, rc_dist[k] = 0): $\beta_k = \mu_k + (L\eta)_k$ has mean $\mu_k$, so the mean WTP is

$$\mathrm{WTP}_k = -\frac{\mu_k}{\alpha}$$

with the linear-ratio gradient of §10.1 applied to $(\mu_k, \alpha)$.

Shifted log-normal RC with estimated location (rc_mean = TRUE, rc_dist[k] = 1): $\beta_k = \exp(\mu_k) + \exp((L\eta)_k)$. Since $\beta_k$ is monotone in the normal draw $(L\eta)_k$ and the median of $\exp((L\eta)_k)$ is $\exp(0) = 1$, the population median is $\exp(\mu_k) + 1$ and the median WTP is

$$\mathrm{WTP}_k^{\text{med}} = -\frac{\exp(\mu_k) + 1}{\alpha}, \qquad \frac{\partial g}{\partial \mu_k} = -\frac{\exp(\mu_k)}{\alpha}, \qquad \frac{\partial g}{\partial \alpha} = \frac{\exp(\mu_k) + 1}{\alpha^2}$$

The median is preferred to the mean, $\exp(\mu_k) + \exp(\sigma_k^2/2)$ with $\sigma_k^2 = (LL^T)_{kk}$, because the mean is highly sensitive to the estimated variance (a thin right tail dominates it).

Shifted log-normal RC without location (rc_mean = FALSE, rc_dist[k] = 1): $\beta_k = \exp((L\eta)_k)$ has median exactly 1, so the median WTP is $-1/\alpha$ with $\partial g / \partial \alpha = 1/\alpha^2$ — uncertainty comes solely from $\hat\alpha$.

Normal RC without mean (rc_mean = FALSE, rc_dist[k] = 0): the mean is 0 by construction; no WTP row is reported.

10.3 Random Price Coefficients Are Rejected

The implementation requires price_var to be a fixed-coefficient column of $X$ and rejects every random price coefficient. For a normally distributed denominator this is substantively important: positive density arbitrarily close to zero generally makes $E|\beta_k/\alpha|$ diverge. That no-moment argument does not cover every possible mixing distribution—for example, inverse moments of a sign-constrained log-normal can be finite—but wtp() still does not simulate the joint ratio distribution, account for its correlation, or report its quantiles. Thus the rejection is both a specification guard for normal price mixing and an implementation boundary for other distributions. WTP-space estimation (Train and Weeks 2005) is not implemented in choicer v0.2.0.

Code reference: R/wtp.R

---

11. Consumer Surplus and the Simulated Logsum

11.1 Expected Logsum

Conditional on a coefficient realization $\beta(\eta)$, expected maximum utility over the iid type-I extreme-value errors is $\log \sum_j \exp(V_{ij}(\eta))+\gamma_E$, where $\gamma_E$ is Euler’s constant. As in the other package methods, logsum() omits this scenario-invariant additive constant. The mixed-logit logsum integrates over the mixing distribution:

$$\mathrm{logsum}_i = \mathbb{E}_\eta\left[\log \sum_{j \in C_i} \exp\big(V_{ij}(\eta)\big)\right] \;\approx\; \frac{1}{S} \sum_{s=1}^{S} \log \sum_{j \in C_i} \exp\big(V_{ij}^s\big)$$

simulated over the fit’s Halton sequence, one block per choice situation (regenerated from the stored draw metadata; for a cross-sectional fit, the estimation draws themselves; a panel fit’s estimation draws are per decision maker instead), with the outside option contributing $\exp(0)$ in each draw’s sum when present, and the max-subtraction trick applied per draw.

Order of operations matters. The log-sum-exp function is convex, so by Jensen’s inequality

$$\frac{1}{S}\sum_s \mathrm{lse}(V^s) \;\geq\; \mathrm{lse}\left(\frac{1}{S}\sum_s V^s\right),$$

with strict inequality when the draws induce non-common variation in utilities. A nondegenerate $\Sigma$ alone is not sufficient for strictness: if the random contribution is the same additive shift for every available alternative, log-sum-exp is affine in that direction and equality holds. A dedicated C++ kernel (mxl_logsum) computes the per-draw logsum and then averages across draws.

11.2 Expected Consumer Surplus

With a fixed price coefficient $\alpha$ (the marginal utility of income $-\alpha$ is then nonrandom):

$$\mathbb{E}[CS_i] = \frac{\mathrm{logsum}_i}{-\alpha}$$

A random price coefficient is rejected for the capability and specification reasons in §10.3; the divergence argument applies directly to normal price mixing but not to every sign-constrained distribution. As in MNL, the formula assumes no income effects. The routine does not enforce $\alpha<0$; the economic interpretation as positive marginal utility of income requires the analyst to verify that sign. CS levels inherit the ASC normalization, so only differences across counterfactual scenarios are meaningful. The implementation reports point estimates for MXL (the delta-method SE of mean CS is currently provided for MNL only); Krinsky–Robb resampling is a practical interval method.

Code reference: R/surplus.R; kernel in src/mxlogit.cpp (mxl_logsum)

---

12. Goodness of Fit

The measures are model-agnostic and identical to §10 of the MNL document, computed from the simulated log-likelihood and simulated choice probabilities:

• McFadden pseudo R²: $R^2 = 1 - \ell(\hat\theta)/\ell_0$, adjusted $R^2_{\text{adj}} = 1 - (\ell(\hat\theta) - K)/\ell_0$, with $K$ counting all estimated parameters ($\beta$, $\mu$, Cholesky block, ASCs).
• Equal-shares null (default; exact for unbalanced choice sets and weights): $\ell_0 = -\sum_i w_i \log(M_i + \mathbf{1}_{\text{outside}})$.
• Market-shares null (constants-only closed form): $\ell_0 = c\sum_j N_j \log(s_j)$, where $c$ is the common weight, valid only for identical choice-set composition and uniform weights. The constants-only model contains no random coefficients, so the MNL closed form applies unchanged.
• Hit rate: weighted share of choice situations whose observed choice has the highest simulated probability $\bar{P}_{ij} = \frac{1}{S}\sum_s P_{ij}^s$; with an outside option, the kernel returns the simulated outside probability directly and the outside good competes for the maximum.

Code reference: R/gof.R

---

13. Conditional (Individual-Level) Tastes

conditional_tastes() recovers, for each likelihood unit, the mean and standard deviation of the random coefficients conditional on the choices observed — the same posterior draw weights $\omega_{ns}$ that build the panel score and Hessian (§2.6, §3.5, §4.5), now applied to the realized taste draws themselves rather than to the score:

$$E[\beta_n \mid y_n, x_n; \hat\theta] \approx \sum_{s=1}^S \omega_{ns}\, \beta_{ns}, \qquad \mathrm{SD}(\beta_n \mid y_n, x_n; \hat\theta) \approx \sqrt{\sum_{s=1}^S \omega_{ns}\,(\beta_{ns} - \bar\beta_n)^2},$$

with $\beta_{ns} = \mu^*_\cdot + \gamma_{ns}^{*}$ the taste at draw $s$ on the scale it enters utility (§1.1–1.3), $\bar\beta_n$ the conditional mean above, and $y_n$ the sequence of choices observed across $n$’s $T_n$ situations (Revelt and Train 2000; Train 2009, ch. 11). $T_n \equiv 1$ defines the same object for the default cross-sectional fit, one choice situation at a time; a single choice reveals comparatively little, so its conditional mean sits closer to the population mean than a panel decision maker’s typically does.

Not $\beta_n$ itself. $E[\beta_n \mid y_n, x_n; \hat\theta]$ is the average taste among decision makers who would make $n$’s choices, not $n$’s own realized $\beta_n$ — it shrinks toward the population mean, more so the fewer situations $n$ contributes. Consequently the variance across decision makers of the conditional means understates $\Sigma$: by the law of total variance, $$\mathrm{Var}(\beta) = \underbrace{\mathrm{Var}_n\bigl(E[\beta_n \mid y_n]\bigr)}_{\texttt{revealed}} \;+\; \underbrace{E_n\bigl[\mathrm{Var}(\beta_n \mid y_n)\bigr]}_{\texttt{residual}},$$ so the revealed share ($w_n$-weighted variance of the conditional means) and the residual share ($w_n$-weighted average of the conditional variances) should sum to about one. These are numerical checks, not specification tests: for normal random coefficients with free means, both this identity and $\texttt{mean\_cond} \approx \texttt{model\_mean}$ are the likelihood’s first-order conditions for $\mu$ and $\Sigma$ — the population mean and variance are exactly the $\omega$-weighted mean and variance the estimating equations set to zero at $\hat\theta$ — so they hold at a converged maximum whatever the true taste distribution. A departure points to too few draws $S$ or an unconverged $\hat\theta$; when a coefficient’s mean is fixed at zero (rc_mean = FALSE and no fixed coefficient on the same variable), a mean_cond far from zero can also flag that restriction. For log-normal coefficients the variance identity holds on the underlying normal scale $\gamma$, not the utility scale $\beta$ reported here, so expect departures of a few percent in either direction at moderate $S$. print.choicer_tastes() reports these alongside model_mean/model_sd (the moments of $f(\beta \mid \hat\theta)$ itself, §1.2–1.3) and mean_cond/sd_cond (the weighted mean and SD of the conditional means across units).

Implementation. mxl_conditional_tastes_parallel() reuses mxl_unit_load() and mxl_unit_simulate() (§3.5) to obtain $\omega_{ns}$ per unit, then forms $\bar\gamma_n = \Gamma_n\,\omega_n$ over the unit’s $K_w \times S$ realized-deviation block $\Gamma_n$ (the same $\gamma_{ns}^{*}$ values already used to build $V_{ij}^s$, §1.1) and adds back the mean block, $\bar\beta_n = \mu^* + \bar\gamma_n$; the SD uses the analogous weighted second moment of $\Gamma_n$ around $\bar\gamma_n$. Both are matrix-vector products — no second pass is needed, since neither quantity is an outer product. Units whose log-likelihood is non-finite — possible only when their utilities overflow (§4.5, §5.4) — return NA and are dropped from the population summary.

Code reference: mxl_conditional_tastes_parallel() in src/mxlogit.cpp; R/conditional_tastes.R.

---

References

• Train, K. E. (2009). Discrete Choice Methods with Simulation. Cambridge University Press.
• Revelt, D., & Train, K. (1998). Mixed logit with repeated choices: households’ choices of appliance efficiency level. Review of Economics and Statistics, 80(4), 647-657.
• Revelt, D., & Train, K. (2000). Customer-specific taste parameters and mixed logit: households’ choice of electricity supplier. Working Paper E00-274, Department of Economics, University of California, Berkeley.
• Louis, T. A. (1982). Finding the observed information matrix when using the EM algorithm. Journal of the Royal Statistical Society: Series B, 44(2), 226-233.
• Train, K., & Weeks, M. (2005). Discrete choice models in preference space and willingness-to-pay space. In R. Scarpa & A. Alberini (Eds.), Applications of Simulation Methods in Environmental and Resource Economics (pp. 1-16). Springer.
• McFadden, D., & Train, K. (2000). Mixed MNL models for discrete response. Journal of Applied Econometrics, 15(5), 447-470.
• Berry, S., Levinsohn, J., & Pakes, A. (1995). Automobile prices in market equilibrium. Econometrica, 63(4), 841-890.
• Bhat, C. R. (2003). Simulation estimation of mixed discrete choice models using randomized and scrambled Halton sequences. Transportation Research Part B, 37(9), 837-855.
• Gouriéroux, C., & Monfort, A. (1996). Simulation-Based Econometric Methods. Oxford University Press.
• Hajivassiliou, V. A., & Ruud, P. A. (1994). Classical estimation methods for LDV models using simulation. In R. F. Engle & D. L. McFadden (Eds.), Handbook of Econometrics, Vol. 4 (pp. 2383-2441). North-Holland.
• Lee, L.-F. (1995). Asymptotic bias in simulated maximum likelihood estimation of discrete choice models. Econometric Theory, 11(3), 437-483.
• Manski, C. F., & Lerman, S. R. (1977). The estimation of choice probabilities from choice based samples. Econometrica, 45(8), 1977-1988.
• Manski, C. F., & McFadden, D. (1981). Alternative estimators and sample designs for discrete choice analysis. In Structural Analysis of Discrete Data with Econometric Applications. MIT Press.
• Cosslett, S. R. (1981). Maximum likelihood estimator for choice-based samples. Econometrica, 49(5), 1289-1316.
• Wooldridge, J. M. (2010). Econometric Analysis of Cross Section and Panel Data (2nd ed.), Section 13.8. MIT Press.