mnp_study <- readRDS("../data/mnp_study.rds")
mnp <- mnp_study$data
X <- model.matrix(~ alt + ram + screen + price, data = mnp)[, -1]
K <- ncol(X)
## rows of each alternative (the design is balanced and ordered: @sec-data contract)
i1 <- which(mnp$alt == "Acer"); i2 <- which(mnp$alt == "Dell")
i3 <- which(mnp$alt == "Apple")
Xd2 <- X[i2, ] - X[i1, ] # attributes of (Dell - Acer), one row per task
Xd3 <- X[i3, ] - X[i1, ] # attributes of (Apple - Acer)
chosen <- as.integer(mnp$alt[mnp$choice == 1]) # 1=Acer, 2=Dell, 3=Apple
n_task <- length(chosen)19 Estimating the MNP
One model, both toolkits. The MNP’s likelihood is a wall of normal orthant integrals (Equation 18.2); this chapter scales it twice — once with Part II’s machinery (simulate the integral, maximize the simulated likelihood) and once with Part III’s (augment the integral away, Gibbs the rest). Beyond estimating the MNP, the chapter is the book’s thesis made concrete: the model is a set of assumptions; the estimator is a separate choice; and a well-built toolkit lets you make that choice on the merits. By the end, the two routes will sit in one table, agreeing with each other to the second decimal.
Load the study saved in Chapter 18, and build the differenced objects that both routes consume — everything here runs on the identification analysis of Chapter 18, so the differencing is not preprocessing, it is the identified model:
And compute the identified truth this chapter must recover — the Chapter 18 arithmetic executed: difference \(\Omega^\ast\) against the base, normalize the covariance by its first diagonal, and rescale \(\boldsymbol{\beta}^\ast\) by the square root of that diagonal — \(\boldsymbol{\beta}^\ast_{\text{id}} = \boldsymbol{\beta}^\ast/\sqrt{(M\Omega^\ast M')_{11}}\) and \(\tilde\Omega^\ast = M\Omega^\ast M'/(M\Omega^\ast M')_{11}\); coefficients live on the utility scale, covariances on the utility-squared scale, so they pay the scale tax at different rates:
M <- rbind(c(-1, 1, 0), c(-1, 0, 1))
Ot <- M %*% mnp_study$Omega_true %*% t(M) # differenced covariance
sc <- sqrt(Ot[1, 1]) # the scale normalization
beta_star <- mnp_study$beta_true / sc
omega_star <- c(w12 = Ot[1, 2], w22 = Ot[2, 2]) / Ot[1, 1]
round(c(beta_star, omega_star), 3) dell apple ram16 ram32 screen15 price w12 w22
0.559 1.118 0.671 1.006 0.335 -1.342 0.500 2.500 Those eight numbers — six scaled coefficients, one differenced covariance, one differenced variance ratio — are everything the data can reveal, and the scoreboard for both routes.
19.1 Route 1: Simulated Maximum Likelihood via GHK
19.1.1 Why the obvious simulator fails
Chapter 1’s frequency simulator applies directly: draw errors, count how often the chosen alternative wins. It is unbiased — and useless inside an optimizer, for two reasons that Chapter 12 equips us to name instantly. It is a step function of the parameters (nudge \(\boldsymbol{\beta}\) slightly and either no draw changes its vote or one flips discretely), so gradients are zero almost everywhere and undefined elsewhere — BFGS starves. And for tasks where the chosen alternative was unlikely, all \(R\) draws can vote no, returning \(\widehat{P} = 0\) and \(\log \widehat{P} = -\infty\) — Chapter 8’s underflow catastrophe, by construction rather than accident. Smoothing the indicator (the logit-kernel AR simulator) patches both at the cost of a tuning parameter and bias (Train 2009, sec. 5.6). The field’s solution is better than a patch.
19.1.2 The GHK simulator
The GHK simulator (Geweke–Hajivassiliou–Keane) computes the orthant probability by recursive conditioning, and it is Chapter 2’s truncated-normal trick wearing its work clothes. The probability that alternative \(j\) was chosen is the probability that a bivariate normal vector \(\mathbf{w} = (w_1, w_2)\) — the utility differences against \(j\) — lands entirely below zero. Factor the joint event through the Cholesky root of its covariance, \(A_j \tilde\Omega A_j' = LL'\), writing \(\mathbf{w} = \boldsymbol{\mu} + L\boldsymbol{\eta}\) with \(\boldsymbol{\eta}\) standard normal. Because \(L\) is lower triangular with positive diagonal, \(w_1 < 0 \iff \eta_1 < b_1 := -\mu_1/\ell_{11}\); and given \(\eta_1\), \(w_2 < 0 \iff \eta_2 < b_2(\eta_1) := (-\mu_2 - \ell_{21}\eta_1)/\ell_{22}\). The law of total probability, plus the independence of \(\eta_2\) from \(\eta_1\), then gives:
\[ \Pr(w_1 < 0, w_2 < 0) = \underbrace{\Pr\left( \eta_1 < b_1 \right)}_{p_1} \times \mathrm{E}\left[ \underbrace{\Pr\left( \eta_2 < b_2(\eta_1) \right)}_{p_2(\eta_1)} \,\middle|\, \eta_1 < b_1 \right] \tag{19.1}\]
where \(b_1 = -\mu_1 / \ell_{11}\) and \(b_2(\eta_1) = (-\mu_2 - \ell_{21}\eta_1)/\ell_{22}\): the first coordinate’s constraint is a univariate normal probability, exactly; the second’s is a univariate probability given the first, averaged over truncated draws of \(\eta_1\). GHK draws \(\eta_1\) from its truncated distribution by the inverse-CDF trick of Equation 2.2 — qnorm(u * p1) maps a uniform into exactly the admissible region — and averages \(p_1 \cdot p_2(\eta_1^{(r)})\) over \(R\) draws. Every ingredient is a pnorm/qnorm call: smooth in the parameters (BFGS eats), strictly positive (each factor is a genuine probability), and unbiased. For \(J\) alternatives the recursion just continues (Train 2009, sec. 5.6.3 writes out the general case) — each coordinate truncated given the previous ones, the triangularity of \(L\) making each step univariate; you can see in that description why Section 2.4 called Cholesky triangularity “the same idea that powers GHK.”
One bookkeeping matter: the differences-against-\(j\) vector for each chosen alternative is a fixed linear map \(A_j\) of the base differences \((d_2, d_3)\) — e.g., for \(j=2\): \((U_1 - U_2, U_3 - U_2) = (-d_2, \; d_3 - d_2)\). All three maps, for the record: \(A_1 = I_2\) (Acer chosen: \((d_2, d_3) < \mathbf{0}\) already); \(A_2 = \begin{pmatrix} -1 & 0 \\ -1 & 1 \end{pmatrix}\), so \(A_2 \mathbf{d} = (-d_2, \; d_3 - d_2)\); \(A_3 = \begin{pmatrix} 0 & -1 \\ 1 & -1 \end{pmatrix}\). Row \(i\) of \(A_j\) simply rewrites \(U_k - U_j\) in terms of the base differences. So means transform per task, but the covariance \(A_j \tilde\Omega A_j'\) (hence its Cholesky) is computed once per chosen-alternative type per likelihood call — three small factorizations, not four thousand.
## GHK simulated log-likelihood for the 3-alternative MNP
## theta = (beta [K], a, log b) with differenced-normalized chol Lt = [1 0; a b]
ghk_loglik <- function(theta, Xd2, Xd3, chosen, U1) {
K <- ncol(Xd2)
beta <- theta[1:K]; a <- theta[K+1]; b <- exp(theta[K+2])
Ot_id <- rbind(c(1, a), c(a, a^2 + b^2)) # identified covariance: first diagonal fixed at 1 (this is Omega-tilde)
m2 <- as.vector(Xd2 %*% beta) # mean of d2 = U_dell - U_acer
m3 <- as.vector(Xd3 %*% beta)
A <- list(diag(2), # chosen = Acer: (d2, d3) < 0
rbind(c(-1, 0), c(-1, 1)), # chosen = Dell: (-d2, d3-d2) < 0
rbind(c(0, -1), c(1, -1))) # chosen = Apple: (-d3, d2-d3) < 0
ll <- 0
for (j in 1:3) {
idx <- which(chosen == j)
mu <- cbind(m2[idx], m3[idx]) %*% t(A[[j]])
Lj <- t(chol(A[[j]] %*% Ot_id %*% t(A[[j]])))
u <- U1[idx, , drop = FALSE] # fixed uniforms (@sec-msl!)
b1 <- (-mu[, 1]) / Lj[1, 1]
p1 <- pmax(pnorm(b1), 1e-12)
e1 <- qnorm(pmax(u * p1, 1e-14)) # truncated draws via the inverse-CDF trick (eq-trunc)
b2 <- (-mu[, 2] - Lj[2, 1] * e1) / Lj[2, 2]
Phat <- pmax(rowMeans(p1 * pnorm(b2)), 1e-300)
ll <- ll + sum(log(Phat))
}
ll
}Note the Chapter 12 discipline embedded: the uniforms U1 are generated once and passed in — fixed draws, rule one — and the covariance enters through its Cholesky elements with the diagonal log-parameterized, so every candidate the optimizer visits implies a valid matrix. Before trusting it, the mandatory audit — GHK against an independent oracle, R’s numerical bivariate-normal routine, on a handful of tasks:
library(mvtnorm)
set.seed(190)
U1 <- matrix(runif(n_task * 100), n_task, 100)
check <- sapply(1:5, function(i) {
mu <- c(sum(Xd2[i, ] * beta_star), sum(Xd3[i, ] * beta_star))
Ot_id <- rbind(c(1, omega_star["w12"]),
c(omega_star["w12"], omega_star["w22"]))
A <- list(diag(2), rbind(c(-1,0), c(-1,1)), rbind(c(0,-1), c(1,-1)))[[chosen[i]]]
exact <- pmvnorm(upper = c(0, 0), mean = as.vector(A %*% mu),
sigma = A %*% Ot_id %*% t(A))
ghk <- { # one-task GHK, 100 draws
Lj <- t(chol(A %*% Ot_id %*% t(A))); mm <- as.vector(A %*% mu)
p1 <- pnorm(-mm[1]/Lj[1,1])
e1 <- qnorm(U1[i, ] * p1)
mean(p1 * pnorm((-mm[2] - Lj[2,1]*e1)/Lj[2,2]))
}
c(exact = exact, ghk = ghk)
})
round(check, 4) [,1] [,2] [,3] [,4] [,5]
exact 0.6359 0.8164 0.6167 0.3432 0.5639
ghk 0.6170 0.8115 0.5975 0.3296 0.5641Agreement in the third decimal with 100 draws — certified. Now the fit, with the covariance parameters box-bounded (the Chapter 10 device) to keep the search in well-conditioned territory:
fit_ghk <- optim(c(rep(0, K), 0, 0), ghk_loglik,
Xd2 = Xd2, Xd3 = Xd3, chosen = chosen, U1 = U1,
method = "L-BFGS-B",
lower = c(rep(-Inf, K), -5, -2.5),
upper = c(rep( Inf, K), 5, 2.5),
control = list(fnscale = -1, maxit = 500), hessian = TRUE)
se_ghk <- sqrt(diag(solve(-fit_ghk$hessian)))
est_ghk <- c(fit_ghk$par[1:K],
w12 = fit_ghk$par[K+1],
w22 = fit_ghk$par[K+1]^2 + exp(fit_ghk$par[K+2])^2)
round(data.frame(truth = c(beta_star, omega_star),
ghk = est_ghk,
se = c(se_ghk[1:K], se_ghk[K+1], NA)), 3) truth ghk se
dell 0.559 0.502 0.065
apple 1.118 1.142 0.066
ram16 0.671 0.805 0.057
ram32 1.006 1.173 0.067
screen15 0.335 0.350 0.039
price -1.342 -1.423 0.072
w12 0.500 0.418 0.121
w22 2.500 2.610 NARecovery of all eight identified quantities, including the covariance parameters that carry the substitution story — the estimated \(\tilde\omega_{12}\) is the data confirming that Acer and Dell share unobserved appeal. (The w22 transformation’s standard error needs the delta method; Chapter 21, or route 2, where it’s free.)
19.2 Route 2: Gibbs with Data Augmentation
Write the differenced system explicitly first, because it is the object the whole route manipulates:
\[\mathbf{d}_{nt} = \tilde{X}_{nt}\boldsymbol{\beta}+ \mathbf{e}_{nt}, \qquad \mathbf{e}_{nt} \sim N(\mathbf{0}, W)\]
where \(\mathbf{d}_{nt} = (d_2, d_3)\) stacks \(d_2 = U_{\text{Dell}} - U_{\text{Acer}}\) and \(d_3 = U_{\text{Apple}} - U_{\text{Acer}}\), the rows of \(\tilde{X}_{nt}\) are the corresponding rows of Xd2 and Xd3, and \(W\) is Chapter 18’s differenced covariance before the \(\tilde\omega_{11} = 1\) convention — the unnormalized object the sampler is free to draw. The identified quantities are \(\tilde\Omega = W/w_{11}\) and \(\boldsymbol{\beta}/\sqrt{w_{11}}\); the sampler never imposes this, and the post-processing recovers it.
Now forget optimization entirely. The Chapter 16 Albert–Chib logic, scaled up per McCulloch and Rossi (1994): conditional on the latent differenced utilities \(\mathbf{d}_{nt} = (d_2, d_3)_{nt}\), the MNP is a two-equation normal regression system — conjugate everywhere; conditional on parameters and the observed choice, each task’s \(\mathbf{d}_{nt}\) is bivariate normal truncated to the region its choice implies. Three blocks:
- Latent utilities. The regions, read off the choice rule: Acer chosen \(\Rightarrow d_2 < 0, d_3 < 0\); Dell \(\Rightarrow d_2 > \max(0, d_3)\); Apple \(\Rightarrow d_3 > \max(0, d_2)\). Sampling a truncated bivariate normal is done by one more level of Gibbs: draw \(d_2\) given \(d_3\) (a univariate truncated normal — the conditional mean and variance are the textbook bivariate-normal formulas, the truncation from the region; concretely, if \((d_2, d_3) \sim N\big((m_2, m_3), W\big)\) then \(d_2 \mid d_3 \sim N\!\big(m_2 + \tfrac{w_{12}}{w_{22}}(d_3 - m_3), \; w_{11} - \tfrac{w_{12}^2}{w_{22}}\big)\), and symmetrically for \(d_3 \mid d_2\) — the code’s
c21/v21andc12/v12), then \(d_3\) given \(d_2\). All 4,000 tasks in two vectorizedrtnorm()calls. - \(\boldsymbol{\beta}\). Given the \(\mathbf{d}\)’s: generalized least squares dressed as Bayes — whiten the two-equation system by \(W^{-1/2}\) (Cholesky again) and apply the conjugate regression update of Chapter 16.
- \(W\). Given residuals: a \(2 \times 2\) inverse-Wishart.
And the McCulloch–Rossi move that makes it elegant: sample the unidentified model, report the identified one. No constraint forces \(w_{11} = 1\) during sampling — the chain wanders the scale ridge freely, conjugacy intact — and every saved draw is normalized afterward: \(\boldsymbol{\beta}/ \sqrt{w_{11}}\), \(W / w_{11}\). Identified functions of unidentified draws are perfectly valid posterior draws of the identified quantities.1
rtnorm <- function(n, mean, sd, a, b) { # @sec-draws
Fa <- pnorm(a, mean, sd); Fb <- pnorm(b, mean, sd)
u <- pmin(pmax(runif(n, Fa, Fb), 1e-12), 1 - 1e-12)
qnorm(u, mean, sd)
}
rinvwishart <- function(nu, S) { # @sec-hb-mnl
K <- nrow(S); L <- t(chol(solve(S)))
A <- matrix(0, K, K)
diag(A) <- sqrt(rchisq(K, nu - (1:K) + 1))
A[lower.tri(A)] <- rnorm(K*(K-1)/2)
solve(L %*% A %*% t(A) %*% t(L))
}
gibbs_mnp <- function(Xd2, Xd3, chosen, S_iter, A0 = 0.01,
nu0 = 5, V0 = diag(2) * 5) {
K <- ncol(Xd2); n <- nrow(Xd2)
beta <- rep(0, K); W <- diag(2)
d2 <- ifelse(chosen == 2, 1, -0.5) # any region-respecting start
d3 <- ifelse(chosen == 3, 1, -0.5)
out <- list(beta = matrix(NA_real_, S_iter, K),
W = matrix(NA_real_, S_iter, 3))
for (s in 1:S_iter) {
## -- latent utilities: two vectorized truncated-normal sweeps --
m2 <- as.vector(Xd2 %*% beta); m3 <- as.vector(Xd3 %*% beta)
c21 <- W[1,2]/W[2,2]; v21 <- W[1,1] - W[1,2]^2/W[2,2]
c12 <- W[1,2]/W[1,1]; v12 <- W[2,2] - W[1,2]^2/W[1,1]
lo2 <- ifelse(chosen == 2, pmax(0, d3), -Inf)
hi2 <- ifelse(chosen == 2, Inf, ifelse(chosen == 1, 0, d3))
d2 <- rtnorm(n, m2 + c21*(d3 - m3), sqrt(v21), lo2, hi2)
lo3 <- ifelse(chosen == 3, pmax(0, d2), -Inf)
hi3 <- ifelse(chosen == 3, Inf, ifelse(chosen == 1, 0, d2))
d3 <- rtnorm(n, m3 + c12*(d2 - m2), sqrt(v12), lo3, hi3)
## -- beta: whiten the 2-equation system, conjugate regression --
Li <- solve(t(chol(W)))
Y1 <- Li[1,1]*d2 + Li[1,2]*d3; Y2 <- Li[2,1]*d2 + Li[2,2]*d3
Xw1 <- Li[1,1]*Xd2 + Li[1,2]*Xd3; Xw2 <- Li[2,1]*Xd2 + Li[2,2]*Xd3
Vb <- solve(crossprod(Xw1) + crossprod(Xw2) + diag(A0, K))
mb <- Vb %*% (crossprod(Xw1, Y1) + crossprod(Xw2, Y2))
beta <- as.vector(mb + t(chol(Vb)) %*% rnorm(K))
## -- W: inverse-Wishart on residuals --
e2 <- d2 - as.vector(Xd2 %*% beta); e3 <- d3 - as.vector(Xd3 %*% beta)
Sres <- rbind(c(sum(e2^2), sum(e2*e3)), c(sum(e2*e3), sum(e3^2)))
W <- rinvwishart(nu0 + n, V0 + Sres)
out$beta[s, ] <- beta
out$W[s, ] <- c(W[1,1], W[1,2], W[2,2])
}
out
}set.seed(191)
g <- gibbs_mnp(Xd2, Xd3, chosen, S_iter = 3000)
post <- 501:3000
scd <- sqrt(g$W[post, 1]) # per-draw normalization
beta_id <- g$beta[post, ] / scd
tab_gibbs <- data.frame(
truth = c(beta_star, omega_star),
gibbs = c(colMeans(beta_id),
mean(g$W[post, 2] / g$W[post, 1]),
mean(g$W[post, 3] / g$W[post, 1])),
post_sd = c(apply(beta_id, 2, sd),
sd(g$W[post, 2] / g$W[post, 1]),
sd(g$W[post, 3] / g$W[post, 1]))
)
rownames(tab_gibbs) <- c(colnames(X), "w12", "w22")
round(tab_gibbs, 3) truth gibbs post_sd
altDell 0.559 0.513 0.059
altApple 1.118 1.148 0.063
ram16 0.671 0.809 0.056
ram32 1.006 1.175 0.063
screen15 0.335 0.353 0.037
price -1.342 -1.428 0.069
w12 0.500 0.436 0.112
w22 2.500 2.683 0.367Recovery again — and notice what this sampler did not contain: no likelihood evaluation (the orthant integral never appears — Chapter 16’s augmentation promise, delivered at multinomial scale), no Metropolis step, no acceptance rate, no tuning constant of any kind, and posterior uncertainty for the covariance parameters included, no delta method required.
19.3 The Two Routes, Head to Head
head2head <- data.frame(
truth = c(beta_star, omega_star),
ghk_msl = est_ghk,
gibbs = tab_gibbs$gibbs
)
rownames(head2head) <- c(colnames(X), "w12", "w22")
round(head2head, 3) truth ghk_msl gibbs
altDell 0.559 0.502 0.513
altApple 1.118 1.142 1.148
ram16 0.671 0.805 0.809
ram32 1.006 1.173 1.175
screen15 0.335 0.350 0.353
price -1.342 -1.423 -1.428
w12 0.500 0.418 0.436
w22 2.500 2.610 2.683Agreement between the routes to roughly the second decimal on every identified quantity — two philosophies, two codebases, two entirely different computational strategies (one never draws a latent variable, the other never evaluates a probability), one answer. The residual daylight between either column and the truth is this dataset’s sampling noise, which the two estimators share because they share the data.
The choice between routes, then, rides on everything around the point estimates. Tallying this chapter’s experience: the Gibbs route was faster here, tuning-free, and delivered exact posterior uncertainty for every identified function — and, decisively for what comes next, it extends: bolting the Chapter 17 hierarchy onto it is Chapter 20, a short chapter precisely because augmentation composes. The GHK route earns its keep where the Bayesian route strains — very large \(J\) (the truncated-utilities block grows with alternatives), classical inference requirements, or embedding in moment-based estimators — and GHK itself is a permanent resident of the choice modeler’s toolkit (Train 2009, Ch. 5). But for the workhorse cases, the augmentation route is the one modern practice reaches for, and having built both, you know exactly why.
19.4 Validation Against bayesm
The reference implementation of route 2 is rmnpGibbs (Rossi, Allenby, and McCulloch 2005). Its conventions differ from ours in one consequential way — it differences against the last alternative (base \(= J\)), where we used the first — so validating requires restating the identified truth in its parameterization: differences against Apple, ASCs for Acer and Dell, normalization by the first differenced variance in that system. Getting a reference package’s conventions right is not busywork; it is most of what “validation” means in practice, and mismatched differencing conventions produce exactly the kind of sign-flipped, rescaled nonsense that looks like a broken sampler.2
library(bayesm)
na <- 4 # alt-varying attributes
Xa <- matrix(0, n_task, 3*na)
for (a_i in 1:na) { # columns: variable-major, alt within
Xa[, (a_i-1)*3 + 1] <- X[i1, 2 + a_i]
Xa[, (a_i-1)*3 + 2] <- X[i2, 2 + a_i]
Xa[, (a_i-1)*3 + 3] <- X[i3, 2 + a_i]
}
Xc <- createX(p = 3, na = na, nd = NULL, Xa = Xa, Xd = NULL,
INT = TRUE, DIFF = TRUE, base = 3)
set.seed(192)
outb <- rmnpGibbs(Data = list(p = 3, y = chosen, X = Xc),
Mcmc = list(R = 3000, keep = 1, nprint = 0))## identified truth in bayesm's (base = Apple) parameterization
M3 <- rbind(c(1, 0, -1), c(0, 1, -1))
Ot3 <- M3 %*% mnp_study$Omega_true %*% t(M3)
sc3 <- sqrt(Ot3[1, 1])
truth_b <- c(int_acer = (0 - 1.0)/sc3, # base = Apple: each ASC is its raw value minus Apple's (1.0)
int_dell = (0.5 - 1.0)/sc3,
mnp_study$beta_true[3:6]/sc3)
## reuse post = 501:3000: same burn-in convention for bayesm's chain
s11 <- outb$sigmadraw[post, 1]
round(rbind(truth = truth_b,
bayesm = colMeans(outb$betadraw[post, ] / sqrt(s11))), 3) int_acer int_dell ram16 ram32 screen15 price
truth -0.707 -0.354 0.424 0.636 0.212 -0.849
bayesm -0.711 -0.394 0.497 0.726 0.218 -0.884round(c(w12 = mean(outb$sigmadraw[post, 2] / s11),
w22 = mean(outb$sigmadraw[post, 4] / s11),
truth_w12 = Ot3[1,2]/Ot3[1,1], truth_w22 = Ot3[2,2]/Ot3[1,1]), 3) w12 w22 truth_w12 truth_w22
0.831 1.057 0.800 1.000 The package agrees — with our estimates and with the (re-parameterized) truth. Same normalize-after-sampling treatment of its draws, note: bayesm also samples unidentified and leaves the division to you.
19.5 Key Learnings
- The frequency simulator fails inside optimizers for nameable reasons — non-smooth, can return zero — and GHK (Equation 19.1) repairs both: recursive conditioning through the Cholesky root turns an orthant probability into a product of exact univariate probabilities averaged over truncated draws. Smooth, positive, unbiased; audited against
pmvnormbefore use. - GHK-MSL is Chapter 12 verbatim: fixed uniforms, log-parameterized Cholesky, box bounds, BFGS — and it recovered all eight identified quantities.
- The augmentation route dissolves the integral instead: truncated-bivariate utility draws (Gibbs-within-Gibbs, two vectorized
rtnormsweeps), a whitened conjugate regression for \(\boldsymbol{\beta}\), an inverse-Wishart for \(W\) — no tuning, no likelihood. Sample the unidentified model, normalize every draw afterward: identified posteriors from an unconstrained chain, per McCulloch and Rossi (1994). - Head to head: second-decimal agreement across routes. The differences that matter are architectural — Gibbs extends to the hierarchy (next chapter) and prices uncertainty exactly; GHK serves large \(J\) and classical pipelines.
- Validation against
rmnpGibbssucceeded only after matching its base-alternative convention — reference-implementation conventions are half the work, and mismatches masquerade as broken code.
This is the resolution of Chapter 18’s warning about samplers on unidentified spaces: the raw traces of \(\boldsymbol{\beta}\) and \(W\) indeed drift along the scale ridge and look unconverged — diagnose the normalized series instead. Try plotting both; it is a five-line experiment and a permanent inoculation.↩︎
Ask me how I know. (The first attempt at this comparison used base \(= 1\) and produced beautifully converged garbage.)↩︎