\((s_1,\dots,s_C)\sim
D_C(\delta)\), where \(D_C(\delta)\) denotes the \(C\)-dimensional Dirichlet distribution with
concentration parameter vector \(\delta =
(\delta_1,\dots,\delta_C)\),
\(\alpha\sim
\text{MVN}_{P_f}(\psi,\Psi)\), where \(\text{MVN}_{P_f}\) denotes the \(P_f\)-dimensional normal distribution with
mean \(\psi\) and covariance \(\Psi\),
\(b_c \sim
\text{MVN}_{P_r}(\xi,\Xi)\), independent for all \(c\),
\(\Omega_c \sim
W^{-1}_{P_r}(\nu,\Theta)\), independent for all \(c\), where \(W^{-1}_{P_r}(\nu,\Theta)\) denotes the
\(P_r\)-dimensional inverse Wishart
distribution with \(\nu\) degrees of
freedom and scale matrix \(\Theta\),
and \(\Sigma \sim
W^{-1}_{J-1}(\kappa,\Lambda)\).
These prior distributions imply the following conditional posterior
distributions:
The class weights are drawn from the Dirichlet distribution \[\begin{equation}
(s_1,\dots,s_C)\mid \delta,z \sim D_C(\delta_1+m_1,\dots,\delta_C+m_C),
\end{equation}\] where for \(c=1,\dots,C\), \(m_c=\#\{n:z_n=c\}\) denotes the current
absolute class size.
Independently for all \(n\), we
update the allocation variables \((z_n)_n\) from their conditional
distribution \[\begin{equation}
\text{Prob}(z_n=c\mid s,\beta,b,\Omega )=\frac{s_c\phi_{P_r}(\beta_n\mid
b_c,\Omega_c)}{\sum_c s_c\phi_{P_r}(\beta_n\mid b_c,\Omega_c)}.
\end{equation}\]
The class means \((b_c)_c\) are
updated independently for all \(c\) via
\[\begin{equation}
b_c\mid \Xi,\Omega,\xi,z,\beta \sim\text{MVN}_{P_r}\left( \mu_{b_c},
\Sigma_{b_c} \right),
\end{equation}\] where \(\mu_{b_c}=(\Xi^{-1}+m_c\Omega_c^{-1})^{-1}(\Xi^{-1}\xi
+m_c\Omega_c^{-1}\bar{b}_c)\), \(\Sigma_{b_c}=(\Xi^{-1}+m_c\Omega_c^{-1})^{-1}\),
\(\bar{b}_c=m_c^{-1}\sum_{n:z_n=c}
\beta_n\).
The class covariance matrices \((\Omega_c)_c\) are updated independently
for all \(c\) via \[\begin{equation}
\Omega_c \mid \nu,\Theta,z,\beta,b \sim
W^{-1}_{P_r}(\mu_{\Omega_c},\Sigma_{\Omega_c}),
\end{equation}\] where \(\mu_{\Omega_c}=\nu+m_c\) and \(\Sigma_{\Omega_c}=\Theta^{-1} + \sum_{n:z_n=c}
(\beta_n-b_c)(\beta_n-b_c)'\).
Independently for all \(n\) and
\(t\) and conditionally on the other
components, the utility vectors \((U_{nt:})\) follow a \(J-1\)-dimensional truncated multivariate
normal distribution, where the truncation points are determined by the
choices \(y_{nt}\). To sample from a
truncated multivariate normal distribution, we apply a sub-Gibbs
sampler, following the approach of Geweke (1998): \[\begin{equation}
U_{ntj} \mid U_{nt(-j)},y_{nt},\Sigma,W,\alpha,X,\beta
\sim \mathcal{N}(\mu_{U_{ntj}},\Sigma_{U_{ntj}}) \cdot \begin{cases}
1(U_{ntj}>\max(U_{nt(-j)},0) ) & \text{if}~ y_{nt}=j\\
1(U_{ntj}<\max(U_{nt(-j)},0) ) & \text{if}~ y_{nt}\neq j
\end{cases},
\end{equation}\] where \(U_{nt(-j)}\) denotes the vector \((U_{nt:})\) without the element \(U_{ntj}\), \(\mathcal{N}\) denotes the univariate normal
distribution, \(\Sigma_{U_{ntj}} =
1/(\Sigma^{-1})_{jj}\) and \[\begin{equation}
\mu_{U_{ntj}} = W_{ntj}'\alpha + X_{ntj}'\beta_n -
\Sigma_{U_{ntj}} (\Sigma^{-1})_{j(-j)} (U_{nt(-j)} -
W_{nt(-j)}'\alpha - X_{nt(-j)}' \beta_n ),
\end{equation}\] where \((\Sigma^{-1})_{jj}\) denotes the \((j,j)\)th element of \(\Sigma^{-1}\), \((\Sigma^{-1})_{j(-j)}\) the \(j\)th row without the \(j\)th entry, \(W_{nt(-j)}\) and \(X_{nt(-j)}\) the coefficient matrices \(W_{nt}\) and \(X_{nt}\), respectively, without the \(j\)th column.
Updating the fixed coefficient vector \(\alpha\) is achieved by applying the
formula for Bayesian linear regression of the regressors \(W_{nt}\) on the regressands \((U_{nt:})-X_{nt}'\beta_n\), i.e. \[\begin{equation}
\alpha \mid \Psi,\psi,W,\Sigma,U,X,\beta \sim
\text{MVN}_{P_f}(\mu_\alpha,\Sigma_\alpha),
\end{equation}\] where \(\mu_\alpha =
\Sigma_\alpha (\Psi^{-1}\psi + \sum_{n=1,t=1}^{N,T} W_{nt} \Sigma^{-1}
((U_{nt:})-X_{nt}'\beta_n) )\) and \(\Sigma_\alpha = (\Psi^{-1} + \sum_{n=1,t=1}^{N,T}
W_{nt}\Sigma^{-1} W_{nt}^{'} )^{-1}\).
Analogously to \(\alpha\), the
random coefficients \((\beta_n)_n\) are
updated independently via \[\begin{equation}
\beta_n \mid \Omega,b,X,\Sigma,U,W,\alpha \sim
\text{MVN}_{P_r}(\mu_{\beta_n},\Sigma_{\beta_n}),
\end{equation}\] where \(\mu_{\beta_n}
= \Sigma_{\beta_n} (\Omega_{z_n}^{-1}b_{z_n} + \sum_{t=1}^{T} X_{nt}
\Sigma^{-1} (U_{nt:}-W_{nt}'\alpha) )\) and \(\Sigma_{\beta_n} = (\Omega_{z_n}^{-1} +
\sum_{t=1}^{T} X_{nt}\Sigma^{-1} X_{nt}^{'} )^{-1}\)
.
The error term covariance matrix \(\Sigma\) is updated by means of \[\begin{equation}
\Sigma \mid \kappa,\Lambda,U,W,\alpha,X,\beta \sim
W^{-1}_{J-1}(\kappa+NT,\Lambda+S), \\
\end{equation}\] where \(S =
\sum_{n=1,t=1}^{N,T} \varepsilon_{nt} \varepsilon_{nt}'\) and
\(\varepsilon_{nt} = (U_{nt:}) -
W_{nt}'\alpha - X_{nt}'\beta_n\).
Parameter normalization
Samples obtained from the updating scheme described above lack
identification (except for \(s\) and
\(z\) draws), compare to the vignette
on the model definition. Therefore, subsequent to the sampling, the
following normalizations are required for the \(i\)th updates in each iterations \(i\):
\(\alpha^{(i)} \cdot
\omega^{(i)}\),
\(b_c^{(i)} \cdot
\omega^{(i)}\), \(c=1,\dots,C\),
\(U_{nt}^{(i)} \cdot
\omega^{(i)}\), \(n =
1,\dots,N\), \(t =
1,\dots,T\),
\(\beta_n^{(i)} \cdot
\omega^{(i)}\), \(n =
1,\dots,N\),
\(\Omega_c^{(i)} \cdot
(\omega^{(i)})^2\), \(c=1,\dots,C\), and
\(\Sigma^{(i)} \cdot
(\omega^{(i)})^2\),
where either \(\omega^{(i)} =
\sqrt{\text{const} / (\Sigma^{(i)})_{jj}}\) with \((\Sigma^{(i)})_{jj}\) the \(j\)th diagonal element of \(\Sigma^{(i)}\), \(1\leq j \leq J-1\), or alternatively \(\omega^{(i)} = \text{const} /
\alpha^{(i)}_p\) for some coordinate \(1\leq p \leq P_f\) of the \(i\)th draw for the coefficient vector \(\alpha\). Here, \(\text{const}\) is any positive constant
(typically 1). The preferences will be flipped if \(\omega^{(i)} < 0\), which only is the
case if \(\alpha^{(i)}_p < 0\).
Burn-in and thinning
The theory behind Gibbs sampling constitutes that the sequence of
samples produced by the updating scheme is a Markov chain with
stationary distribution equal to the desired joint posterior
distribution. It takes a certain number of iterations for that
stationary distribution to be approximated reasonably well. Therefore,
it is common practice to discard the first \(B\) out of \(R\) samples (the so-called burn-in period).
Furthermore, correlation between nearby samples should be expected. In
order to obtain independent samples, we consider only every \(Q\)th sample when computing Gibbs sample
statistics like expectation and standard deviation. The independence of
the samples can be verified by computing the serial correlation and the
convergence of the Gibbs sampler can be checked by considering trace
plots, see below.
