13  The Mixed MNL Model by MSL

This chapter assembles the summit model of the likelihood track. The mixed multinomial logit — MNL choice behavior with coefficients that vary continuously across people — is where Part II has been heading since heterogeneity entered in Chapter 11: the discrete classes melt into a continuous distribution, the random intercept of Chapter 12 grows into a random vector, and every tool we have built gets used at once. It is also, not coincidentally, the single most-used model in applied choice work, from transportation demand to virtually every commercial conjoint study run today (McFadden and Train 2000; K. E. Train 2009, Ch. 6).

By the end of the chapter the model will be simulated, estimated by MSL, validated against mlogit, and audited for simulation error — and we will have felt, concretely, where the approach strains. That strain is data for Chapter 14.

13.1 Random Coefficients, Two Ways to Read Them

The model adds one letter to the MNL. Utility is still linear with Gumbel errors, but the coefficient vector now belongs to the person:

\[ U_{ntj} = \mathbf{x}_{ntj}'\boldsymbol{\beta}_n + \varepsilon_{ntj}, \qquad \boldsymbol{\beta}_n \sim N(\bar{\boldsymbol{\beta}}, \Sigma), \qquad \varepsilon_{ntj} \stackrel{iid}{\sim}\text{Gumbel} \tag{13.1}\]

Each person draws their taste vector \(\boldsymbol{\beta}_n\) once from a multivariate normal population — mean \(\bar{\boldsymbol{\beta}}\), covariance \(\Sigma\) — and carries it through all their choice tasks. The estimation targets are the population parameters \(\boldsymbol{\theta}= (\bar{\boldsymbol{\beta}}, \Sigma)\): what tastes are like on average, how much they vary, and (through \(\Sigma\)’s off-diagonals) which tastes travel together.

Two readings of Equation 13.1 coexist, and both matter:

• The taste-heterogeneity reading. People genuinely differ — some care intensely about price, some about brand — and \(\Sigma\) is the market’s heterogeneity, the quantity segmentation and targeting monetize.
• The substitution-pattern reading. Integrate the person-level MNL over the population and the resulting aggregate choice probabilities no longer satisfy IIA: people who like Apples are the same people across tasks, so alternatives similar in attribute space draw from overlapping customer pools — correlation in utilities, achieved through coefficients rather than Chapter 10’s nests, and estimated rather than pre-specified. This is the data-driven substitution flexibility that Section 10.8 promised.

The second reading has a remarkable theoretical guarantee: McFadden and Train (2000) prove that a mixed logit with appropriately chosen variables and mixing distribution can approximate any random utility model’s choice probabilities arbitrarily well. Nested logit, probit, anything — all live in the mixed logit’s closure. The theorem is an existence result, not a recipe (it doesn’t say the approximating specification is parsimonious or findable), but it explains the model’s status: this is not one more family member; it is a universal chassis.

13.2 Specification Choices

Equation 13.1 as written leaves the modeler three decisions.

Which coefficients vary? Nothing requires all \(K\). Our specification — motivated the way an applied researcher would — lets tastes vary where taste variation is behaviorally plausible and commercially interesting: the Apple premium (Mac people versus the indifferent), the top-memory premium (power users versus browsers), and price sensitivity (budgets differ). The dell, ram16, and screen15 coefficients stay fixed.¹ In the notation below, rc indexes the random subset.

What mixing distribution? Normal, here. Its known weakness: it puts mass on both signs of every random coefficient, so our specification implies some fraction of consumers with positive price coefficients — about \(\Phi(-1.2/0.7) \approx 4\%\) prefer expensive, all else equal, under the truth below. For this chapter we accept the blemish (it is small, and normality keeps every layer transparent); the standard fix — a lognormal for the price coefficient, built from the Chapter 2 transformation — is deployed in Chapter 17, and more flexible mixing distributions are Chapter 22 material (K. Train 2016).

Correlated or independent? A full \(\Sigma\) has \(K_r(K_r+1)/2\) free elements for \(K_r\) random coefficients, kept positive-definite by optimizing over its Cholesky factor \(L\) (Section 2.4 — estimate \(L\)’s elements, log-parameterizing the diagonal, and \(\Sigma = LL'\) is valid by construction). We estimate the diagonal (independent) version: three standard deviations, log-parameterized as in Chapter 12. The reason is worth stating baldly now and remembering at the chapter’s end: with MSL, every step toward a full covariance matrix costs real optimization pain, and we will tally that cost rather than pay it — then watch Chapter 17 get the full matrix essentially for free.

The truth for our study:

\[ \bar{\boldsymbol{\beta}}^\ast = (0.5, \; 1.0, \; 0.6, \; 0.9, \; 0.3, \; -1.2), \qquad \sigma^\ast_{\texttt{apple}} = 0.8, \;\; \sigma^\ast_{\texttt{ram32}} = 0.6, \;\; \sigma^\ast_{\texttt{price}} = 0.7 \tag{13.2}\]

Means as always (Equation 7.6); spreads large enough that heterogeneity is the story, not a garnish — a price-coefficient standard deviation of \(0.7\) around \(-1.2\) makes the market’s price sensitivity range enormously.

13.3 The Likelihood

Assemble it by the now-familiar hierarchy logic — third time, so briskly. Conditional on \(\boldsymbol{\beta}_n\): the person’s sequence probability is a product of MNL terms over their tasks. Unconditional: integrate over the population,

\[ \Pr(\mathbf{y}_n \mid \boldsymbol{\theta}) = \int \left[ \prod_{t=1}^{T} P_{nt}(\boldsymbol{\beta}) \right] f(\boldsymbol{\beta}\mid \bar{\boldsymbol{\beta}}, \Sigma) \, d\boldsymbol{\beta} \tag{13.3}\]

— the integral outside the product (the person keeps their draw), no closed form (logit kernel against normal density, now in \(K_r\) dimensions), and therefore MSL (Chapter 12): fixed standard-normal draws \(\mathbf{z}_n^{(r)}\), mapped through candidate parameters, averaged on the log scale. The only novelty versus Chapter 12 is dimension — and dimension is precisely where simulation’s advantage over quadrature lives (Chapter 12’s footnote, now cashed).

One bookkeeping note: Equation 13.3 is written over the full coefficient vector, but with dell, ram16, and screen15 held fixed their “distribution” is a point mass — those dimensions integrate out trivially, and the integral we actually simulate is \(K_r = 3\)-dimensional, over the random subset only.

13.4 Simulating the Study

The simulator is Section 7.6’s with one insertion — person-level coefficient draws, via the location-scale map (diagonal case) that Chapter 2 built:

sim_mixed_data <- function(N, T, J, bar_beta, sd_beta, rc) {
    df <- sim_laptop_design(N, T, J)
    X  <- model.matrix(~ brand + ram + screen + price, data = df)[, -1]
    K  <- ncol(X)
    task_id <- cumsum(!duplicated(df[, c("id", "task")]))

    betas <- matrix(rep(bar_beta, each = N), N, K)      # start everyone at the mean
    for (k in rc) betas[, k] <- betas[, k] + sd_beta[k] * rnorm(N)

    V   <- rowSums(X * betas[df$id, ])                  # each row: its person's beta
    eps <- -log(-log(runif(nrow(df))))
    U   <- V + eps
    df$choice <- as.integer(ave(U, task_id, FUN = function(u) u == max(u)))
    attr(df, "betas") <- betas                          # answer key, person level
    df
}

set.seed(130)
bar_true <- c(dell = 0.5, apple = 1.0, ram16 = 0.6, ram32 = 0.9,
              screen15 = 0.3, price = -1.2)
sd_true  <- c(0, 0.8, 0, 0.6, 0, 0.7)
rc       <- c(2, 4, 6)                                  # apple, ram32, price

d <- sim_mixed_data(N = 500, T = 8, J = 3, bar_true, sd_true, rc)
X <- model.matrix(~ brand + ram + screen + price, data = d)[, -1]
y <- d$choice
task_id <- cumsum(!duplicated(d[, c("id", "task")]))
id <- d$id; n_id <- 500; K <- ncol(X)

The betas[df$id, ] indexing is Chapter 11’s row-specific-coefficients idiom with classes replaced by persons. And note we stored the true individual \(\boldsymbol{\beta}_n\)’s as an attribute — Chapter 14 will grade individual-level estimates against them.

A heterogeneity-specific sanity check joins the usual battery: within-person behavior should be more consistent than the population is. Cheap version — among people who ever bought the most expensive laptop in a task, what fraction of their other choices were also top-priced, versus the population rate?

chosen <- d[d$choice == 1, ]
top_price <- as.integer(ave(d$price, task_id, FUN = max)[d$choice == 1] == chosen$price)
person_rate <- tapply(top_price, chosen$id, mean)
c(population_rate = round(mean(top_price), 2),
  sd_across_people = round(sd(person_rate), 2))

 population_rate sd_across_people 
            0.21             0.15 

The across-person spread in that rate is far larger than binomial noise from eight tasks would produce — persistent individual differences, visible before any model is fit. (The formal version of this check is the LM test of McFadden and Train (2000); the informal version usually suffices to know heterogeneity is worth modeling.)

13.5 The Simulated Likelihood, With Its Gradient

Here is the estimation engine — the chapter’s centerpiece function. It computes the simulated log-likelihood of Equation 13.3 and, optionally, its analytic gradient; the parameter vector stacks the \(K\) means and the log standard deviations of the random subset:

msl_mixed <- function(theta, y, X, task_id, id, n_id, rc, Z,
                      want_grad = FALSE) {
    K   <- ncol(X)
    R   <- dim(Z)[2]
    bar <- theta[1:K]
    s   <- exp(theta[(K+1):(K+length(rc))])
    id_of_task <- id[!duplicated(task_id)]

    lseq <- matrix(0, n_id, R)          # log P(sequence_n | draw r)
    if (want_grad) {
        Gbar <- array(0, c(n_id, R, K))
        Gs   <- array(0, c(n_id, R, length(rc)))
    }
    for (r in 1:R) {
        betas <- matrix(rep(bar, each = n_id), n_id, K)
        for (jj in seq_along(rc))
            betas[, rc[jj]] <- betas[, rc[jj]] + s[jj] * Z[, r, jj]
        V    <- rowSums(X * betas[id, ])
        Vmax <- ave(V, task_id, FUN = max)
        eV   <- exp(V - Vmax)
        logP <- V - Vmax - log(rowsum(eV, task_id)[task_id])
        lseq[, r] <- as.vector(rowsum(logP[y == 1], id_of_task))
        if (want_grad) {
            p  <- eV / rowsum(eV, task_id)[task_id]
            Gp <- rowsum(X * (y - p), id)              # per-person MNL score
            Gbar[, r, ] <- Gp                       # d l_nr / d bar_beta
            for (jj in seq_along(rc))
                Gs[, r, jj] <- Gp[, rc[jj]] * Z[, r, jj] * s[jj]   # d l_nr / d log sigma_jj
        }
    }
    M  <- apply(lseq, 1, max)                          # log-mean-exp, per person
    W  <- exp(lseq - M)
    Sw <- rowSums(W)
    ll <- sum(M + log(Sw / R))
    if (!want_grad) return(ll)

    w    <- W / Sw                                     # person-by-draw weights
    gbar <- sapply(1:K, function(k) sum(w * Gbar[, , k]))
    gs   <- sapply(seq_along(rc), function(jj) sum(w * Gs[, , jj]))
    list(ll = ll, grad = c(gbar, gs))
}

Read it in layers, because every layer is a callback with a job:

• The draw loop (for (r in 1:R)) is Chapter 9’s tolerable loop — \(R\) iterations, each one fully vectorized over all 12,000 rows via the Chapter 8 kernel (max-shift, rowsum(), broadcast). Per-person sequence log-probabilities come from one more rowsum(..., id_of_task).
• The log-mean-exp across draws is Chapter 12’s device: sequence probabilities are underflow bait, so the average over draws happens under a per-person max-shift.
• The gradient is where the mathematics rewards attention. Write \(\ell_{nr} = \log \prod_t P_{nt}(\boldsymbol{\beta}_n^{(r)})\) for person \(n\)’s log sequence probability under draw \(r\) — the code’s lseq[n, r], a building block of the log-likelihood, not the log-likelihood \(\ell(\boldsymbol{\theta})\) itself. The simulated probability is \(\widehat{\Pr}(\mathbf{y}_n) = \frac{1}{R}\sum_r e^{\ell_{nr}}\), and differentiating its log is one chain-rule line: \[\frac{\partial}{\partial \boldsymbol{\theta}} \log \frac{1}{R}\sum_r e^{\ell_{nr}} = \frac{\sum_r e^{\ell_{nr}} \; \partial \ell_{nr}/\partial \boldsymbol{\theta}}{\sum_{r'} e^{\ell_{nr'}}}\] (the \(1/R\) cancels top and bottom) — a weighted average of per-draw gradients, each draw’s gradient counted by how well that draw explains the person:

\[ \frac{\partial}{\partial \boldsymbol{\theta}} \log \widehat{\Pr}(\mathbf{y}_n) = \sum_{r} w_{nr} \, \frac{\partial \ell_{nr}}{\partial \boldsymbol{\theta}}, \qquad w_{nr} = \frac{e^{\ell_{nr}}}{\sum_{r'} e^{\ell_{nr'}}} \tag{13.4}\]

And the per-draw gradient is an old friend: with respect to the means it is exactly the MNL score \(\mathbf{X}'(\mathbf{y}- \mathbf{p})\) of Equation 8.4, summed within person; with respect to \(\log\sigma_k\), the chain rule appends a factor \(z_{nk}^{(r)} \sigma_k\). Look hard at those weights \(w_{nr}\): they are the conditional distribution of “which draw is this person” given their choices — Bayes’ rule, materializing inside a frequentist gradient, and the draw-level twin of Chapter 11’s class weights \(w_{nc}\) (Equation 11.3): prior density times sequence likelihood, normalized, with “which class” replaced by “which draw.” Chapter 14 picks up precisely this object.

Per standing policy, the analytic gradient faces the numerical audit before we trust it:

set.seed(131)
R <- 50
Z <- array(0, c(n_id, R, length(rc)))
bases <- c(2, 3, 5)                     # one prime base per random dimension
for (jj in seq_along(rc)) {
    h <- halton(n_id*R + 50, bases[jj])[-(1:50)]
    Z[, , jj] <- qnorm(matrix((h + runif(1)) %% 1, n_id, R, byrow = TRUE))
}

theta_test <- c(bar_true, log(c(0.8, 0.6, 0.7)))
g <- msl_mixed(theta_test, y, X, task_id, id, n_id, rc, Z, want_grad = TRUE)
num <- sapply(1:9, function(k) {
    e <- replace(numeric(9), k, 1e-5)
    (msl_mixed(theta_test + e, y, X, task_id, id, n_id, rc, Z) -
     msl_mixed(theta_test - e, y, X, task_id, id, n_id, rc, Z)) / 2e-5
})
round(rbind(analytic = g$grad, numeric = num), 4)

            [,1]   [,2]   [,3]    [,4]   [,5]    [,6]     [,7]    [,8]     [,9]
analytic 14.5077 6.6827 8.2526 -2.1064 21.196 -9.0848 -15.7974 -6.0929 -13.5267
numeric  14.5077 6.6827 8.2526 -2.1064 21.196 -9.0848 -15.7974 -6.0929 -13.5267

The draws live in an \(N \times R \times K_r\) array: Z[n, r, jj] is person \(n\)’s \(r\)-th standard-normal draw for the jj-th random coefficient — jj runs over the random subset rc, not over all \(K\) columns of \(X\); when the math writes \(z_{nk}^{(r)}\), \(k\) ranges over that subset only.

Agreement to four decimals. (The draws: Halton with a different prime base per dimension — the multivariate practice promised in Chapter 12 — person \(n\) taking segment \(n\) of each sequence, randomized by a common shift.)

13.6 Estimation

Starting values from the model hierarchy, per Chapter 6’s doctrine: the plain MNL supplies the means (it estimates a taste-blurred but well-centered version of \(\bar{\boldsymbol{\beta}}\)), and moderate spreads start the \(\sigma\)’s:

loglik_mnl <- function(beta, y, X, task_id) {
    V <- as.vector(X %*% beta); Vmax <- ave(V, task_id, FUN = max)
    eV <- exp(V - Vmax)
    sum((V - Vmax - log(rowsum(eV, task_id)[task_id]))[y == 1])
}
mnl <- optim(rep(0, K), loglik_mnl, y = y, X = X, task_id = task_id,
             method = "BFGS", control = list(fnscale = -1))

start <- c(mnl$par, rep(log(0.5), 3))
fit50 <- optim(start,
               fn = function(th, ...) msl_mixed(th, ...),
               gr = function(th, ...) msl_mixed(th, ..., want_grad = TRUE)$grad,
               y = y, X = X, task_id = task_id, id = id, n_id = n_id,
               rc = rc, Z = Z,
               method = "BFGS", control = list(fnscale = -1, maxit = 300))

The Chapter 12 doubling diagnostic — same fit at \(R = 100\), warm-started from the \(R = 50\) solution:

set.seed(131)
R2 <- 100
Z2 <- array(0, c(n_id, R2, length(rc)))
for (jj in seq_along(rc)) {
    h <- halton(n_id*R2 + 50, bases[jj])[-(1:50)]
    Z2[, , jj] <- qnorm(matrix((h + runif(1)) %% 1, n_id, R2, byrow = TRUE))
}
fit100 <- optim(fit50$par,
                fn = function(th, ...) msl_mixed(th, ...),
                gr = function(th, ...) msl_mixed(th, ..., want_grad = TRUE)$grad,
                y = y, X = X, task_id = task_id, id = id, n_id = n_id,
                rc = rc, Z = Z2,
                method = "BFGS", control = list(fnscale = -1, maxit = 300),
                hessian = TRUE)

comp <- rbind(truth  = c(bar_true, log(sd_true[rc])),
              R_50   = fit50$par,
              R_100  = fit100$par)
colnames(comp) <- c(colnames(X), paste0("log_sd_", colnames(X)[rc]))
round(comp, 3)

      brandDell brandApple ram16 ram32 screen15  price log_sd_brandApple
truth      0.50      1.000 0.600 0.900    0.300 -1.200            -0.223
R_50       0.55      1.037 0.614 0.892    0.336 -1.203            -0.402
R_100      0.55      1.034 0.615 0.891    0.337 -1.217            -0.388
      log_sd_ram32 log_sd_price
truth       -0.511       -0.357
R_50        -0.660       -0.620
R_100       -0.611       -0.656

Stable across the doubling — \(R = 100\) with Halton draws is enough here. The results table, with standard errors from the Hessian and the standard deviations mapped back from the log scale:

se <- sqrt(diag(solve(-fit100$hessian)))
tab <- data.frame(
    truth    = c(bar_true, sd_true[rc]),
    estimate = c(fit100$par[1:6], exp(fit100$par[7:9])),
    se       = c(se[1:6], exp(fit100$par[7:9]) * se[7:9])   # delta method
)
rownames(tab) <- c(colnames(X), paste0("sd_", colnames(X)[rc]))
round(tab, 3)

              truth estimate    se
brandDell       0.5    0.550 0.056
brandApple      1.0    1.034 0.066
ram16           0.6    0.615 0.056
ram32           0.9    0.891 0.062
screen15        0.3    0.337 0.044
price          -1.2   -1.217 0.060
sd_brandApple   0.8    0.679 0.080
sd_ram32        0.6    0.543 0.085
sd_price        0.7    0.519 0.099

Read the table with Section 8.8 calibration in mind. The six means recover cleanly. The three standard deviations land within roughly one-and-a-half standard errors of the truth — respectable, and note their standard errors are several times larger than the means’: spreads are simply harder to learn than centers, because they live in second moments of behavior that eight tasks per person illuminate dimly. (The delta method converts log-scale standard errors to the \(\sigma\) scale — Chapter 21 treats it properly.)

13.7 Validation Against mlogit

mlogit estimates mixed logits through its rpar argument (naming which coefficients are random and their distribution), with Halton draws and panel weighting — the same estimator, independently implemented:

library(dfidx)
library(mlogit)

md <- d
md$chid <- task_id; md$choice <- as.logical(md$choice)
md <- dfidx(md, idx = list(c("chid", "id"), "alt"))

fit_ml <- mlogit(choice ~ brand + ram + screen + price | 0, data = md,
                 rpar = c(brandApple = "n", ram32 = "n", price = "n"),
                 R = 100, halton = NA, panel = TRUE)

round(rbind(ours   = c(fit100$par[1:6], exp(fit100$par[7:9])),
            mlogit = coef(fit_ml)[c(1:6, 7:9)]), 3)

       brandDell brandApple ram16 ram32 screen15  price sd.brandApple sd.ram32
ours       0.550      1.034 0.615 0.891    0.337 -1.217         0.679    0.543
mlogit     0.551      1.036 0.614 0.889    0.337 -1.214         0.681    0.537
       sd.price
ours      0.519
mlogit    0.526

Agreement to the second decimal across all nine parameters — not exact to five decimals as in Chapter 8, and the reason is instructive rather than troubling: two MSL implementations share a model but not a simulated likelihood. Different Halton details (sequence handling, burn-in, person assignment) mean slightly different objective surfaces whose maximizers differ by simulation noise. Which raises exactly the right question: how big is that noise?

13.8 The Simulation-Error Audit

Chapter 12 promised an accounting of the error component that Hessian standard errors ignore: variation in \(\hat{\boldsymbol{\theta}}\) caused by the draws, not the data. Measuring it is refreshingly direct — re-estimate on the same data with several independent draw sets (different randomizing shifts) and watch the estimates move:

set.seed(132)
audit <- sapply(1:5, function(a) {
    Za <- array(0, c(n_id, R2, length(rc)))
    for (jj in seq_along(rc)) {
        h <- halton(n_id*R2 + 50, bases[jj])[-(1:50)]
        Za[, , jj] <- qnorm(matrix((h + runif(1)) %% 1, n_id, R2, byrow = TRUE))
    }
    optim(fit100$par,
          fn = function(th, ...) msl_mixed(th, ...),
          gr = function(th, ...) msl_mixed(th, ..., want_grad = TRUE)$grad,
          y = y, X = X, task_id = task_id, id = id, n_id = n_id,
          rc = rc, Z = Za,
          method = "BFGS", control = list(fnscale = -1, maxit = 300))$par
})
round(rbind(sim_error_sd  = apply(audit, 1, sd),
            statistical_se = se), 4)

                 [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]
sim_error_sd   0.0012 0.0034 0.0009 0.0025 0.0010 0.0020 0.0069 0.0219 0.0197
statistical_se 0.0560 0.0660 0.0563 0.0618 0.0444 0.0599 0.1174 0.1573 0.1910

The verdict at \(R = 100\) Halton draws: simulation-induced spread runs an order of magnitude (or more) below the statistical standard errors, for every parameter. The draws are not the binding source of uncertainty — which is precisely what “enough draws” means, now measured rather than assumed. Keep this audit in your repertoire; it costs a few refits and answers the referee’s question before it is asked.

13.9 Feeling the Strain

Now the honest tally this chapter owes Chapter 14. Everything above worked — and consider what it took and how it scales.

Time. Each \(R = 100\) fit above runs in minutes — entirely tolerable — but the multiplication table is merciless: the multiple-starts discipline (Chapter 11; mixture-like surfaces lurk here too²), the draw-doubling diagnostic, the simulation audit, and any replication study each multiply that base cost. Our modest audit already consumed five fits.

Dimension. Three random coefficients, independent, cost 9 parameters. The full menu — all six random, freely correlated — costs \(6 + 21 = 27\): a Cholesky factor’s worth of curvature directions for BFGS to learn, each additional random dimension another prime base of Halton draws, another axis the fixed-\(R\) simulator must cover. Published experience and our own arithmetic agree: MSL’s practical frontier sits near “a handful of random coefficients, correlations chosen sparingly” , Chapter 17.

And the deliverable gap. After all of it, the output is \((\hat{\bar{\boldsymbol{\beta}}}, \hat\Sigma)\) — the population. The individual \(\boldsymbol{\beta}_n\)’s that applications consume are not in the output; extracting them takes a second, bolt-on computation whose uncertainty accounting is awkward at best. That extraction — what it can do, and exactly where it creaks — is Chapter 14, the short chapter that closes the likelihood track and opens the Bayesian one.

One housekeeping step first — save everything the next chapters need (data with its answer key, estimates, draws):

saveRDS(list(data = d, betas_true = attr(d, "betas"),
             bar_true = bar_true, sd_true = sd_true, rc = rc,
             theta_hat = fit100$par, Z = Z2),
        file = "../data/mixed_study.rds")

13.10 Key Learnings

• The mixed MNL gives each person a coefficient vector drawn from \(N(\bar{\boldsymbol{\beta}}, \Sigma)\) (Equation 13.1): taste heterogeneity and data-driven, IIA-breaking substitution in one model — with the McFadden–Train theorem certifying the family can mimic any RUM.
• Specification is three real choices: which coefficients vary (identification from \(T\) tasks is finite), what mixing distribution (normal’s wrong-sign mass; lognormal fix deferred to Chapter 17), and whether correlations enter (Cholesky-parameterized when they do; diagonal here, with the full matrix deliberately deferred).
• The likelihood is the Chapter 11 hierarchy with an integral (Equation 13.3), estimated by Chapter 12’s MSL: fixed multivariate Halton draws (one prime base per dimension), log-mean-exp across draws, reparameterized scales.
• The MSL gradient (Equation 13.4) is a conditionally-weighted average of per-draw MNL scores — \(\mathbf{X}'(\mathbf{y}-\mathbf{p})\) yet again, with Bayes-rule weights foreshadowing Chapter 14. Verified numerically, it cut fit times several-fold.
• Recovery succeeded (spreads harder than means, visibly), mlogit agreed to the second decimal — with the residual gap correctly attributed to simulation noise, which the five-refit audit then measured: an order of magnitude below statistical error at \(R = 100\) Halton.
• The strain is real and structural: costs multiply across starts, draws, audits, and dimensions (\(27\) parameters for the full covariance), and the output is still only the population. The person-level deliverable — and likelihood’s last-mile problem with it — is next.

McFadden, Daniel, and Kenneth Train. 2000. “Mixed MNL Models for Discrete Response.” Journal of Applied Econometrics 15 (5): 447–70.

Train, Kenneth. 2016. “Mixed Logit with a Flexible Mixing Distribution.” Journal of Choice Modelling 19: 40–53.

Train, Kenneth E. 2009. Discrete Choice Methods with Simulation. 2nd ed. Cambridge University Press. https://eml.berkeley.edu/books/choice2.html.

---

1. This is a genuine modeling choice, not a computational dodge — though it is also computationally kind, and packages’ “make everything random” default deserves more skepticism than it gets: each additional random coefficient must be identified from within-person choice consistency, and \(T = 8\) tasks support only so much.

2. Non-concavity is milder here than in latent class models, but real: sign flips in \(\sigma\) (resolved by our log-parameterization) and near-flat directions when a \(\sigma\) approaches zero. Warm starts from the MNL, as above, are the standard hygiene.