AppReferences

# Ultimate Pólya Gamma Samplers – Efficient MCMC for possibly imbalanced
binary and categorical data

Sylvia Frühwirth-Schnatter (WU Wien)
Gregor Zens (WU Wien)
Helga Wagner (JKU Linz)

###### Abstract

Modeling binary and categorical data is one of the most commonly encountered
tasks of applied statisticians and econometricians. While Bayesian methods in
this context have been available for decades now, they often require a high
level of familiarity with Bayesian statistics or suffer from issues such as
low sampling efficiency. To contribute to the accessibility of Bayesian models
for binary and categorical data, we introduce novel latent variable
representations based on Pólya Gamma random variables for a range of commonly
encountered discrete choice models. From these latent variable
representations, new Gibbs sampling algorithms for binary, binomial and
multinomial logistic regression models are derived. All models allow for a
conditionally Gaussian likelihood representation, rendering extensions to more
complex modeling frameworks such as state space models straight-forward.
However, sampling efficiency may still be an issue in these data augmentation
based estimation frameworks. To counteract this, MCMC boosting strategies are
developed and discussed in detail. The merits of our approach are illustrated
through extensive simulations and a real data application.

Keywords: Bayesian, Data augmentation, Discrete choice, Gibbs sampling, MCMC
boosting.

##  1 Introduction

Applied statisticians and econometricians commonly have to deal with modelling
binary or categorical outcome variables. Widely used tools for analyzing such
data include logit and probit models, binomial and multinomial logit (MNL)
regression-type models as well as more flexible models such as binary state
space or random effects models. Bayesian approaches toward inference are very
useful in this context, in particular for uncertainty quantification, but as
opposed to Gaussian outcome variables their implementation is demanding from a
computational viewpoint.

One strategy to implement sampling-based inference relies on importance
sampling (Zellner & Rossi, 1984) or Metropolis-Hastings (MH) algorithms (Rossi
et al., 2005), exploiting directly the non-Gaussian likelihood. However, these
algorithms require careful tuning and substantial experience with Bayesian
computation for more complex models such as state space models.

Routine Bayesian computation for these type of data more often relies on
Markov Chain Monte Carlo (MCMC) algorithms based on data augmentation. As
shown by the seminal paper of Albert & Chib (1993), the binary probit model
admits a latent variable representation where the latent variable equation is
linear in the unknown parameters with an error term following a standard
normal distribution. As the latent variables are easily sampled for known
parameters, the latent variable representation admits a straightforward Gibbs
sampler using one level of data augmentation, where the unknown parameters are
sampled from a conditionally Gaussian model. This strategy works also for more
complex models such as probit state space or random effects models.

However, MCMC estimation based on data augmentation is less straightforward
for a logit model which still admits a latent variable representation that is
linear in the unknown parameters, but exhibits an error following a logistic
distribution. Related latent variable representations with non-Gaussian errors
exist for MNL models (see e.g. Frühwirth-Schnatter & Frühwirth (2010)) and
binomial models (Fussl et al., 2013). While the latent variables can be
sampled straightforwardly for any of these latent variable representations,
sampling the unknown parameters is far from straightforward due to the non-
Gaussian errors.

A common resolution relies on a scale-mixture representation of the non-
Gaussian error distribution and introduces the scales as a second level of
data augmentation. Very conveniently, the unknown model parameters are then
sampled from a conditionally Gaussian model. Examples include a representation
of the logistic distribution involving the Kolmogoroff-Smirnov distribution
(Holmes & Held, 2006) and very accurate finite scale-mixture approximations
(Frühwirth-Schnatter & Frühwirth, 2007, 2010; Frühwirth-Schnatter et al.,
2009). A seminal paper in this context which avoids any explicit latent
variable representation and works instead with a single level of data
augmentation is Polson et al. (2013). They derive the Pólya-Gamma sampler,
using a mixture representation of the non-Gaussian likelihood of the marginal
model involving the Pólya-Gamma distribution. In the present paper, we propose
a new sampler involving the Pólya-Gamma distribution introduced by Polson et
al. (2013). However, as opposed to working with the marginal model, we
introduce a new mixture representation of the logistic distribution based on
the Pólya-Gamma distribution in the latent variable representation of the
logit model. Similarly as Holmes & Held (2006) and Frühwirth-Schnatter &
Frühwirth (2010), we expand data augmentation and introduce the Pólya-Gamma
mixing variables as a second set of latent variables. Our new Pólya-Gamma
mixture representation of the logit model has the advantage that the joint
(conditional) posterior distribution of all data-augmented variables is easy
to sample from as the Pólya-Gamma mixing variable conditional on the latent
variable follows a tilted Pólya-Gamma distribution. Conditional on all data-
augmented variables, this allows to sample the unknown model parameters from a
conditionally Gaussian model, even for complex settings such as state space or
random effects models.

One challenge when working with any MCMC method based on data augmentation is
poor mixing, in particular in the case of imbalanced data, where the success
probability is either close to zero or one for the majority of the observed
data points, see the excellent work of Johndrow et al. (2019). Neither the
original Pólya-Gamma sampler of Polson et al. (2013) with a single level of
data augmentation nor our new Pólya-Gamma sampler with two levels of data
augmentation is an exception to this rule.

In the present paper we introduce imbalanced marginal data augmentation (iMDA)
as a boosting strategy to make not only our new sampler, but also the original
probit sampler of Albert & Chib (1993) robust to possibly imbalanced data.
This strategy is inspired by earlier work on marginal data augmentation for
binary and categorical data (van Dyk & Meng, 2001; Imai & van Dyk, 2005;
McCulloch et al., 2000). Starting from a latent variable representation of the
binary model, we move to an unidentified version of the latent variable
representation with the help of two unidentified parameters. One parameter is
a global scale parameter for the latent variable which usually takes the value
one for the identified model. Introducing such an unidentified parameter has
been shown to improve mixing considerably, among others by Imai & van Dyk
(2005). However, this strategy does not resolve the situation of highly
imbalanced data. To this aim, we introduce in addition an unknown threshold
parameter which takes the value zero for the identified model. As will be
shown, this parameter restores balance in the unidentified model and in this
way improves mixing considerably for unbalanced data, both for the sampler of
Albert & Chib (1993) as well as for our new Pólya-Gamma samplers. As iMDA only
works in the context of a latent variable representation, this strategy cannot
be applied to the original Pólya-Gamma sampler of Polson et al. (2013) due to
the lack of such a representation.

The new Pólya-Gamma representation of the logit model is very generic and is
easily combined with iMDA, not only for binary regression models, but also for
more flexible models such as binary state space models, as the resulting full
conditional of the state process is Gaussian which allows multi-move sampling
of the entire state process (Frühwirth-Schnatter, 1994; Carter & Kohn, 1994).
We name any sampling strategy that combines a Pólya-Gamma mixture
representation with iMDA an ultimate Pólya-Gamma sampler due to its
efficiency.

A further contribution of the present paper is to show that such an ultimate
Pólya-Gamma sampler can be derived for more general non-Gaussian problems,
including multinomial modelling of categorical data and modelling binomial
data. For the MNL model usually a logit model based on a partial dRUM
representation is applied to sample the category specific parameters, see e.g.
Holmes & Held (2006); Frühwirth-Schnatter & Frühwirth (2010); Polson et al.
(2013). In the present paper, we derive a new latent variable representation
of the MNL model involving a partial dRUM representation with two latent
variables instead of a single one. One latent variable equation is linear in
the unknown parameters and involves a logistic error distribution. As for
binary logit models, we use the Pólya-Gamma mixture representation of the
logistic distribution and introduce the mixing variables as additional latent
variables. For binomial models, a latent variable representation has been
introduced earlier by Fussl et al. (2013), but this representation does not
provide an explicit choice equation. However, since an explicit choice
equation is needed to apply iMDA, we derive a new latent variable
representation for binomial data involving two latent variables. Both latent
variable equations are linear in the unknown parameters and involve a
generalized logistic error distribution. Extending the strategy introduced for
binary logit models, we use Pólya-Gamma mixture representations of the
generalized logistic distribution and introduce the Pólya-Gamma mixing
variables as additional latent variables.

Both for MNL models as well as for binomial models, this data augmentation
scheme leads to a conditionally Gaussian posterior and allows to sample all
unknowns through highly efficient block moves. Again, we apply iMDA to derive
an ultimate Pólya-Gamma sampler which is shown to mix well also for highly
unbalanced data.

The rest of the paper is structured as follows. Section 2 introduces the
ultimate Pólya-Gamma sampler for binary data. This sampling strategy is
extended to multinomial data in Section 3 and to binomial data in Section 4.
For all three types of models, the ultimate Pólya-Gamma sampler is compared to
alternative algorithms for simulated data in Section 5. Section 6 considers
binary state space models. Section 7 concludes.

##  2 Ultimate Pólya-Gamma samplers for binary data

###  2.1 Latent variable representations for binary data

Binary models for a set of $N$ binary data ${\mathbf{y}}=(y_{1},\ldots,y_{N})$
are defined by

$\mbox{\rm Pr}(y_{i}=1|\lambda_{i})=F_{\varepsilon}(\log\lambda_{i}),$ (1)

where $\lambda_{i}$ depends on exogenous variables and unknown parameters
$\boldsymbol{\beta}$, e.g.
$\log\lambda_{i}={\mathbf{x}}_{i}\boldsymbol{\beta}$ in a standard binary
regression model. Choosing the cdf
$F_{\varepsilon}(\varepsilon)=\Phi(\varepsilon)$ of the standard normal
distribution leads to the probit model $\mbox{\rm
Pr}(y_{i}=1|\lambda_{i})=\Phi(\log\lambda_{i})$, whereas the cdf
$F_{\varepsilon}(\varepsilon)=e^{\varepsilon}/(1+e^{\varepsilon})$ of the
logistic distribution leads to the logit model

$\mbox{\rm Pr}(y_{i}=1|\lambda_{i})=\lambda_{i}/(1+\lambda_{i}).$

A latent variable representation of model (1) involving the latent utility
$z_{i}$ is given by:

$\displaystyle y_{i}=I\\{z_{i}>0\\},$ $\displaystyle\displaystyle
z_{i}=\log\lambda_{i}+\varepsilon_{i},\quad\varepsilon_{i}\sim
f_{\varepsilon}(\varepsilon_{i}),$ (2)

where
$f_{\varepsilon}(\varepsilon)=F_{\varepsilon}^{\prime}(\varepsilon)=\phi(\varepsilon)$
is equal to the standard normal pdf for a probit model and equal to
$f_{\varepsilon}(\varepsilon)=e^{\varepsilon}/(1+e^{\varepsilon})^{2}$ for a
logit model.

While MCMC estimation based on (2) is straightforward for the probit model
using one level of data augmentation involving the latent utilities
$(z_{1},\ldots,z_{N})$ (Albert & Chib, 1993), for the logit model a second
level of data augmentation is required, based on a mixture representation of
the logistic distribution. In the present paper, we apply a new mixture
representation of the logistic distribution,

$\displaystyle
f_{\varepsilon}(\varepsilon)=e^{\varepsilon}/(1+e^{\varepsilon})^{2}=\frac{1}{4}\int
e^{-\omega\,\varepsilon^{2}/2}p(\omega)d\,\omega,$ (3)

where $\omega\sim\mathcal{PG}\left(2,0\right)$ follows a Pólya-Gamma
distribution (Polson et al., 2013) with parameter $b=2$ and $\kappa=0$, see
Appendix A.1.1 and A.1.2 for details. This representation is very convenient,
as the conditional posterior of $\omega_{i}$ given $\varepsilon_{i}$ is equal
to a tilted Pólya-Gamma distribution, i.e.
$\omega_{i}|\varepsilon_{i}\sim\mathcal{PG}\left(2,|\varepsilon_{i}|\right)$
which is easy to sample from, see Polson et al. (2013). Hence, whenever
$\log\lambda_{i}$ is linear the unknown parameters $\boldsymbol{\beta}$, such
as $\log\lambda_{i}={\mathbf{x}}_{i}\boldsymbol{\beta}$ for a binary
regression model, this new representation allows to construct a Pólya-Gamma
sampler that (a) samples the unknown parameters $\boldsymbol{\beta}$
conditional on the latent variables $(z_{1},\ldots,z_{N})$ and
$(\omega_{1},\ldots,\omega_{N})$ from a conditionally Gaussian model, and (b)
samples $z_{i},\omega_{i}|\lambda_{i},y_{i}$ independently for each $i$ by
sampling the latent variable $z_{i}$ from $z_{i}|\lambda_{i},y_{i}$ in the
original latent variable model (2) in the usual manner (see Appendix A.3.1)
and sampling the parameter $\omega_{i}$ conditional on $z_{i}$ from
$\omega_{i}|z_{i},\lambda_{i}\sim\mathcal{PG}\left(2,|z_{i}-\log\lambda_{i}|\right)$.

###  2.2 MCMC estimation

While the sampler outlined above is easy to implement, such data augmentation
schemes are known to be slowly mixing for imbalanced data sets where only a
few cases with $y_{i}=1$ or $y_{i}=0$ among the $N$ data points are observed,
see e.g. Johndrow et al. (2019). To deal with this issue, we combine our new
Pólya-Gamma mixture representation with marginal data augmentation (van Dyk &
Meng, 2001; Imai & van Dyk, 2005; McCulloch et al., 2000).

Given a global scale parameter $\delta$ and a threshold parameter $\gamma$, we
transform the latent variable $z_{i}$ such that

$\displaystyle\tilde{z}_{i}=\sqrt{\delta}z_{i}+\gamma,$ (4)

to define the following expanded model:

$\displaystyle\displaystyle
y_{i}=I\\{\tilde{z}_{i}>\gamma\\},\qquad\tilde{z}_{i}=\log\tilde{\lambda}_{i}+\sqrt{\delta}\varepsilon_{i}.$
(5)

Transformation (4) implies that
$\log\tilde{\lambda}_{i}=\sqrt{\delta}\log\lambda_{i}+\gamma$. For a binary
regression model with $x_{i,d}\equiv 1$, e.g.,
$\log\tilde{\lambda}_{i}={\mathbf{x}}_{i}\tilde{\boldsymbol{\beta}}$, where
$\tilde{\beta}_{j}=\sqrt{\delta}\beta_{j}+I\\{j=d\\}\gamma$.

The global scale parameter $\delta$ and the threshold parameter $\gamma$,
however, are not identified, since the marginal version of the expanded model
(5), where the latent variables $(\tilde{z}_{1},\ldots,\tilde{z}_{N})$ are
integrated out, is independent both of $\gamma$ and $\delta$:

$\displaystyle\mbox{\rm Pr}(y_{i}=1|\lambda_{i},\gamma,\delta)=\mbox{\rm
Pr}(\tilde{z}_{i}>\gamma|\lambda_{i},\gamma,\delta)=F_{\varepsilon}(\frac{\log\tilde{\lambda}_{i}-\gamma}{\sqrt{\delta}})=\mbox{\rm
Pr}(y_{i}=1|\lambda_{i}).$ (6)

Various MCMC schemes based on the expanded model (5) are possible. We consider
here marginal data augmentation where the following prior
$p(\gamma,\delta)=p(\gamma|\delta)p(\delta)$ is defined on the parameters
$\gamma$ and $\delta$ with fixed hyper parameters $d_{0},D_{0}$ and $G_{0}$:

$\displaystyle\gamma|\delta\sim\mathcal{N}\left(0,\delta
G_{0}\right),\quad\delta\sim\mathcal{G}^{-1}\left(d_{0},D_{0}\right).$ (7)

Note that the joint prior $p(\gamma,\delta)$ can also be written as
$p(\delta|\gamma)p(\gamma)$ where

$\displaystyle\gamma\sim t_{2d_{0}}\left(0,G_{0}^{\star}\right),\quad
G_{0}^{\star}=G_{0}D_{0}/d_{0},\quad\delta|\gamma\sim\mathcal{G}^{-1}\left(d_{0}+1/2,D_{0}+\gamma^{2}/(2G_{0})\right).$
(8)

The prior for the model parameters is defined by specifying a prior for the
parameter $\tilde{\boldsymbol{\beta}}$ in the expanded model (5) which is
allowed to depend on the global scale parameter $\delta$.

Let $\tilde{\boldsymbol{z}}$ summarize all latent variables in the model i.e.
$\tilde{\boldsymbol{z}}=(\tilde{z}_{1},\ldots,\tilde{z}_{N})$ for a probit
model and
$\tilde{\boldsymbol{z}}=(\tilde{z}_{1},\ldots,\tilde{z}_{N},\omega_{1},\ldots,\omega_{N})$
for a logit model. Marginal data augmentation is based on (a) sampling the
latent variables $\tilde{\boldsymbol{z}}$ from
$p(\tilde{\boldsymbol{z}}|\boldsymbol{\beta},{\mathbf{y}})$ and (b) sampling
the unknown parameters $\boldsymbol{\beta}$ from
$p(\boldsymbol{\beta}|\tilde{\boldsymbol{z}},{\mathbf{y}})$ without
conditioning on $\delta$ and $\gamma$.

Step (a) is carried out by sampling from the augmented posterior
$p(\tilde{\boldsymbol{z}},\delta,\gamma|\boldsymbol{\beta},{\mathbf{y}})=p(\tilde{\boldsymbol{z}}|\delta,\gamma,\boldsymbol{\beta},{\mathbf{y}})p(\delta,\gamma|\boldsymbol{\beta},{\mathbf{y}})$
by first drawing $(\delta,\gamma)$ from the conditional prior
$p(\delta,\gamma|\boldsymbol{\beta})$111Note that
$p(\delta,\gamma|\boldsymbol{\beta},{\mathbf{y}})\propto
p({\mathbf{y}}|\delta,\gamma,\boldsymbol{\beta})p(\delta,\gamma|\boldsymbol{\beta})\propto
p(\delta,\gamma|\boldsymbol{\beta})$ is independent of ${\mathbf{y}}$, see
again (6), but may depends on the model parameters $\boldsymbol{\beta}$., then
drawing $\tilde{\boldsymbol{z}}$ conditional on $(\delta,\gamma)$, and finally
discarding $(\delta,\gamma)$. To sample from
$\tilde{\boldsymbol{z}}|\delta,\gamma$, we sample $z_{i}$ for each $i$ from
$p(z_{i}|\lambda_{i},y_{i})$ in the identified model (2), sample
$\omega_{i}|z_{i},\lambda_{i}\sim\mathcal{PG}\left(2,|z_{i}-\log\lambda_{i}|\right)$
and define $\tilde{z}_{i}$ according to transformation (4). In each sweep of
the sampler, we are dealing with a different expanded model. This allows much
larger moves than in conventional sampling where $\delta\equiv 1$ and
$\gamma\equiv 0$ for all MCMC draws.

Step (b) is carried out by joint sampling from
$p(\boldsymbol{\beta},\delta,\gamma|\tilde{\boldsymbol{z}},{\mathbf{y}})$ and
dropping $\delta$ and $\gamma$. To this aim, we sample
$(\tilde{\boldsymbol{\beta}},\delta^{\mathrm{new}},\gamma^{\mathrm{new}})\sim
p(\tilde{\boldsymbol{\beta}},\delta^{\mathrm{new}},\gamma^{\mathrm{new}}|\tilde{\boldsymbol{z}},{\mathbf{y}})$
in the expanded model (5) and match back to (2) through the inverse of (4),
$\log\lambda_{i}=(\log\tilde{\lambda}_{i}-\gamma^{\mathrm{new}})/\sqrt{\delta^{\mathrm{new}}}$.
Using representation (8) of $p(\delta,\gamma)$, the joint posterior
$p(\tilde{\boldsymbol{\beta}},\delta,\gamma|\tilde{\boldsymbol{z}},{\mathbf{y}})$
factorizes for any binary model as
$p(\tilde{\boldsymbol{\beta}},\delta,\gamma|\tilde{\boldsymbol{z}},{\mathbf{y}})\propto
p({\mathbf{y}}|\gamma,\tilde{\boldsymbol{z}})p(\gamma)p(\delta|\gamma)p(\tilde{\boldsymbol{z}}|\tilde{\boldsymbol{\beta}},\delta)p(\tilde{\boldsymbol{\beta}}|\delta)$.
Thus, we first sample $\gamma$ from the marginal posterior
$p(\gamma|{\mathbf{y}},\tilde{\boldsymbol{z}})\propto
p({\mathbf{y}}|\gamma,\tilde{\boldsymbol{z}})p(\gamma)$ which is equal to the
Student-$t$ prior given in (8), however, truncated to a region depending on
the latent variables $\tilde{z}_{i}$:

$\displaystyle\displaystyle
p(\gamma|{\mathbf{y}},\tilde{\boldsymbol{z}})\propto
p(\gamma)\prod_{i:y_{i}=0}I\\{\gamma>\tilde{z}_{i}\\}\prod_{i:y_{i}=1}I\\{\gamma<\tilde{z}_{i}\\}\propto
p(\gamma)I\\{L(\tilde{\boldsymbol{z}})<\gamma<O(\tilde{\boldsymbol{z}})\\},$
(9)

where $L(\tilde{\boldsymbol{z}})=\max_{i:y_{i}=0}\tilde{z}_{i}$ and
$O(\tilde{\boldsymbol{z}})=\min_{i:y_{i}=1}\tilde{z}_{i}$, if
$0<\\#\\{i:y_{i}=1\\}<N$. If $\\#\\{i:y_{i}=0\\}=0$, then
$L(\tilde{\boldsymbol{z}})=-\infty$; if $\\#\\{i:y_{i}=1\\}=0$, then
$O(\tilde{\boldsymbol{z}})=+\infty$. Hence, $\gamma$ is sampled as:

$\displaystyle\gamma=\sqrt{G_{0}^{\star}}T_{2d_{0}}^{-1}\left(U\times
T_{2d_{0}}(\frac{O(\tilde{\boldsymbol{z}})}{\sqrt{G_{0}^{\star}}})+(1-U)\times
T_{2d_{0}}(\frac{L(\tilde{\boldsymbol{z}})}{\sqrt{G_{0}^{\star}}})\right),$
(10)

where $U\sim\mathcal{U}\left[0,1\right]$ and $T_{\nu}(\cdot)$ and
$T_{\nu}^{-1}(\cdot)$, are, the cdf and the quantile function of a
$t_{\nu}$-distribution.

Algorithm 1 The ultimate Pólya-Gamma sampler for binary data.

Choose starting values for
$\boldsymbol{\lambda}=(\lambda_{1},\ldots,\lambda_{N})$ and repeat the
following steps:

  * (a-1)

For each $i=1,\ldots,N$, sample
$z_{i}=\log\lambda_{i}+F_{\varepsilon}^{-1}(y_{i}+U_{i}(1-y_{i}-\pi_{i}))$ in
model (2), where $U_{i}\sim\mathcal{U}\left[0,1\right]$,
$\pi_{i}=F_{\varepsilon}(\log\lambda_{i})$, and
$F_{\varepsilon}^{-1}(p)=\Phi^{-1}(p)$ for the probit and
$F_{\varepsilon}^{-1}(p)=\log p-\log(1-p)$ for the logit model. For a logit
model, sample
$\omega_{i}|z_{i},\log\lambda_{i}\sim\mathcal{PG}\left(2,|z_{i}-\log\lambda_{i}|\right)$.

  * (a-2)

Sample $\delta^{\star}\sim\mathcal{G}^{-1}\left(d_{0},D_{0}\right)$ and
$\gamma^{\star}\sim\gamma|\delta^{\star},\boldsymbol{\lambda}$. For a binary
regression model, $\gamma|\delta^{\star},\beta_{d}$ is given in (15). Move to
the unidentified model by defining
$\tilde{z}_{i}=\sqrt{\delta^{\star}}z_{i}+\gamma^{\star}$ for $i=1,\ldots,N$.

  * (b-1)

Sample $\gamma^{\mathrm{new}}$ from the truncated $t$-prior $p(\gamma)$ using
(10). Sample
$\delta^{\mathrm{new}}\sim\delta^{\mathrm{new}}|\gamma^{\mathrm{new}},\tilde{\boldsymbol{z}}$.
For a binary regression model,
$\delta^{\mathrm{new}}|\gamma^{\mathrm{new}},\tilde{\boldsymbol{z}}\sim\mathcal{G}^{-1}\left(d_{N},D_{N}(\gamma^{\mathrm{new}})\right)$
with $d_{N}$ and $D_{N}(\gamma^{\mathrm{new}})$ given by (17).

  * (b-2)

Sample $\tilde{\boldsymbol{\beta}}$ from
$\tilde{\boldsymbol{\beta}}|\tilde{\boldsymbol{z}},\delta^{\mathrm{new}}$. For
a binary regression model,
$\tilde{\boldsymbol{\beta}}|\tilde{\boldsymbol{z}},\delta^{\mathrm{new}}\sim\mathcal{N}_{d}\left({\mathbf{b}}_{N},\delta^{\mathrm{new}}{\mathbf{B}}_{N}\right)$
where ${\mathbf{b}}_{N}$ and ${\mathbf{B}}_{N}$ are given by (18). Move back
to the identified model using
$\log\lambda_{i}=(\log\tilde{\lambda}_{i}-\gamma^{\mathrm{new}})/\sqrt{\delta^{\mathrm{new}}}$,
e.g. transformation (11) for a binary regression model and define
$z_{i}^{\mathrm{new}}=(\tilde{z}_{i}-\gamma^{\mathrm{new}})/\sqrt{\delta^{\mathrm{new}}}$.

Conditional on $\gamma$, we sample $(\tilde{\boldsymbol{\beta}},\delta)$ from
$p(\tilde{\boldsymbol{\beta}},\delta|\gamma,\tilde{\boldsymbol{z}})=p(\tilde{\boldsymbol{\beta}}|\tilde{\boldsymbol{z}},\delta)p(\delta|\tilde{\boldsymbol{z}},\gamma)$.
The specifics depend on the chosen model for $\lambda_{i}$. We consider now in
detail a binary regression model, where
$\log\lambda_{i}={\mathbf{x}}_{i}\boldsymbol{\beta}$ and $x_{i,d}\equiv 1$
corresponds to the intercept. In the expanded model (5), the linear predictor
is $\log\tilde{\lambda}_{i}={\mathbf{x}}_{i}\tilde{\boldsymbol{\beta}}$, where
$\tilde{\beta}_{j}=\sqrt{\delta}\beta_{j}+I\\{j=d\\}\gamma$. We specify the
prior
$\tilde{\boldsymbol{\beta}}|\delta\sim\mathcal{N}_{d}\left({\mathbf{0}},\delta\boldsymbol{A}_{0}\right)$
in the expanded model. Using the inverse of (4),

$\displaystyle\beta_{j}=\tilde{\beta}_{j}/\sqrt{\delta},j=1,\ldots,d-1,\quad\beta_{d}=(\tilde{\beta}_{d}-\gamma)/\sqrt{\delta},$
(11)

the conditional prior $p(\boldsymbol{\beta}|\gamma,\delta)$ can be derived:
222The marginal prior in the identified model reads
$\boldsymbol{\beta}\sim\mathcal{N}_{d}\left(0,\boldsymbol{A}_{0}+G_{0}\boldsymbol{e}_{d}\boldsymbol{e}_{d}^{\top}\right)$
where
$\boldsymbol{e}_{d}=\left(\begin{array}[]{cc}{\mathbf{0}}_{1\times(d-1)}&\end{array}\right)^{\top}$.

$\displaystyle\boldsymbol{\beta}|\gamma,\delta\sim\mathcal{N}_{d}\left(\boldsymbol{a}_{0},\boldsymbol{A}_{0}\right),\qquad\boldsymbol{a}_{0}=\left(\begin{array}[]{c}{\mathbf{0}}_{(d-1)\times
1}\\\ -\frac{\gamma}{\sqrt{\delta}}\end{array}\right).$ (14)

From (14) we obtain that $p(\delta,\gamma|\boldsymbol{\beta})\propto
p(\beta_{d}|\delta,\gamma)p(\gamma|\delta)p(\delta)$ is proportional to
$p(\gamma|\delta,\beta_{d})p(\delta)$ where $p(\delta)$ is given by (7) and
$p(\gamma|\delta,\beta_{d})$ is equal to the following normal distribution:

$\displaystyle\gamma|\delta,\beta_{d}\sim\mathcal{N}\left(\sqrt{\delta}g_{1},\delta
G_{1}\right),\quad G_{1}=G_{0}A_{0,dd}/(G_{0}+A_{0,dd}),\quad
g_{1}=-\beta_{d}\times G_{0}/(G_{0}+A_{0,dd}).$ (15)

Hence, to sample from the conditional prior
$p(\delta,\gamma|\boldsymbol{\beta})$, $\delta$ is sampled from
$\mathcal{G}^{-1}\left(d_{0},D_{0}\right)$ and $\gamma|\delta,\beta_{d}$ is
sampled from $\mathcal{N}\left(\sqrt{\delta}g_{1},\delta G_{1}\right)$.

Next, for each $i$ the latent variables $(z_{i},\omega_{i})$ are sampled from
the identified model (2) and transformation (4) yields the latent variables
$(\tilde{z}_{i},\omega_{i})$ in the expanded model

$\displaystyle
y_{i}=I\\{\tilde{z}_{i}>\gamma\\},\qquad\tilde{z}_{i}={\mathbf{x}}_{i}\tilde{\boldsymbol{\beta}}+\sqrt{\delta}\varepsilon_{i}.$
(16)

Discarding $\gamma$ and $\delta$ completes Step (a).

Step (b) is carried out in the expanded model (16) conditional on all
$(\tilde{z}_{i},\omega_{i})$s. A new value of $\gamma$ is sampled from the
marginal posterior (10). The conditionally conjugate prior for
$\tilde{\boldsymbol{\beta}}|\delta$ allows joint sampling of
$\tilde{\boldsymbol{\beta}},\delta|\tilde{\boldsymbol{z}},\gamma$. The
marginal posterior
$\delta|\tilde{\boldsymbol{z}},\gamma\sim\mathcal{G}^{-1}\left(d_{N},D_{N}(\gamma)\right)$
is inverted Gamma with following moments:

$\displaystyle d_{N}=d_{0}+{\frac{1}{2}}+\frac{N}{2},$ (17) $\displaystyle
D_{N}(\gamma)=D_{0}+{\frac{\gamma^{2}}{2G_{0}}}+\frac{1}{2}\sum_{i=1}^{N}\omega_{i}(\tilde{z}_{i}-{\mathbf{x}}_{i}{\mathbf{b}}_{N})^{2}+\frac{1}{2}{\mathbf{b}}_{N}^{\top}\boldsymbol{A}_{0}^{-1}{\mathbf{b}}_{N},$

where

$\displaystyle{\mathbf{b}}_{N}={\mathbf{B}}_{N}(\sum_{i=1}^{N}\omega_{i}{\mathbf{x}}_{i}^{\top}\tilde{z}_{i}),\quad{\mathbf{B}}_{N}=(\boldsymbol{A}_{0}^{-1}+\sum_{i=1}^{N}\omega_{i}{\mathbf{x}}_{i}^{\top}{\mathbf{x}}_{i})^{-1},$
(18)

while the conditional posterior
$\tilde{\boldsymbol{\beta}}|\delta,\tilde{\boldsymbol{z}},\gamma\sim\mathcal{N}_{d}\left({\mathbf{b}}_{N},\delta{\mathbf{B}}_{N}\right)$
is multivariate normal. The posterior draw of $\tilde{\boldsymbol{\beta}}$ is
transformed back to $\boldsymbol{\beta}$ using (11).

The ultimate Pólya-Gamma (UPG) sampler for binary data is summarized in
Algorithm 1.333All algorithms outlined in this work will soon be made
available in R package UPG. Note that for the probit model the marginal data
augmentation scheme of Imai & van Dyk (2005) would result, if $\gamma\equiv 0$
was fixed to zero. Our data augmentation scheme is much more general than
that. First of all, it works for any binary data model with a latent variable
representation. Second, freeing the location of the threshold $\gamma$ in
model (5) restores balance, with $\gamma$ being considerably larger than zero,
if only a few cases with $y_{i}=1$ are present among the $N$ binary
observations, as will be illustrated in Section 2.3. The imbalanced marginal
data augmentation (iMDA) scheme introduced in this section leads to an MCMC
scheme that is well mixing, even in cases of extremely unbalanced data, see
Section 5 for further illustration.

###  2.3 Illustration & intuition

Figure 1: MCMC draws and autocorrelation function of $\beta_{d}$ for the UPG
logit sampler (top) and the original Pólya-Gamma logit sampler (bottom); 2 out
of $N=10,000$ observations are non-zero.

Figure 2: Posterior densities for latent utilities in the unboosted model
(first column), $\gamma$ (second column), $\delta$ (third column) and boosted
utilities (fourth column). The rows correspond to different data sets with
$N=1,000$ and success probabilities of 1% (first row), 50% (second row) and
99% (third row). The figure demonstrates how the auxiliary parameters $\gamma$
and $\delta$ transform the latent utilities $z_{i}$ to pseudo-balanced,
boosted utilities $\tilde{z}_{i}$, allowing for efficient estimation.

As a first illustration of the potential merits of the UPG approach, we
compare estimation efficiency of the popular Pólya-Gamma sampler from Polson
et al. (2013) with UPG in an imbalanced data setting in Figure 1. It is
clearly visible that, while both samplers simulate from the same posterior
distribution, the UPG approach outperforms the conventional Pólya-Gamma
sampler considerably in terms of efficiency. Note that this efficiency gain is
realized in spite of an additional layer of data augmentation to sample the
latent utilities $z_{i}$. While this second layer of data augmentation allows
for straight-forward model extensions due to a conditionally Gaussian
likelihood representation, it usually increases autocorrelation in the
posterior draws significantly. This is counteracted by our novel MCMC boosting
strategy based on the working parameters $\gamma$ and $\delta$.

The interplay of these two working parameters is illustrated for the case of a
logistic regression model in Figure 2. In case of imbalanced data where only a
few data points are a success, most of the latent utilities $z_{i}$ will be
negative. On the contrary, $z_{i}$ will mostly be greater than zero if the
majority of data points are successes. In the balanced case, $z_{i}$ is
concentrated around zero. The role of $\gamma$ is to shift the mean of the
latent utilities towards zero in order to restore balance in the unidentified
model. At the same time, $z_{i}$ is rescaled using $\delta$. Ultimately, this
results in an estimation setup where the boosted utilities
$\tilde{z_{i}}=z_{i}\sqrt{\delta}+\gamma$ are re-centered and re-scaled such
that the resulting distribution exhibits fat tails and a zero mean. This leads
to an improvement in the signal-to-noise ratio in the boosted model when
compared to the unboosted model, effectively allowing for more efficient MCMC
estimation.

In addition to this intuitive explanation of our boosting strategy, the role
of $\delta$ and $\gamma$ can be explored in a more technical context as well.
Variants of the working parameter $\delta$ have been discussed for instance in
Imai & van Dyk (2005). Put briefly, $\delta$ is used to ensure that the
posterior of the boosted utilities $\tilde{z_{i}}$ becomes more diffuse than
the posterior of the unboosted utilities $z_{i}$. Through this diffusion,
autocorrelation in the posterior draws will decrease as $\tilde{z_{i}}$
becomes less dependent on the current values of the remaining model parameters
(e.g. the regression coefficients $\boldsymbol{\beta}$). This allows the
sampler to move through the parameter space more quickly. A formal proof that
this strategy improves convergence of data augmentation algorithms can be
found in Meng & Van Dyk (1999). Theoretical results on MCMC sampling
efficiency from the perspective of computational complexity are provided in
Johndrow et al. (2019). One key insight is that data augmentation strategies
in imbalanced categorical data sets results in a sampling behavior where the
average step size of the sampler is too small to explore the whole posterior
distribution efficiently. In our boosting strategy, this issue is effectively
offset by restoring pseudo-balance through the working parameter $\gamma$,
allowing the UPG sampler to efficiently explore the entire posterior
distribution.

##  3 Ultimate Pólya-Gamma samplers for multinomial data

Let $\left\\{y_{i}\right\\}$ be a sequence of categorical data,
$i=1,\ldots,N$, where $y_{i}$ is equal to one of at least three unordered
categories. The categories are labeled by $L=\\{0,\ldots,m\\}$, and for any
$k$ the set of all categories but $k$ is denoted by
$L_{-k}=L\setminus\\{k\\}$. We assume that the observations are mutually
independent and that for each $k\in L$ the probability of $y_{i}$ taking the
value $k$ depends on covariates ${\mathbf{x}}_{i}$ in the following way:

$\displaystyle\mbox{\rm
Pr}(y_{i}=k|\boldsymbol{\beta}_{0},\ldots,\boldsymbol{\beta}_{m})=\pi_{ki}(\boldsymbol{\beta}_{0},\ldots,\boldsymbol{\beta}_{m})=\displaystyle\frac{\exp({\mathbf{x}}_{i}\boldsymbol{\beta}_{k})}{\displaystyle\sum_{l=0}^{m}\exp({\mathbf{x}}_{i}\boldsymbol{\beta}_{l})},$
(19)

where $\boldsymbol{\beta}_{0},\ldots,\boldsymbol{\beta}_{m}$ are category
specific unknown parameters of dimension $d$. To make the model identifiable,
the parameter $\boldsymbol{\beta}_{k_{0}}$ of a baseline category $k_{0}$ is
set equal to ${\mathbf{0}}$: $\boldsymbol{\beta}_{k_{0}}={\mathbf{0}}$. Thus
the parameter $\boldsymbol{\beta}_{k}$ is in terms of the change in log-odds
relative to the baseline category $k_{0}$. In the following, we assume without
loss of generality that $k_{0}=0$. A more general version of the multinomial
logit (MNL) model (again with baseline $k_{0}=0$) reads:

$\displaystyle\mbox{\rm
Pr}(y_{i}=k|\boldsymbol{\beta})=\displaystyle\lambda_{ki}/(1+\sum_{l=1}^{m}\lambda_{li}),$
(20)

where $\lambda_{1i},\ldots,\lambda_{mi}$ depend on unknown parameters
$\boldsymbol{\beta}$. For the standard MNL regression model (19), for
instance, $\log\lambda_{ki}={\mathbf{x}}_{i}\boldsymbol{\beta}_{k}$ for
$k=1,\ldots,m$.

In this section, we consider a new random utility representation of a MNL
model with at least 3 categories for the general formulation (20). Starting
point is writing the MNL model as a random utility model (RUM), see McFadden
(1974):

$\displaystyle\displaystyle{u}_{ki}=\log\lambda_{ki}+\epsilon_{ki},\quad
k=0,\ldots,m,$ (21) $\displaystyle\displaystyle
y_{i}=k\Leftrightarrow{u}_{ki}=\max_{l\in L}{u}_{li}.$ (22)

Thus the observed category is equal to the category with maximal utility. If
the errors $\epsilon_{0i},\ldots,\epsilon_{mi}$ in (21) are i.i.d. random
variables from an extreme value ($\mathrm{EV}$) distribution, then the MNL
model (20) results as marginal distribution of the categorical variable
$y_{i}$.

To sample the category specific parameters in $\log\lambda_{ki}$ conditional
on the remaining parameters, it is common practice to reduce the RUM (22) to a
binary logit model with observations $d_{ki}=I\\{y_{i}=k\\}$ (Frühwirth-
Schnatter & Frühwirth, 2010). However, the UPG sampler introduced in Section 2
cannot be applied, since the parameters for all categories in $L_{-k}$ appear
as an off-set in the latent variable representation.

Our new data augmentation scheme resolves this problem by reducing the data to
three categories: category $k$, the baseline category, and a category which
collapses the remaining categories in $A=\\{l\in L|l\neq\\{k,0\\}\\}$. For all
categories in $A$, we define an aggregated utility ${u}_{a,i}=\max_{l\in
A}{u}_{li}$ as the maximal utility in $A$. It can be shown that

$\displaystyle{u}_{a,i}=\log\lambda_{a,i}+\epsilon_{a,i},\quad\lambda_{a,i}=\sum_{l\in
A}\lambda_{li},$

where $\epsilon_{a,i}$ follows an $\mathrm{EV}$ distribution. The proof is
straightforward since
$\exp(-{u}_{li})\sim\mathcal{E}\left(\lambda_{li}\right)$ and therefore
$\exp\left(-{u}_{a,i}\right)=\min_{l\in
A}\exp(-{u}_{li})\sim\mathcal{E}\left(\lambda_{a,i}\right)$. Evidently, for
$m=1$, $A$ is the empty set as there are no alternative categories and we
obtain the binary logit model.

For each specific category $k=1,\ldots,m$, the full RUM (21) can be
represented as following aggregated RUM, based on the latent utilities
$({u}_{ki},{u}_{0i},{u}_{a,i})$:

$\displaystyle\displaystyle{u}_{ki}=\log\lambda_{ki}+\epsilon_{ki},$ (23)
$\displaystyle\displaystyle{u}_{0i}=\epsilon_{0i},$
$\displaystyle\displaystyle{u}_{a,i}=\log\lambda_{a,i}+\epsilon_{a,i},$ (24)
$\displaystyle\displaystyle
y_{i}=\left\\{\begin{array}[]{ll}k,&{u}_{ki}\geq\max({u}_{0i},{u}_{a,i}),\\\
0,&{u}_{0i}\geq\max({u}_{ki},{u}_{a,i}),\\\ a\in
A,&{u}_{a,i}\geq\max({u}_{ki},{u}_{0i}),\end{array}\right.$ (28)

where the errors $\epsilon_{ki},\epsilon_{0i},\epsilon_{a,i}$ are iid from an
$\mathrm{EV}$ distribution.

To derive an efficient sampler to estimate the category specific parameters in
$\log\lambda_{ki}$ conditional on ${u}_{a,i}$, it is useful to rewrite (23) in
the following way:

$\displaystyle{u}_{ki}-{u}_{0i}=\log\lambda_{ki}+\varepsilon_{ki},$ (29)

where $\varepsilon_{ki}\sim\mathcal{LO}$ follows a logistic distribution,
independent of $\epsilon_{a,i}$, and the choice equation is the same as in
(28). We term this model the aggregated RUM representation.

Conditional on having observed $y_{i}$, the posterior distribution
$p({\mathbf{{u}}}_{i}|\lambda_{ki},y_{i})$ of the latent variables
${\mathbf{{u}}}_{i}=({u}_{ki},{u}_{0i},{u}_{a,i})$ is of closed form and easy
to sample from, see Theorem 1 which is proven in Appendix A.2.

###### Theorem 1

Draws from the posterior distribution
$p({\mathbf{{u}}}_{i}|\lambda_{ki},y_{i})$ can be represented as:

$\displaystyle\left(\begin{array}[]{c}e^{-{u}_{ki}}\\\ e^{-{u}_{0i}}\\\
e^{-{u}_{a,i}}\\\
\end{array}\right)=-\frac{\log(U_{i})}{1+\sum_{l=1}^{m}\lambda_{li}}\left(\begin{array}[]{c}1\\\
1\\\ 1\\\ \end{array}\right)-\left(\begin{array}[]{c}\displaystyle\frac{\log
V_{ki}}{\lambda_{ki}}I\\{y_{i}\neq k\\}\\\\[2.84526pt] \log
V_{0i}\,I\\{y_{i}\neq 0\\}\\\\[2.84526pt]
\displaystyle\frac{\log(V_{a,i})}{\lambda_{a,i}}I\\{y_{i}\notin A\\}\\\
\end{array}\right),$ (39)

where $U_{i},V_{ki},V_{0i},V_{a,i}$, $i=1,\ldots,N$, are iid uniform random
numbers.

For the standard MNL regression model, the parameters $\boldsymbol{\beta}_{k}$
corresponding to category $k$ appears only in a linear regression model with
logistic errors, given by (29). The term $\log\lambda_{a,i}$ in (24) depends
only on the regression parameters $\boldsymbol{\beta}_{-k}$ of all other
categories in $A$ and, as opposed to the original partial dRUM,
$\log\lambda_{a,i}$ does not appear as an offset in the latent equation (23),
from which we estimate the category specific parameters
$\boldsymbol{\beta}_{k}$. Rather, the information of all other categories is
sufficiently summarized in the latent variable ${u}_{a,i}$ which affects only
the choice equation. Conditional on the aggregated utility ${u}_{a,i}$, the
logit model (29) is independent of $\boldsymbol{\beta}_{-k}$ and we can
proceed similarly as in Section 2 to derive an ultimate Pólya-Gamma sampler by
exploiting the Pólya-Gamma mixture representation of the logistic distribution
in (29) with latent variables $\omega_{ki}$, $i=1,\dots,N$.

To protect our sampler against imbalanced data, we apply iMDA with category
specific parameters $\gamma_{k}$ and $\delta_{k}$, using the transformations

$\displaystyle\tilde{u}_{ki}=\sqrt{\delta_{k}}{u}_{ki}+\gamma_{k},\quad\tilde{u}_{0i}=\sqrt{\delta_{k}}{u}_{0i},\quad\tilde{u}_{a,i}=\sqrt{\delta_{k}}{u}_{a,i}.$
(40)

Transformation (40) implies that
$\log\tilde{\lambda}_{ki}=\sqrt{\delta_{k}}\log\lambda_{ki}+\gamma_{k}$. For a
MNL regression model, e.g., it follows that
$\log\tilde{\lambda}_{ki}={\mathbf{x}}_{i}\tilde{\boldsymbol{\beta}}_{k}$,
where
$\tilde{\boldsymbol{\beta}}_{k,j}=\sqrt{\delta_{k}}\boldsymbol{\beta}_{k,j}+I\\{j=d\\}\gamma_{k}$.

Transformation (40) leads to following unidentified model:

$\displaystyle\tilde{u}_{ki}-\tilde{u}_{0i}=\log\tilde{\lambda}_{ki}+\sqrt{\delta_{k}}\varepsilon_{ki},\quad\varepsilon_{ki}\sim\mathcal{LO},$
(41) $\displaystyle\displaystyle
y_{i}=\left\\{\begin{array}[]{lll}k,&\tilde{u}_{ki}\geq\tilde{u}_{a,i}+\gamma_{k},&\tilde{u}_{ki}\geq\tilde{u}_{0i}+\gamma_{k},\\\
0,&\tilde{u}_{0i}\geq\tilde{u}_{ki}-\gamma_{k},&\tilde{u}_{0i}\geq\tilde{u}_{a,i},\\\
a\in
A,&\tilde{u}_{a,i}\geq\tilde{u}_{ki}-\gamma_{k},&\tilde{u}_{a,i}\geq\tilde{u}_{0i}.\\\
\end{array}\right.$ (45)

The priors for the parameters $\delta_{k}$ and $\gamma_{k}$ are chosen as for
the binary model:

$\displaystyle\gamma_{k}|\delta_{k}\sim\mathcal{N}\left(0,\delta_{k}G_{0}\right),\qquad\delta_{k}\sim\mathcal{G}^{-1}\left(d_{0},D_{0}\right).$

With $G_{0}^{\star}=G_{0}D_{0}/d_{0}$, the joint prior
$p(\delta_{k},\gamma_{k})$ of $\gamma_{k}$ and $\delta_{k}$ factors as:

$\displaystyle\gamma_{k}\sim
t_{2d_{0}}\left(0,G_{0}^{\star}\right),\quad\delta_{k}|\gamma_{k}\sim\mathcal{G}^{-1}\left(d_{0}+1/2,D_{0}+\gamma_{k}^{2}/(2G_{0})\right).$
(46)

Comparable to the binary case, we define a prior for
$\tilde{\boldsymbol{\beta}}_{k}$ in the boosted model (41). For a MNL
regression model, e.g., we assume the prior
$\tilde{\boldsymbol{\beta}}_{k}|\delta_{k}\sim\mathcal{N}_{d}\left({\mathbf{0}},\delta_{k}\boldsymbol{A}_{0}\right)$.
It is possible to derive explicitly the conditional prior
$\boldsymbol{\beta}_{k}|\gamma_{k},\delta_{k}\sim\mathcal{N}_{d}\left(-\gamma_{k}/\sqrt{\delta_{k}}\boldsymbol{e}_{d},\boldsymbol{A}_{0}\right)$
and the marginal prior
$\boldsymbol{\beta}_{k}\sim\mathcal{N}_{d}\left(0,\boldsymbol{A}_{0}+G_{0}\boldsymbol{e}_{d}\boldsymbol{e}_{d}^{\top}\right)$
in, respectively, model (41) and (29).

Algorithm 2 The ultimate Pólya-Gamma sampler for multinomial data.

Choose starting values for
$\boldsymbol{\lambda}=(\boldsymbol{\lambda}_{1},\ldots,\boldsymbol{\lambda}_{m})$,
$\boldsymbol{\lambda}_{k}=(\lambda_{k1},\ldots,\lambda_{kN})$ and repeat for
$k=1,\ldots,m$:

  * (a-1)

For each $i=1,\ldots,N$, sample the latent variables
${\mathbf{{u}}}_{ki}=({u}_{ki},{u}_{0i},{u}_{a,i})\sim
p({\mathbf{{u}}}_{ki}|\lambda_{ki},y_{i})$ in the aggregated RUM model (23) of
category $k$ using Theorem 1. Sample
$\omega_{ki}|\lambda_{ki},{\mathbf{{u}}}_{ki}\sim\mathcal{PG}\left(2,|{u}_{ki}-{u}_{0i}-\log\lambda_{ki}|\right)$
in the model (29).

  * (a-2)

Sample $\delta_{k}^{\star}\sim\mathcal{G}^{-1}\left(d_{0},D_{0}\right)$ and
$\gamma_{k}^{\star}\sim\gamma_{k}|\delta_{k}^{\star},\boldsymbol{\lambda}_{k}$.
For a MNL regression model,
$\gamma_{k}|\delta_{k}^{\star},\beta_{k,d}\sim\mathcal{N}\left(\sqrt{\delta^{\star}_{k}}\cdot
g_{1},\delta_{k}^{\star}G_{1}\right)$ with
$g_{1}=-\beta_{k,d}G_{0}/(G_{0}+A_{0,dd})$ and
$G_{1}=G_{0}A_{0,dd}/(G_{0}+A_{0,dd}).$ Define
$\tilde{u}_{ki}=\sqrt{\delta_{k}^{\star}}{u}_{ki}+\gamma_{k}$,
$\tilde{u}_{0i}=\sqrt{\delta_{k}^{\star}}{u}_{0i}$, and
$\tilde{u}_{a,i}=\sqrt{\delta_{k}^{\star}}{u}_{a,i}$ for $i=1,\ldots,N$, to
move to the unidentified model (41).

  * (b-1)

Sample $\gamma_{k}^{\mathrm{new}}$ from the truncated $t$-prior
$p(\gamma_{k})$ using (49). Sample
$\delta_{k}^{\mathrm{new}}\sim\delta_{k}^{\mathrm{new}}|\gamma_{k}^{\mathrm{new}},\tilde{\boldsymbol{z}}$.
For the MNL regression model,
$\delta_{k}^{\mathrm{new}}|\gamma_{k}^{\mathrm{new}},\tilde{\boldsymbol{z}}\sim\mathcal{G}^{-1}\left(d_{0}+(N+1)/2,D_{k}(\gamma_{k}^{\mathrm{new}})\right)$,
where

$\displaystyle
D_{k}(\gamma_{k})=D_{0}+\frac{\gamma_{k}^{2}}{2G_{0}}+\frac{1}{2}\sum_{i=1}^{N}\omega_{ki}(\tilde{u}_{ki}-\tilde{u}_{0i}-{\mathbf{x}}_{i}{\mathbf{b}}_{k})^{2}+\frac{1}{2}{\mathbf{b}}_{k}^{\top}\boldsymbol{A}_{0}^{-1}{\mathbf{b}}_{k},$
$\displaystyle{\mathbf{b}}_{k}={\mathbf{B}}_{k}(\sum_{i=1}^{N}\omega_{ki}{\mathbf{x}}_{i}^{\top}(\tilde{u}_{ki}-\tilde{u}_{0i})),\quad{\mathbf{B}}_{k}=(\boldsymbol{A}_{0}^{-1}+\sum_{i=1}^{N}\omega_{ki}{\mathbf{x}}_{i}^{\top}{\mathbf{x}}_{i})^{-1}.$
(47)

  * (b-2)

Sample from
$\tilde{\boldsymbol{\beta}}_{k}|\tilde{\boldsymbol{z}}_{k},\delta_{k}^{\mathrm{new}}$.
For the MNL regression model,
$\tilde{\boldsymbol{\beta}}_{k}|\tilde{\boldsymbol{z}}_{k},\delta_{k}^{\mathrm{new}}\sim\mathcal{N}_{d}\left({\mathbf{b}}_{k},\delta_{k}^{\mathrm{new}}{\mathbf{B}}_{k}\right)$
where ${\mathbf{b}}_{k}$ and ${\mathbf{B}}_{k}$ are given by (47). Move back
to the full RUM (22) using
$\log\lambda_{ki}=(\log\tilde{\lambda}_{ki}-\gamma_{k}^{\mathrm{new}})/\sqrt{\delta_{k}^{\mathrm{new}}}$
and discard $\tilde{\boldsymbol{z}}$. For the MNL regression model, e.g.
$\beta_{kj}=(\tilde{\beta}_{kj}-I\\{j=d\\}\gamma_{k}^{\mathrm{new}})/\sqrt{\delta_{k}^{\mathrm{new}}}$,
$j=1,\ldots,d$.

Based on (41), Algorithm 1 can be extended in a fairly straightforward manner.
The ultimate Pólya-Gamma sampler for multinomial data is summarized in
Algorithm 2, with latent variables
$\tilde{\boldsymbol{z}}_{k}=(\tilde{\mathbf{{u}}}_{1},\ldots,\tilde{\mathbf{{u}}}_{N},\omega_{k1},\ldots,\omega_{kN})$,
where
$\tilde{\mathbf{{u}}}_{i}=(\tilde{u}_{ki},\tilde{u}_{0i},\tilde{u}_{a,i})$.
The choice equation (45) implies following constraint for $\gamma_{k}$,
conditional on $\tilde{\boldsymbol{z}}$:

$\displaystyle L(\tilde{\boldsymbol{z}})<\gamma_{k}\leq
O(\tilde{\boldsymbol{z}}),\quad
O(\tilde{\boldsymbol{z}})=\min_{i:y_{i}=k}(\tilde{u}_{ki}-\max(\tilde{u}_{0i},\tilde{u}_{a,i}))$
(48) $\displaystyle\displaystyle
L(\tilde{\boldsymbol{z}})=\max\left(\max_{i:y_{i}=0}(\tilde{u}_{ki}-\tilde{u}_{0i}),\max_{i:y_{i}\in
A}(\tilde{u}_{ki}-\tilde{u}_{a,i})\right).$

If $\\#\\{y_{i}=k\\}=0$, then $O(\tilde{\boldsymbol{z}})=+\infty$; if
$\\#\\{y_{i}=k\\}=N$, then $L(\tilde{\boldsymbol{z}})=-\infty$. Hence, the
marginal posterior
$p(\gamma_{k}|{\mathbf{y}},\tilde{\boldsymbol{z}}_{k})\propto
p({\mathbf{y}}|\gamma_{k},\tilde{\boldsymbol{z}}_{k})p(\gamma_{k})$ is equal
to the Student-$t$ prior given in (46), however, truncated to a region
depending on $\tilde{\boldsymbol{z}}$:

$\displaystyle\displaystyle
p(\gamma_{k}|{\mathbf{y}},\tilde{\boldsymbol{z}}_{k})$ $\displaystyle\propto$
$\displaystyle
p(\gamma_{k})\prod_{i:y_{i}=k}I\\{\gamma_{k}<\tilde{u}_{ki}-\max(\tilde{u}_{0i},\tilde{u}_{a,i})\\}\prod_{i:y_{i}=0}I\\{\gamma_{k}>\tilde{u}_{ki}-\tilde{u}_{0i}\\}$
$\displaystyle\qquad\qquad\times\prod_{i:y_{i}\in
A}I\\{\gamma_{k}>\tilde{u}_{ki}-\tilde{u}_{a,i}\\}\propto
p(\gamma_{k})I\\{L(\tilde{\boldsymbol{z}})\leq\gamma_{k}\leq
O(\tilde{\boldsymbol{z}})\\}.$

Given $U\sim\mathcal{U}\left[0,1\right]$, $\gamma_{k}$ is sampled as:

$\displaystyle\gamma_{k}^{\mathrm{new}}=\sqrt{G_{0}^{\star}}T_{2d_{0}}^{-1}\left(U\times
T_{2d_{0}}\left(\frac{O}{\sqrt{G_{0}^{\star}}}\right)+(1-U)\times
T_{2d_{0}}\left(\frac{L}{\sqrt{G_{0}^{\star}}}\right)\right).$ (49)

Conditional on $\gamma_{k}$, we sample from
$\delta_{k}|\tilde{\boldsymbol{z}}_{k},\gamma_{k}$ and
$\tilde{\boldsymbol{\beta}}_{k}|\tilde{\boldsymbol{z}}_{k},\delta_{k}$, see
Step (b) in Algorithm 2.

##  4 Ultimate Pólya-Gamma samplers for binomial data

In this section, we consider models with binomial outcomes, i.e. models of the
form

$\displaystyle y_{i}\sim\mbox{\rm
BiNom}\left(N_{i},\pi_{i}\right),\qquad\mbox{\rm
logit}\,\pi_{i}=\log\lambda_{i},\qquad i=1,\ldots,N,$ (50)

with $\log\lambda_{i}={\mathbf{x}}_{i}\boldsymbol{\beta}$ for a binomial
regression model. As shown by Johndrow et al. (2019), Bayesian inference for
binomial regression models based on the Pólya Gamma sampler (Polson et al.,
2013) is sensitive to imbalanced data. A latent variable representation of
binomial models introduced by Fussl et al. (2013) is also is sensitive to
unbalanced data, as shown in Section 5. As for a logit model, which results
for $N_{i}\equiv 1$, iMDA would be an option to boost MCMC, however, Fussl et
al. (2013) provide no explicit choice equation, which is needed for iMDA. The
goal of this section is to define an ultimate Pólya-Gamma sampler for binary
outcomes by combing a new latent variable representation of binomial models
based on a Pólya-Gamma representation of the generalized logistic distribution
with iMDA to protect the algorithm against imbalanced data.

###  4.1 A new latent variable representations for binomial data

We introduce in Theorem 2 a new latent variable representation for binomial
outcomes involving latent variables $\boldsymbol{z}_{i}$ such that the latent
variable equation is linear in $\log\lambda_{i}$ and the errors follows a
generalized logistic distribution. An explicit choice equation is provided
which relates the latent variables
$\boldsymbol{z}_{1},\ldots,\boldsymbol{z}_{N}$ to observations
$y_{1},\ldots,y_{N}$. We show in Theorem 3 that, conditional on $y_{i}$, the
posterior distribution $\boldsymbol{z}_{i}|y_{i},\lambda_{i}$ of the latent
variables $\boldsymbol{z}_{i}$ is of closed form and easy to sample from. See
Appendix A.2 for a proof of Theorem 2 and 3.

###### Theorem 2

Latent variable representation of a binomial model For $0<y_{i}<N_{i}$, the
binomial model has following random utility representation:

$\displaystyle{w}_{i}=\log\lambda_{i}+\varepsilon_{{w},i},\quad\varepsilon_{{w},i}\sim{\mathcal{GL}_{\mbox{\tiny
II}}\left(k\right)},$ (51)
$\displaystyle{v}_{i}=\log\lambda_{i}+\varepsilon_{{v},i},\quad\varepsilon_{{v},i}\sim{\mathcal{GL}_{\mbox{\tiny
I}}\left(N_{i}-k\right)},$ $\displaystyle
y_{i}=k\Leftrightarrow{w}_{i}>0,{v}_{i}<0,$

where $\mathcal{GL}_{\mbox{\tiny I}}\left(\nu\right)$ and
$\mathcal{GL}_{\mbox{\tiny II}}\left(\nu\right)$ are, respectively, the
generalized logistic distributions of type I and type II. For $y_{i}=0$, the
model reduces to

$\displaystyle{v}_{i}=\log\lambda_{i}+\varepsilon_{{v},i},\quad\varepsilon_{{v},i}\sim{\mathcal{GL}_{\mbox{\tiny
I}}\left(N_{i}\right)},\quad y_{i}=0\Leftrightarrow{v}_{i}<0.$

For $y_{i}=N_{i}$, the model reduces to

$\displaystyle{w}_{i}=\log\lambda_{i}+\varepsilon_{{w},i},\quad\varepsilon_{{w},i}\sim{\mathcal{GL}_{\mbox{\tiny
II}}\left(N_{i}\right)},\quad y_{i}=N_{i}\Leftrightarrow{w}_{i}>0.$

For $N_{i}=1$, the logistic model results, as both $\mathcal{GL}_{\mbox{\tiny
I}}\left(\nu\right)$ and $\mathcal{GL}_{\mbox{\tiny II}}\left(\nu\right)$
reduce to a logistic distribution for $\nu=1$. For $y_{i}=0$, $z_{i}={v}_{i}$,
whereas for $y_{i}=1$, $z_{i}={w}_{i}$, and the choice equation reduces to
$y_{i}=I\\{z_{i}>0\\}$.

###### Theorem 3

Sampling the utilities in the binomial RUM Given $y_{i}=k$ and holding all
model parameters fixed, the latent variables ${w}_{i}|k>0$ and
${v}_{i}|k<N_{i}$ are conditionally independent. The distributions of
${w}_{i}|k>0$ and ${v}_{i}|k<N_{i}$ are equal in distribution to

$\displaystyle{w}_{i}=\log\left((1+\lambda_{i})\frac{1}{W_{i}^{1/y_{i}}}-\lambda_{i}\right),\quad
y_{i}>0,$ (52)
$\displaystyle{v}_{i}=-\log\left(\frac{1+\lambda_{i}}{\lambda_{i}}\frac{1}{V_{i}^{1/(N_{i}-y_{i})}}-\frac{1}{\lambda_{i}}\right),\quad
y_{i}<N_{i},$ (53)

where $W_{i}$ and $V_{i}$ are iid uniform random numbers.

###  4.2 Ultimate Pólya-Gamma samplers for binomial data

The following steps are the main building stones for the ultimate Pólya-Gamma
sampler for binomial data. We use a Gaussian mixture representation of the
generalized logistic distribution based on the Pólya-Gamma distribution and
apply iMDA to achieve robustness against imbalanced data.

A random variable $\varepsilon$ following the generalized logistic
distributions of type I or type II can be represented as mixture of normals
using the Polya-Gamma distribution as mixing measure,

$f_{\varepsilon}(\varepsilon)=c(a,b)\frac{(e^{\varepsilon})^{a}}{(1+e^{\varepsilon})^{b}}=\frac{c(a,b)}{2^{b}}\exp(\frac{\omega\xi^{2}}{2})\int_{0}^{\infty}\exp(-\frac{\omega(\varepsilon-\xi)^{2}}{2})p(\omega)d\omega$
(54)

where $\xi=\kappa/\omega$ and $\kappa=a-b/2$, and
$\omega\sim\mathcal{PG}\left(b,0\right)$ follows the Polya-Gamma distribution
introduced by Polson et al. (2013), see Appendix A.1.1 to A.1.3.

For $y_{i}>0$, the type II generalized logistic distribution
$\varepsilon_{{w},i}\sim\mathcal{GL}_{\mbox{\tiny II}}\left(y_{i}\right)$ in
(51) has such a representation with, see (A.11):

$\displaystyle\kappa_{{w},i}=\frac{1-y_{i}}{2},\qquad\xi_{{w},i}=\frac{\kappa_{{w},i}}{\omega_{{w},i}},\qquad\omega_{{w},i}\sim\mathcal{PG}\left(y_{i}+1,0\right).$

Similarly, for $y_{i}<N_{i}$, the type I generalized logistic distribution
$\varepsilon_{{v},i}\sim\mathcal{GL}_{\mbox{\tiny I}}\left(N_{i}-y_{i}\right)$
in (51) has such a representation with, see (A.7):

$\displaystyle\kappa_{{v},i}=\frac{N_{i}-y_{i}-1}{2},\qquad\xi_{{v},i}=\frac{\kappa_{{v},i}}{\omega_{{v},i}},\qquad\omega_{{v},i}\sim\mathcal{PG}\left(N_{i}-y_{i}+1,0\right).$

Note that $\kappa_{{w},i}=0$ for $y_{i}=1$ and $\kappa_{{v},i}=0$ for
$y_{i}=N_{i}-1$. Hence, for $N_{i}=1$, the Pólya-Gamma mixture approximation
of a logistic model involving $\mathcal{PG}\left(2,0\right)$ results. For
$N_{i}>1$, $\kappa_{{v},i}>0$ for $0\leq y_{i}\leq N_{i}-2$ and
$\kappa_{{w},i}<0$ for $2\leq y_{i}\leq N_{i}$.

Algorithm 3 The ultimate Pólya-Gamma sampler for binomial data.

Choose starting values for
$\boldsymbol{\lambda}=(\lambda_{1},\dots\lambda_{N})$ and repeat the following
steps:

  * (a1)

For each $i=1,\ldots,N$, sample ${w}_{i}|\lambda_{i},(y_{i}>0)$ and
${v}_{i}|\lambda_{i},(y_{i}<N_{i})$ using Theorem 3 and sample
$\omega_{{w},i}|{w}_{i},\lambda_{i},(y_{i}>0)$ and
$\omega_{{v},i}|{v}_{i},\lambda_{i},(y_{i}<N_{i})$ using (4.2).

  * (a2)

Sample $\delta^{\star}\sim\mathcal{G}^{-1}\left(d_{0},D_{0}\right)$ and
$\gamma^{\star}\sim\gamma|\delta^{\star},\boldsymbol{\lambda}$. For a binomial
regression model, $\gamma|\delta^{\star},\beta_{d}$ is given in (A.3.3).
Define $\tilde{w}_{i}=\sqrt{\delta^{\star}}{w}_{i}+\gamma^{\star}$, if
$y_{i}>0$, and $\tilde{v}_{i}=\sqrt{\delta^{\star}}{v}_{i}+\gamma^{\star}$, if
$y_{i}<N_{i}$, for $i=1,\ldots,N$, to move to the unidentified model.

  * (b-1)

Sample $\gamma^{\mathrm{new}}$ from the truncated Student-$t$ prior using
(62). Sample $\delta^{\mathrm{new}}\sim
p(\delta|\gamma^{\mathrm{new}},\tilde{\boldsymbol{z}},{\mathbf{y}})\propto\left(\frac{1}{\delta}\right)^{d_{I}+1}e^{-D_{I}/\delta}e^{B_{I}/\sqrt{\delta}}$
using the resampling described in Appendix A.3.3. For a binomial regression
model, $d_{I},D_{I}$ and $B_{I}$ are given by (A.36).

  * (b-2)

Sample $\tilde{\boldsymbol{\beta}}\sim
p(\tilde{\boldsymbol{\beta}}|\tilde{\boldsymbol{z}},\delta^{\mathrm{new}})$.
For a binomial regression model,
$\tilde{\boldsymbol{\beta}}|\tilde{\boldsymbol{z}},\delta^{\mathrm{new}}\sim\mathcal{N}_{d}\left({\mathbf{b}}_{N}(\delta^{\mathrm{new}}),\delta^{\mathrm{new}}{\mathbf{B}}_{N}\right)$
where ${\mathbf{b}}_{N}(\delta)$ and ${\mathbf{B}}_{N}$ are given by (A.35).
Move back to the identified model using
$\log\lambda_{i}=(\log\tilde{\lambda}_{i}-\gamma^{\mathrm{new}})/\sqrt{\delta^{\mathrm{new}}}$
and discard $\tilde{\boldsymbol{z}}$. For a binomial regression model, e.g.,
use transformation (11).

For each $i=1,\ldots,N$, we introduce the latent variables
$\boldsymbol{z}_{i}=({w}_{i},\omega_{{w},i},{v}_{i},\omega_{{v},i})$, if
$0<y_{i}<N_{i}$, $\boldsymbol{z}_{i}=({w}_{i},\omega_{{w},i})$, if
$y_{i}=N_{i}$, and $\boldsymbol{z}_{i}=({v}_{i},\omega_{{v},i})$, if
$y_{i}=0$. Conditional on holding all unknown parameters in $\log\lambda_{i}$
fixed, the latent variables ${w}_{i}|\lambda_{i},(y_{i}>0)$ and
${v}_{i}|\lambda_{i},(y_{i}<N_{i})$ are sampled from Theorem 3 without
conditioning on $\omega_{{w},i}$ and $\omega_{{v},i}$. Given the latent
variables ${w}_{i}$ and ${v}_{i}$, the parameters
$\omega_{{w},i}|{w}_{i},(y_{i}>0),\lambda_{i}$ and
$\omega_{{v},i}|{v}_{i},(y_{i}<N_{i}),\lambda_{i}$ are independent and follow
(tilted) Pólya-Gamma distributions:

$\displaystyle\omega_{{w},i}|{w}_{i},y_{i},\lambda_{i}\sim\mathcal{PG}\left(y_{i}+1,|{w}_{i}-\log\lambda_{i}|\right),\quad
y_{i}>0,$
$\displaystyle\omega_{{v},i}|{v}_{i},y_{i},\lambda_{i}\sim\mathcal{PG}\left(N_{i}-y_{i}+1,|{v}_{i}-\log\lambda_{i}|\right),\quad
y_{i}<N_{i}.$

To protect our sampler against imbalanced data, we apply iMDA and transform
the latent variables $\tilde{w}_{i}$ and $\tilde{v}_{i}$ as for a logit model:

$\displaystyle\tilde{w}_{i}=\sqrt{\delta}{w}_{i}+\gamma,\quad\tilde{v}_{i}=\sqrt{\delta}{v}_{i}+\gamma,$
(56)

which implies $\log\tilde{\lambda}_{i}=\sqrt{\delta}\log\lambda_{i}+\gamma$.
Transformation (56) leads to an unidentified, boosted version of model (51)
with a choice equation involving $\gamma$:

$\displaystyle\tilde{w}_{i}=\log\tilde{\lambda}_{i}+\sqrt{\delta}\varepsilon_{{w},i},\quad
y_{i}>0,$ (57)
$\displaystyle\tilde{v}_{i}=\log\tilde{\lambda}_{i}+\sqrt{\delta}\varepsilon_{{v},i},\quad
y_{i}<N_{i},$ $\displaystyle
y_{i}=k\Leftrightarrow\left\\{\begin{array}[]{ll}\tilde{v}_{i}<\gamma<\tilde{w}_{i},&0<k<N_{i},\\\
\gamma>\tilde{v}_{i},&k=0,\\\ \gamma<\tilde{w}_{i},&k=N_{i}.\\\
\end{array}\right.$ (61)

The parameters $\delta$ and $\gamma$ have the same prior as in the logit
model, see (7). The choice equation (61) implies following constraint for
$\gamma$, conditional on
$\tilde{\boldsymbol{z}}=(\tilde{\boldsymbol{z}}_{1},\ldots,\tilde{\boldsymbol{z}}_{N})$,
where
$\tilde{\boldsymbol{z}}_{i}=(\tilde{w}_{i},\omega_{{w},i},\tilde{v}_{i},\omega_{{v},i})$:

$\displaystyle
L(\tilde{\boldsymbol{z}})<\gamma<O(\tilde{\boldsymbol{z}}),\quad
L(\tilde{\boldsymbol{z}})=\max(\tilde{v}_{i}),\quad
O(\tilde{\boldsymbol{z}})=\min(\tilde{w}_{i}).$ (62)

However, due to the presence of the (fixed) location parameters $\xi_{{w},i}$
and $\xi_{{v},i}$ in the likelihood (54), as opposed to a binary or a MNL
model, the posterior $p(\delta|\gamma,\tilde{\boldsymbol{z}})$ does not take a
closed form. Only if $N_{i}=1$ for all $i$, then all $\xi_{{w},i}$s and
$\xi_{{v},i}$s are zero and an inverted gamma posterior for $\delta$ results.
If $\log\lambda_{i}$ is a linear function of the unknown model parameters
$\boldsymbol{\beta}$, then completing squares yields a conditionally Gaussian
distribution for
$p(\tilde{\boldsymbol{\beta}}|\gamma,\delta,\tilde{\boldsymbol{z}})$ and a
closed form expression for $p(\delta|\gamma,\tilde{\boldsymbol{z}})$ is
available which, however, is not longer an inverted Gamma posterior as for the
logit model, but a more general distribution taking the form

$p(\delta|\gamma,\tilde{\boldsymbol{z}})\propto\left(\frac{1}{\delta}\right)^{d_{I}+1}e^{-D_{I}/\delta}e^{B_{I}/\sqrt{\delta}}.$

Details how to sample from such a distribution are provided in Appendix A.3.3.
The ultimate Pólya-Gamma sampler for binomial data is summarized in Algorithm
3.

##  5 Illustrative applications to simulated data

In this section, we discuss a set of systematic simulations carried out to
compare the efficiency of the UPG approach to other popular Bayesian samplers
for discrete choice models. The main results are based on simulations with
varying intercepts in the true data generating process. Large negative and a
large positive intercepts correspond to simulation setups resulting in highly
imbalanced data.444In the main simulation results, we fix $N=1000$ for probit,
logit and multinomial logit models. For the multinomial case, we consider
three categories. Finally, for binomial models, we fix $N=500$ and $N_{i}=10$.
For logit and multinomial logit specifications, we also look at a varying
sample size when the number of observed successes is kept constant at one. For
the multinomial case, we furthermore investigate the scenario of a growing
number of categories/choices. Finally, the binomial logit case allows us to
vary the number of trials $N_{i}$ observed in the data. These additional
results as well as results for probit models and a number of numerical results
can be found in the Appendix A.5.

Each Markov chain was run for 10,000 iterations after a burn in phase of 2,000
iterations. To gain some robustness with respect to the computed inefficiency
factors, each simulation is repeated 50 times. Below, the median results from
these 50 runs are discussed. To compare the efficiency of different samplers,
we provide inefficiency factors based on an estimate of the spectral density
of the posterior chain evaluated at zero.555Estimating the spectral density at
0 is accomplished via the R package coda (Plummer et al., 2006) and is based
on fitting an autoregressive process to the posterior draws.

For binary logit models, we compare the UPG sampling scheme, the Pólya-Gamma
sampler of Polson et al. (2013) (P&S), the adaptive MH step proposed in
Johndrow et al. (2019) (AMH) and the auxiliary mixture data augmentation
scheme outlined in Frühwirth-Schnatter & Frühwirth (2010) (FS&F). To assess
sampling efficiency for the multinomial logit model (MNL), we compare the MNL
implementation from the technical appendix of Polson et al. (2013) (P&S), a
version of the AMH, the partial dRUM sampler from Frühwirth-Schnatter &
Frühwirth (2010) (FS&F) as well as the UPG approach outlined in Section 3.
Finally, to illustrate the efficiency gains in the case of regression analysis
of binomial data, we compare the UPG approach outlined in Section 4 to a
sampler based on an AMH algorithm and the auxiliary mixture sampler for
binomial logit models introduced in Fussl et al. (2013) (AuxMix).

(a) Binary logit.

(b) Multinomial logit.

(c) Binomial logit.

Figure 3: Inefficiency factors across 50 simulation runs for varying
intercepts $\beta_{d}$.

The results of the main simulation exercise are summarized in Figure 3. The
results show that standard data augmentation techniques exhibit extremely
inefficient sampling behavior when increasing data sparsity as pointed out in
Johndrow et al. (2019). While this general pattern can be observed for the UPG
approach as well, our novel boosting strategy clearly alleviates this issue
and allows for rather efficient estimation also in imbalanced data
settings.666Interestingly, all samplers exhibit asymmetric inefficiency
factors with respect to baseline sparsity in the MNL case. That is, one has to
differentiate between cases with an increasingly sparse baseline (positive
intercept range) and increasingly sparse non-baseline groups (negative
intercept range). A sparse baseline seems to increase inefficiency in all
observed cases. This is rather intuitive given that intercepts are the log of
the probability of a specific category divided by the log of the probability
of the baseline/reference. Generally the reference category should be a
”common” rather than a ”rare” category. However, this problem can be avoided
in applied statistical modelling as the baseline can be chosen by the
researcher, effectively eliminating this issue.

A direct comparison between the UPG approach and the original Pólya-Gamma
sampler shows that the former clearly outperforms the latter in terms of
efficiency in settings with imbalanced data. That is, UPG samplers are the
preferable data augmentation sampler for imbalanced data cases. However, in
data sets that are very balanced, the original Pólya-Gamma sampler exhibits
the highest observed efficiency levels. This is most likely due to the
additional layer of latent variables that the UPG approach utilizes compared
to P&S. Nevertheless, note that all data augmentation algorithms under
consideration even surpass the performance of the efficiency ”gold standard”
AMH in balanced scenarios. In other words, these settings correspond to data
sets where inefficient estimation plays a minor role to begin with.

In summary, UPG as well as Metropolis-Hastings based estimation allow for
efficient sampling across a wide range of settings when compared to standard
data augmentation algorithms. This can also be seen through extensive,
additional simulations that are summarized in subsection A.5. It is worth to
mention that the UPG approach has some benefits compared to Metropolis-
Hastings based samplers. First, no tuning is necessary. Second, one of the
major advantages of the UPG approach is that, for all models discussed in this
paper, the conditional posterior of the regression effects is multivariate
normal. This enables the researcher to readily rely on the large statistical
toolbox of modifications and extensions of the linear regression model that
has been developed over the past decades. We will proceed to demonstrate this
convenient fact through applying the UPG approach to one of the most general
modeling frameworks in Section 6, where we discuss binary state space modeling
and present a brief real data application.

##  6 Application to a binary state space model

Let $\left\\{y_{t}\right\\}$ be a time series of binary observations, observed
for $t=1,\ldots,T$. Each $y_{t}$ is assumed to take one of two possible
values, labelled by $\\{0,1\\}$. The probability that $y_{t}$ takes the value
$1$ depends on covariates ${\mathbf{x}}_{t}$, including a constant (i.e.
$x_{td}\equiv 1$), through time-varying parameters $\boldsymbol{\beta}_{t}$ in
the following way:

$\displaystyle\mbox{\rm
Pr}(y_{t}=1|\boldsymbol{\beta}_{1},\ldots,\boldsymbol{\beta}_{T})=\frac{\exp({\mathbf{x}}_{t}\boldsymbol{\beta}_{t})}{1+\exp({\mathbf{x}}_{t}\boldsymbol{\beta}_{t})}.$
(63)

We assume that conditional on knowing
$\boldsymbol{\beta}_{1},\ldots,\boldsymbol{\beta}_{T}$, the observations are
mutually independent. A commonly used model for describing the time-variation
of $\boldsymbol{\beta}_{t}$ reads:

$\displaystyle\boldsymbol{\beta}_{t}=\boldsymbol{\beta}_{t-1}+{\mathbf{w}}_{t},\quad{\mathbf{w}}_{t}\sim\mathcal{N}_{d}\left(\bf
0,{\mathbf{{\mathbf{Q}}}}\right),$ (64)

where ${\mathbf{{\mathbf{Q}}}}=\mbox{\rm
Diag}\left(\theta_{1},\ldots,\theta_{d}\right)$. The expected initial values
$\boldsymbol{\beta}=(\beta_{1},\ldots,\beta_{d})^{\top}$ and the variances
$\theta_{1},\ldots,\theta_{d}$ are unknown. To complete the model
specification, the distribution of the initial value $\boldsymbol{\beta}_{0}$
has to be specified. We will discuss this choice in connection with our
boosting strategy.

Markov chain Monte Carlo estimation of such logit-type state space models has
been considered by various authors, in particular by Shephard & Pitt (1997)
and Gamerman (1998). A characteristic feature of several MCMC approaches for
binary state space models is that practical implementation requires the use of
an MH algorithm for sampling the state process, which in turn makes it
necessary to define a suitable proposal density in a rather high-dimensional
parameter space. Single-move sampling for this type of models is known to be
potentially very inefficient; see e.g. Shephard & Pitt (1997). Frühwirth-
Schnatter & Frühwirth (2007) implemented an auxiliary mixture sampler for a
binary regression model with time-varying parameters.

We show in this section that Algorithm 1 can be quite easily extended to
binary state space models, as Step (a) is generic conditional on knowing
$\log\lambda_{t}={\mathbf{x}}_{t}\boldsymbol{\beta}_{t}$.

###  6.1 The UPG sampler for binary state space models

A latent utility $z_{t}$ of choosing category 1 is introduced for each
$y_{t}$:

$\displaystyle\boldsymbol{\beta}_{t}=\boldsymbol{\beta}_{t-1}+{\mathbf{w}}_{t},\quad{\mathbf{w}}_{t}\sim\mathcal{N}_{d}\left(\bf
0,{\mathbf{{\mathbf{Q}}}}\right),$ (65) $\displaystyle y_{t}=1\Leftrightarrow
z_{t}>0,\qquad z_{t}={\mathbf{x}}_{t}\boldsymbol{\beta}_{t}+\varepsilon_{t}.$

For a probit link, (65) is a conditionally Gaussian state space model (SSM),
whereas for a logit link, $\varepsilon_{t}$ follows a logistic distribution
and (65) is a conditionally non-Gaussian. However, using the Pólya Gamma
mixture representation of the logistic distribution as in Section 2 yields a
(high-dimensional) Gaussian posterior for the state process
$\boldsymbol{\beta}=\\{\boldsymbol{\beta}_{0},\boldsymbol{\beta}_{1},\ldots,\boldsymbol{\beta}_{T}\\}$,
which allows multi-move sampling of the entire state process using FFBS
(Frühwirth-Schnatter, 1994; Carter & Kohn, 1994) in a similar fashion as for a
probit SSM.

To achieve robustness again imbalance, we apply iMDA as introduced in Section
2:

$\displaystyle\tilde{z}_{t}=\sqrt{\delta}z_{t}+\gamma,$ (66)

to transform the SSM (65) to the following unidentified SSM

$\displaystyle\tilde{\beta}_{jt}=\tilde{\beta}_{j,t-1}+w_{jt},\quad
w_{jt}\sim\mathcal{N}\left(0,\delta\theta_{j}\right),\quad j=1,\ldots,d,$
$\displaystyle
y_{t}=1\Leftrightarrow\tilde{z}_{t}>\gamma,\quad\tilde{z}_{t}={\mathbf{x}}_{t}\tilde{\boldsymbol{\beta}}_{t}+\sqrt{\delta}\varepsilon_{t},\quad$
(67)

where transformation (66) implies

$\displaystyle\tilde{\beta}_{j}=\sqrt{\delta}\beta_{j}+I\\{j=d\\}\gamma,\quad\tilde{\beta}_{jt}=\sqrt{\delta}\beta_{{jt}}+I\\{j=d\\}\gamma,\quad
t=0,\ldots,T,\quad j=1,\ldots,d.$ (68)

The ultimate Pólya-Gamma sampler for binary SSMs is summarized in Appendix
A.4, Algorithm 4.

###  6.2 Illustration using pandemic outbreak data

Figure 4: Posterior of local level model fitted to global pandemic data. Solid
line is posterior median, dashed lines are 90/10 posterior quantiles. Figure
5: Inefficiency factors for boosted and unboosted logit local level models
fitted to global pandemic data.

To illustrate the gains in sampling efficiency for binary state space models,
we apply the framework outlined in the previous subsections to a dataset on
global pandemics. The data covers 220 years from 1800 to 2020. We define a
pandemic as a disease episode that can be characterized by worldwide spread of
the disease and a relatively high death toll ($>$ 75,000 deaths). In addition,
we focus on diseases that are characterized by relatively short periods of
activity, hence excluding global pandemics such as AIDS.

This leaves us with a total of eight pandemic events, starting with a bubonic
plague outbreak between 1855 and 1860. With more than 12,000,000 deaths in
China and India alone, this is considered the third major outbreak of the
plague in human history, following prior pandemics in the 6th and 14th
century. As second event in the 19th century, we consider a major influenza
outbreak from 1889-1890, also known as ”Russian flu” or ”Asiatic flu”. The
beginning of the 20th century was hit especially hard by pandemics, starting
with an outbreak of Encephalitis lethargica from 1915-1926. Within the same
time period, the well known Spanish flu pandemic took place between 1918 and
1920. Being one of the deadliest recorded pandemics, the Spanish flu
ultimately eradicated around 25% of the world’s population. After the second
world war, we identify four major global pandemic events: Still in the 20th
century, two influenza A strains known as ”Asian flu” (1957-1958) and ”Hong
Kong flu” (1968 - 1969) killed around 3,000,000 people. In more recent
history, a H1N1 virus (”swine flu”) spread globally in 2009 and 2010, leading
to up to 575,000 deaths. Finally, the most recent event is the ongoing global
outbreak of COVID-19, starting in 2019 and claiming around 1,294,000 victims
as of November, 13th 2020.777Sourced from
https://en.wikipedia.org/wiki/List_of_epidemics and sources therein.

For years featuring a global pandemic, the binary indicator $y_{t}$ equals 1.
In years without major pandemic events, $y_{t}=0$. A pandemic is observed in
roughly 1 out of 8 years with extremely high autocorrelation, rendering the
data set relatively imbalanced. To illustrate the efficiency gains arising
from our boosting algorithm, a logistic local level model is fitted to the
data twice, once with and once without boosting, enabling us to directly
compare these samplers in terms of efficiency. We iterate the Gibbs sampler
500,000 times after a burn-in period of 20,000 iterations. The resulting
posterior distribution is shown in Figure 4. The resulting inefficiency
factors are provided in Figure 5. Clearly, long periods of persistence in the
binary outcome time series lead to an aggravated efficiency loss in the
sampler. On the contrary, during periods where $y_{t}$ is likely to change
(corresponding to a predicted probability around 0.5), sampling efficiency is
rather high. Summarizing, the boosted sampling algorithm described in Section
6.1 is able to drastically reduce sampling inefficiency, even during prolonged
periods of no change in the outcome variable. Hence, the UPG approach can be
considered a useful tool for efficient discrete choice analysis, not only in
simple regression models, but also in more complex settings such as state
space frameworks.

##  7 Concluding Remarks

Due to a wide range of applications in many areas of applied science, much
attention has been dedicated towards the development of estimation algorithms
for generalized linear models. In the past decades, various data augmentation
algorithms have been brought forward that have steadily increased
accessibility and popularity of Bayesian estimation techniques in the context
of discrete choice models.

In this article, we introduce a number of sampling algorithms based on the
data augmentation algorithm of Polson et al. (2013). These algorithms may be
used to estimate parameters in the most commonly encountered discrete choice
models and are suited to model binary, categorical and binomial outcomes. The
samplers are easily implemented, intuitively appealing and allow for a
conditionally Gaussian posterior of the regression effects in logit,
multinomial logit and binomial logit regression type models.

While Gibbs sampling algorithms are well-known and require only minimal
tuning, it has been recently noted that estimation based on Metropolis-
Hastings algorithms exhibit higher estimation efficiency, especially in
applications with imbalanced data (Johndrow et al., 2019). To counteract
efficiency loss in such settings, we develop a novel MCMC boosting strategy
based on marginal data augmentation and apply it to the introduced sampling
algorithms. The proposed UPG algorithms combine the appeal of Gibbs samplers
with a comparative level of sampling efficiency, even in unbalanced cases.
This is demonstrated via extensive simulation studies and a real data
application.

## References

  * (1)
  * Albert & Chib (1993) Albert, J. H. & Chib, S. (1993), ‘Bayesian analysis of binary and polychotomous response data’, Journal of the American Statistical Association 88, 669–679.
  * Carter & Kohn (1994) Carter, C. K. & Kohn, R. (1994), ‘On Gibbs sampling for state space models’, Biometrika 81, 541–553.
  * Frühwirth-Schnatter (1994) Frühwirth-Schnatter, S. (1994), ‘Data augmentation and dynamic linear models’, Journal of Time Series Analysis 15, 183–202.
  * Frühwirth-Schnatter & Frühwirth (2007) Frühwirth-Schnatter, S. & Frühwirth, R. (2007), ‘Auxiliary mixture sampling with applications to logistic models’, Computational Statistics and Data Analysis 51, 3509–3528.
  * Frühwirth-Schnatter & Frühwirth (2010) Frühwirth-Schnatter, S. & Frühwirth, R. (2010), Data augmentation and MCMC for binary and multinomial logit models, in T. Kneib & G. Tutz, eds, ‘Statistical Modelling and Regression Structures – Festschrift in Honour of Ludwig Fahrmeir’, Physica-Verlag, Heidelberg, pp. 111–132.
  * Frühwirth-Schnatter et al. (2009) Frühwirth-Schnatter, S., Frühwirth, R., Held, L. & Rue, H. (2009), ‘Improved auxiliary mixture sampling for hierarchical models of non-Gaussian data’, Statistics and Computing 19, 479–492.
  * Fussl et al. (2013) Fussl, A., Frühwirth-Schnatter, S. & Frühwirth, R. (2013), ‘Efficient MCMC for binomial logit models’, ACM Transactions on Modeling and Computer Simulation 22, 3:1–3:21.  doi: 10.1145/2414416.2414419.
  * Gamerman (1998) Gamerman, D. (1998), ‘Markov chain Monte Carlo for dynamic generalized linear models’, Biometrika 85, 215–227.
  * Holmes & Held (2006) Holmes, C. C. & Held, L. (2006), ‘Bayesian auxiliary variable models for binary and multinomial regression’, Bayesian Analysis 1, 145–168.
  * Imai & van Dyk (2005) Imai, K. & van Dyk, D. A. (2005), ‘A Bayesian analysis of the multinomial probit model using marginal data augmentation’, Journal of Econometrics 124, 311–334.
  * Johndrow et al. (2019) Johndrow, J. E., Smith, A., Pillai, N. & Dunson, D. B. (2019), ‘MCMC for imbalanced categorical data’, Journal of the American Statistical Association 114, 1394–1403.
  * McCulloch et al. (2000) McCulloch, R. E., Polson, N. G. & Rossi, P. E. (2000), ‘A Bayesian analysis of the multinomial probit model with fully identified parameters’, Journal of Econometrics 99, 173–193.
  * McFadden (1974) McFadden, D. (1974), Conditional logit analysis of qualitative choice behaviour, in P. Zarembka, ed., ‘Frontiers of Econometrics’, Academic, New York, pp. 105–142.
  * Meng & Van Dyk (1999) Meng, X.-L. & Van Dyk, D. (1999), ‘Seeking efficient data augmentation schemes via conditional and marginal data augmentation’, Biometrika 86, 301–320.
  * Plummer et al. (2006) Plummer, M., Best, N., Cowles, K. & Vines, K. (2006), ‘CODA: Convergence diagnosis and output analysis for MCMC’, R News 6(1), 7–11.
  * Polson et al. (2013) Polson, N. G., Scott, J. G. & Windle, J. (2013), ‘Bayesian inference for logistic models using Pólya-Gamma latent variables’, Journal of the American Statistical Association 108, 1339–49.
  * Robert & Casella (1999) Robert, C. P. & Casella, G. (1999), Monte Carlo Statistical Methods, Springer Series in Statistics, Springer, New York/Berlin/Heidelberg.
  * Rossi et al. (2005) Rossi, P. E., Allenby, G. M. & McCulloch, R. (2005), Bayesian Statistics and Marketing, Wiley, Chichester.
  * Shephard & Pitt (1997) Shephard, N. & Pitt, M. K. (1997), ‘Likelihood analysis of non-Gaussian measurement time series’, Biometrika 84, 653–667.
  * van Dyk & Meng (2001) van Dyk, D. & Meng, X.-L. (2001), ‘The art of data augmentation’, Journal of Computational and Graphical Statistics 10, 1–50.
  * Zellner & Rossi (1984) Zellner, A. & Rossi, P. E. (1984), ‘Bayesian analysis of dichotomous quantal response models’, Journal of Econometrics 25, 365–393.

##  Appendix A Appendix

###  A.1 Mixture representations

####  A.1.1 The Pólya-Gamma mixture representation

For all latent variable representations derived in this paper for binary,
multinomial or categorical data, the error term in the latent equations arises
from a distribution $\varepsilon\sim F_{\varepsilon}(\varepsilon)$, for which
the density $f_{\varepsilon}(\varepsilon)$ can be represented as a mixture of
normals using the Pólya-Gamma distribution as mixing measure:

$\displaystyle
f_{\varepsilon}(\varepsilon)=c(a,b)\frac{(e^{\varepsilon})^{a}}{(1+e^{\varepsilon})^{b}}=\frac{c(a,b)}{2^{b}}{e^{\kappa\varepsilon}}\int
e^{-\frac{\omega\varepsilon^{2}}{2}}p(\omega)d\omega,$ (A.1)

where $\kappa=a-b/2$ and $\omega\sim\mathcal{PG}\left(b,0\right)$ follows the
Pólya-Gamma distribution introduced by Polson et al. (2013) with parameter
$b$.

Conditional on $\omega$ the likelihood contribution of $\varepsilon$ is
proportional to that of a $\mathcal{N}\left(\kappa/\omega,1/\omega\right)$
observation. 888Based on rewriting
$\displaystyle\kappa\varepsilon-\frac{\omega\varepsilon^{2}}{2}=\frac{-\omega(\varepsilon-\kappa/\omega)^{2}}{2}+d,$
(A.2) where $d$ is a constant not depending on $\varepsilon$. This new
representation is very convenient, as the conditional posterior of
$\omega|\varepsilon$ can be derived from following (tilted) Pólya-Gamma
distribution with the same parameter $b$:

$\displaystyle\quad\omega|\varepsilon\sim\mathcal{PG}\left(b,|\varepsilon|\right).$
(A.3)

To simulate from a (tilted) Pólya-Gamma distribution
$\mathcal{PG}\left(q,c\right)$ distribution, the following convolution
property is exploited:

$\displaystyle X_{1}\sim\mathcal{PG}\left(q_{1},c\right),\quad
X_{2}\sim\mathcal{PG}\left(q_{2},c\right)\Rightarrow
X_{1}+X_{2}\sim\mathcal{PG}\left(q_{1}+q_{2},c\right).$

Hence, to simulate from $Y\sim\mathcal{PG}\left(q,c\right)$, use
$Y=\sum_{j=1}^{q}X_{j}$ , where $X_{j}\sim\mathcal{PG}\left(1,c\right)$ are
$q$ independent r.v. from the $\mathcal{PG}\left(1,c\right)$ distribution.

####  A.1.2 The logistic and the type I generalized logistic distribution

For the type I generalized logistic distribution
$\varepsilon\sim\mathcal{GL}_{\mbox{\tiny I}}\left(\nu\right)$ with parameter
$\nu>0$, the density reads

$\displaystyle f_{\varepsilon}(\varepsilon)=\frac{\nu
e^{-\varepsilon}}{(1+e^{-\varepsilon})^{\nu+1}}=\frac{\nu
e^{\nu\varepsilon}}{(1+e^{\varepsilon})^{\nu+1}}.$ (A.4)

$\mathcal{GL}_{\mbox{\tiny I}}\left(\nu\right)$ reduces to the logistic
distribution for $\nu=1$. The c.d.f. of a type I generalized logistic
distribution takes a simple form:

$\displaystyle
F_{\varepsilon}(\varepsilon)=\frac{1}{(1+e^{-\varepsilon})^{\nu}}=\frac{e^{\nu\varepsilon}}{(1+e^{\varepsilon})^{\nu}}.$
(A.5)

Hence, the quantiles are available in closed form:

$\displaystyle\varepsilon_{p}=F_{\varepsilon}^{-1}(p)=-\log\left(\frac{1}{p^{1/\nu}}-1\right).$
(A.6)

The type I generalized logistic distribution
$\varepsilon\sim\mathcal{GL}_{\mbox{\tiny I}}\left(\nu\right)$ can be
represented as a mixture of normals with a Pólya-Gamma distribution serving as
mixing measure, where

$\displaystyle\omega\sim\mathcal{PG}\left(\nu+1,0\right),\quad\kappa=\frac{\nu-1}{2},$
(A.7)

see Section A.1.1. For the logistic distribution,
$\omega\sim\mathcal{PG}\left(2,0\right)$ and $\kappa=0$.

####  A.1.3 The type II generalized logistic distribution

For the type II generalized logistic distribution
$\varepsilon\sim\mathcal{GL}_{\mbox{\tiny II}}\left(\nu\right)$ with parameter
$\nu>0$, the density reads

$\displaystyle f_{\varepsilon}(\varepsilon)=\frac{\nu
e^{-\nu\varepsilon}}{(1+e^{-\varepsilon})^{\nu+1}}=\frac{\nu
e^{\varepsilon}}{(1+e^{\varepsilon})^{\nu+1}}.$ (A.8)

Also $\mathcal{GL}_{\mbox{\tiny II}}\left(\nu\right)$ reduces to the logistic
distribution for $\nu=1$.

The c.d.f. of a type II generalized logistic distribution takes a simple form:

$\displaystyle
F_{\varepsilon}(\varepsilon)=1-\frac{1}{(1+e^{\varepsilon})^{\nu}}=1-\frac{e^{-\nu\varepsilon}}{(1+e^{-\varepsilon})^{\nu}}.$
(A.9)

Hence, the quantiles are available in closed form:

$\displaystyle\varepsilon_{p}=F_{\varepsilon}^{-1}(p)=\log\left(\frac{1}{(1-p)^{1/\nu}}-1\right).$
(A.10)

The type II generalized logistic distribution
$\varepsilon\sim\mathcal{GL}_{\mbox{\tiny II}}\left(\nu\right)$ can be
represented as a mixture of normals with a Pólya-Gamma distribution serving as
mixing measure, where

$\displaystyle\omega\sim\mathcal{PG}\left(\nu+1,0\right),\quad\kappa=\frac{1-\nu}{2},$
(A.11)

see Section A.1.1. Again, for the logistic distribution,
$\omega\sim\mathcal{PG}\left(2,0\right)$ and $\kappa=0$ results.

###  A.2 Proofs

#### Proof of Theorem 1

The proof of Theorem 1 is straightforward. Depending on the observed category
$y_{i}$, the corresponding utility is the maximum among all latent utilities
and the posterior distribution
$p({u}_{ki},{u}_{0i},{u}_{a,i}|\lambda_{ki},y_{i})$ is equal to one of the
following distributions:

$\displaystyle e^{-{u}_{ki}}|y_{i}=k$ $\displaystyle\sim$
$\displaystyle\mathcal{E}\left(\lambda^{\star}_{i}\right)\,\Rightarrow\,e^{-{u}_{0i}}=e^{-{u}_{ki}}+\mathcal{E}\left(1\right),\,e^{-{u}_{a,i}}=e^{-{u}_{ki}}+\mathcal{E}\left(\lambda_{a,i}\right),$
$\displaystyle e^{-{u}_{0i}}|y_{i}=0$ $\displaystyle\sim$
$\displaystyle\mathcal{E}\left(\lambda^{\star}_{i}\right)\,\Rightarrow\,e^{-{u}_{ki}}=e^{-{u}_{0i}}+\mathcal{E}\left(\lambda_{ki}\right),\,e^{-{u}_{a,i}}=e^{-{u}_{0i}}+\mathcal{E}\left(\lambda_{a,i}\right),$
$\displaystyle e^{-{u}_{a,i}}|y_{i}\in A$ $\displaystyle\sim$
$\displaystyle\mathcal{E}\left(\lambda^{\star}_{i}\right)\,\Rightarrow\,e^{-{u}_{ki}}=e^{-{u}_{a,i}}+\mathcal{E}\left(\lambda_{ki}\right),\,e^{-{u}_{0i}}=e^{-{u}_{a,i}}+\mathcal{E}\left(1\right),$

where $\lambda^{\star}_{i}=1+\sum_{l=1}^{m}\lambda_{li}$. For efficient joint
sampling of all utilities for all $i=1,\ldots,N$, this can be rewritten as in
(39).

#### Proof of Theorem 2

A binomial observation $y_{i}$ can be regarded as the aggregated number of
successes among $N_{i}$ independent binary outcomes
$z_{1i},\ldots,z_{N_{i},i}$, labelled $\\{0,1\\}$ and each following the
binary logit model $\mbox{\rm Pr}(z_{ni}=1|\pi_{i})=\pi_{i}$. For each
individual binary observation $z_{ni}$, the logit model can be written as a
RUM (McFadden 1974):

$\displaystyle\displaystyle{u}_{ni}=\log\lambda_{i}+\epsilon_{ni},\qquad\epsilon_{ni}\sim\mathcal{LO},$
$\displaystyle z_{ni}=I\\{{u}_{ni}>0\\},$

involving a latent variable ${u}_{ni}$, where $\epsilon_{ni}$ are i.i.d.
errors following a logistic distribution. Among the $N_{i}$ binary experiment,
$y_{i}$ outcomes $z_{ni}$ choose the category 1, whereas the remaining
$N_{i}-y_{i}$ outcomes $z_{ni}$ choose the category 0. The challenge is to
aggregate the latent variables ${u}_{ni}$ to a few latent variables in such a
way that an explicit choice equation is available. As it turns out, such an
aggregation can be based on the order statistics
${u}_{(1),i}<\ldots<{u}_{(N_{i}),i}$ of ${u}_{1i},\ldots,{u}_{N_{i},i}$.

Consider first the case that $y_{i}=0$. Such an outcome is observed, iff
$z_{ni}=0$ or, equivalently, the latent utility is negative (${u}_{ni}<0$) for
all $n=1,\ldots,N_{i}$. Hence, a necessary and sufficient condition for
$y_{i}=0$ is that the maximum of all utilities is negative, or equivalently,

$\displaystyle y_{i}=0\,\,\Leftrightarrow\,\,{u}_{(N_{i}),i}<0.$ (A.12)

Next, consider the case that $y_{i}=N_{i}$. Such an outcome is observed, iff
$z_{ni}=1$ or, equivalently, the latent utility is positive (${u}_{ni}>0$ )
for all $n=1,\ldots,N_{i}$. Hence, a necessary and sufficient condition for
$y_{i}=N_{i}$ is that the minimum of all utilities is positive, or
equivalently,

$\displaystyle y_{i}=N_{i}\,\,\Leftrightarrow\,\,{u}_{(1),i}>0.$ (A.13)

Also for outcomes $0<y_{i}<N_{i}$, the order statistics ${u}_{(N_{i}-k),i}$
and ${u}_{(N_{i}-k+1),i}$ provide necessary and sufficient conditions:

$\displaystyle
y_{i}=k,0<k<N_{i}\,\,\Leftrightarrow\,\,{u}_{(N_{i}-k),i}<0,{u}_{(N_{i}-k+1),i}>0.$
(A.14)

Note that (A.12) – (A.14) are choice equations involving either a single or
two order statistics. Hence, we introduce the corresponding order statistics
as aggregated latent variables. Given $y_{i}$, we define

$\displaystyle{v}_{i}={u}_{(N_{i}-k),i},\quad y_{i}>0,$
$\displaystyle{w}_{i}={u}_{(N_{i}-k+1),i},\quad y_{i}<N_{i},$ $\displaystyle
y_{i}=k\Leftrightarrow{v}_{i}<0,{w}_{i}>0,$

with obvious modifications for $y_{i}=0$ and $y_{i}=N_{i}$. It remains to
prove that the latent variables can be represented as in the aggregated model
(51):

$\displaystyle{w}_{i}=\log\lambda_{i}+\varepsilon_{{w},i},\quad\varepsilon_{{w},i}\sim\mathcal{GL}_{\mbox{\tiny
II}}\left(y_{i}\right),\quad y_{i}>0$ (A.15)
$\displaystyle{v}_{i}=\log\lambda_{i}+\varepsilon_{{w},i},\quad\varepsilon_{{v},i}\sim\mathcal{GL}_{\mbox{\tiny
I}}\left(N_{i}-y_{i}\right),\quad y_{i}<N_{i}.$

Note that the order statistics can be represented for $j=1,\ldots,N_{i}$ as
${u}_{(j),i}=\log\lambda_{i}+\varepsilon_{(j),i}$, involving the order
statistics $\varepsilon_{(1),i},\ldots,\varepsilon_{(N_{i}),i}$ are of $N_{i}$
iid realisations $\varepsilon_{1,i},\ldots,\varepsilon_{N_{i},i}$ of a
logistic distribution. Their distribution can be derived from the order
statistics $X_{(1),i},\ldots,X_{(N_{i}),i}$ of $N_{i}$ uniform random numbers
$X_{1i},\ldots,X_{N_{i},i}$ using:

$\displaystyle\varepsilon_{(j),i}=F^{-1}(X_{(j),i})=\log\frac{X_{(j),i}}{1-X_{(j),i}}\quad\Leftrightarrow\quad
X_{(j),i}=F(\varepsilon_{(j),i}),$ (A.16)

where $F$ is the cdf of the logistic distribution.

First, for the special cases where $y_{i}=0$ or $y_{i}=N_{i}$, we use that
$X_{(j),i}\sim\mathcal{B}\left(j,N_{i}-j+1\right)$. Using (A.16), we can
derive the density of $\varepsilon_{(N_{i}),i}$:

$\displaystyle
p(\varepsilon_{(N_{i}),i})=N_{i}F(\varepsilon_{(N_{i}),i})^{N_{i}-1}f(\varepsilon_{(N_{i}),i})=\frac{N_{i}\exp(\varepsilon_{(N_{i}),i})^{N_{i}}}{(1+\exp(\varepsilon_{(N_{i}),i}))^{N_{i}+1}},$

which is the density of a $\mathcal{GL}_{\mbox{\tiny I}}\left(N_{i}\right)$
distribution, see (A.4). Hence, for $y_{i}=0$,

$\displaystyle{v}_{i}=\log\lambda_{i}+\varepsilon_{{v},i},\quad\varepsilon_{{v},i}=\varepsilon_{(N_{i}),i}\sim\mathcal{GL}_{\mbox{\tiny
I}}\left(N_{i}\right).$

Using (A.16), we can derive the density of $\varepsilon_{(1),i}$:

$\displaystyle
p(\varepsilon_{(1),i})=N_{i}(1-F(\varepsilon_{(1),i}))^{N_{i}-1}f(\varepsilon_{(1),i})=\frac{N_{i}\exp(\varepsilon_{(1),i})}{(1+\exp(\varepsilon_{(1),i}))^{N_{i}+1}},$

which is the density of a $\mathcal{GL}_{\mbox{\tiny II}}\left(N_{i}\right)$
distribution, see (A.8). Hence, for $y_{i}=N_{i}$:

$\displaystyle{w}_{i}=\log\lambda_{i}+\varepsilon_{{w},i},\quad\varepsilon_{{w},i}=\varepsilon_{(1),i}\sim\mathcal{GL}_{\mbox{\tiny
II}}\left(N_{i}\right).$

Second, for $0<y_{i}<N_{i}$ we need the joint distribution of
$(\varepsilon_{(N_{i}-k),i},\varepsilon_{(N_{i}-k+1),i})^{\top}$, where
$\varepsilon_{(N_{i}-k+1),i}=\varepsilon_{(N_{i}-k),i}+\Delta\varepsilon_{i}$
with $\Delta\varepsilon_{i}>0$. Using that
$(X_{(j),i},X_{(j+1),i}-X_{(j),i},1-X_{(j+1),i})\sim\mathcal{D}\left(j,1,N_{i}-j\right)$
follows a Dirichlet distribution, see e.g. Robert & Casella (1999), we obtain
that
$(X_{(N_{i}-k),i},1-X_{(N_{i}-k+1),i},X_{(N_{i}-k+1),i}-X_{(N_{i}-k),i})\sim\mathcal{D}\left(N_{i}-k,k,1\right)$.
To derive $p(\varepsilon_{(N_{i}-k),i},\varepsilon_{(N_{i}-k+1),i})$, we
consider the transformations

$\displaystyle\varepsilon_{(N_{i}-k),i}=F^{-1}(X_{(N_{i}-k),i}),\quad\varepsilon_{(N_{i}-k+1),i}=F^{-1}(X_{(N_{i}-k+1),i}),$
(A.17)

and their inverse, $X_{(N_{i}-k),i}=F(\varepsilon_{(N_{i}-k),i})$ and
$X_{(N_{i}-k+1),i}=F(\varepsilon_{(N_{i}-k+1),i})$. We determine

$\displaystyle\left|\frac{\partial(X_{(N_{i}-k),i},X_{(N_{i}-k+1),i})}{\partial(\varepsilon_{(N_{i}-k),i}),\varepsilon_{(N_{i}-k+1),i}))}\right|$
$\displaystyle=$
$\displaystyle\left|\begin{array}[]{cc}f(\varepsilon_{(N_{i}-k),i})&0\\\
0&f(\varepsilon_{(N_{i}-k+1),i})\\\
\end{array}\right|=f(\varepsilon_{(N_{i}-k),i})f(\varepsilon_{(N_{i}-k+1),i}),$

where $f$ is the pdf of the logistic distribution. Therefore,

$\displaystyle p(\varepsilon_{(N_{i}-k),i},\varepsilon_{(N_{i}-k+1),i})=$
$\displaystyle{\frac{\Gamma(N_{i}+1)}{\Gamma(k)\Gamma(N_{i}-k)}}F(\varepsilon_{(N_{i}-k),i})^{N_{i}-k-1}(1-F(\varepsilon_{(N_{i}-k+1),i}))^{k-1}\cdot
f(\varepsilon_{(N_{i}-k),i})f(\varepsilon_{(N_{i}-k+1),i})=$
$\displaystyle\frac{\exp(\varepsilon_{(N_{i}-k),i})^{N_{i}-k}}{(1+\exp(\varepsilon_{(N_{i}-k),i}))^{N_{i}-k+1}}\frac{\exp(\varepsilon_{(N_{i}-k+1),i})}{(1+\exp(\varepsilon_{(N_{i}-k+1),i}))^{k+1}}\cdot{\frac{\Gamma(N_{i}+1)}{\Gamma(k)\Gamma(N_{i}-k)}}I\\{\varepsilon_{(N_{i}-k+1),i}>\varepsilon_{(N_{i}-k),i}\\}.$

This density can be expressed as

$\displaystyle p(\varepsilon_{(N_{i}-k),i},\varepsilon_{(N_{i}-k+1),i})$
$\displaystyle=$ $\displaystyle
p(\varepsilon_{(N_{i}-k),i})\,p(\varepsilon_{(N_{i}-k+1),i})\cdot C\,\cdot
I\\{\varepsilon_{(N_{i}-k+1),i}>\varepsilon_{(N_{i}-k),i}\\},$

where

$\displaystyle
p(\varepsilon_{(N_{i}-k),i})=\frac{{(N_{i}-k)}\exp(\varepsilon_{(N_{i}-k),i})^{N_{i}-k}}{(1+\exp(\varepsilon_{(N_{i}-k),i}))^{N_{i}-k+1}}$
(A.19)

is the density of a $\mathcal{GL}_{\mbox{\tiny I}}\left(N_{i}-k\right)$
distribution, see (A.4),

$\displaystyle
p(\varepsilon_{(N_{i}-k+1),i})=\frac{k\exp(\varepsilon_{(N_{i}-k+1),i})}{(1+\exp(\varepsilon_{(N_{i}-k+1),i}))^{k+1}}$
(A.20)

is the density of a $\mathcal{GL}_{\mbox{\tiny II}}\left(k\right)$
distribution, see (A.8), and

$\displaystyle C={\frac{\Gamma(N_{i}+1)}{\Gamma(k+1)\Gamma(N_{i}-k+1)}}$

is a normalising constant. It is possible to verify that

$\displaystyle\int_{-\infty}^{+\infty}\int_{\varepsilon_{(N_{i}-k),i}}^{+\infty}p(\varepsilon_{(N_{i}-k),i},\varepsilon_{(N_{i}-k+1),i})\,\text{d}\,\varepsilon_{(N_{i}-k+1),i}\,\text{d}\,\varepsilon_{(N_{i}-k),i}=1.$

Defining
$(\varepsilon_{{v},i},\varepsilon_{{w},i})^{\top}=(\varepsilon_{(N_{i}-k),i},\varepsilon_{(N_{i}-k+1),i})^{\top}$,
yields (A.15).

#### Proof of Theorem 3

Knowing that $y_{i}=k$, with $0<k<N_{i}$, ${w}_{i}|y_{i}$ and ${v}_{i}|y_{i}$
are conditionally independent and the following holds:

$\displaystyle{w}_{i}=\log\lambda_{i}+\varepsilon_{{w},i},\qquad{v}_{i}=\log\lambda_{i}+\varepsilon_{{v},i},$

where $\varepsilon_{{w},i}|y_{i}=k\sim\mathcal{GL}_{\mbox{\tiny
II}}\left(k\right)$ is truncated to $[-\log\lambda_{i},+\infty)$, since
${w}_{i}>0$ and $\varepsilon_{{v},i}|y_{i}=k\sim\mathcal{GL}_{\mbox{\tiny
I}}\left(N_{i}-k\right)$ is truncated to $(-\infty,-\log\lambda_{i}]$, since
${v}_{i}<0$.

Since both $F_{\varepsilon}$ and $F_{\varepsilon}^{-1}$ are available in
closed form for both types of generalized logistic distributions, we obtain:

$\displaystyle\varepsilon_{{w},i}=\log\left(\frac{1}{(1-p)^{1/y_{i}}}-1\right),\quad
1-p=W_{i}(1-F_{\varepsilon}(-\log\lambda_{i}))=W_{i}\left(\frac{\lambda_{i}}{1+\lambda_{i}}\right)^{y_{i}},$

where $W_{i}$ is a uniform random number, see (A.10). This proves equation
(52):

$\displaystyle\varepsilon_{{w},i}=\log\left(\frac{1+\lambda_{i}}{\lambda_{i}}\frac{1}{W_{i}^{1/y_{i}}}-1\right)\Rightarrow{w}_{i}=\log\lambda_{i}+\varepsilon_{{w},i}=\log\left((1+\lambda_{i})\frac{1}{W_{i}^{1/y_{i}}}-\lambda_{i}\right).$

Furthermore,

$\displaystyle\varepsilon_{{v},i}=-\log\left(\frac{1}{p^{1/(N_{i}-y_{i})}}-1\right),\quad
p=V_{i}F_{\varepsilon}(-\log\lambda_{i})=V_{i}\frac{1}{(1+\lambda_{i})^{N_{i}-y_{i}}},$

where $V_{i}$ is a uniform random number, see (A.6). This proves equation
(53):

$\displaystyle\varepsilon_{{v},i}=-\log\left(\frac{1+\lambda_{i}}{V_{i}^{1/(N_{i}-y_{i})}}-1\right)\Rightarrow{v}_{i}=\log\lambda_{i}+\varepsilon_{{v},i}=-\log\left(\frac{1+\lambda_{i}}{\lambda_{i}}\frac{1}{V_{i}^{1/(N_{i}-y_{i})}}-\frac{1}{\lambda_{i}}\right).$

For $y_{i}=0$, only ${v}_{i}|y_{i}$ is sampled; for $y_{i}=N_{i}$, only
${w}_{i}|y_{i}$ is sampled. It is easy to verify that indeed ${w}_{i}>0$ and
${v}_{i}<0$.

###  A.3 Computational details

####  A.3.1 Sampling the utilities in a binary model

Consider the latent variable representation of a binary model

$\mbox{\rm Pr}(y_{i}=1|\lambda_{i})=F_{\varepsilon}(\log\lambda_{i}),$ (A.21)

involving the latent variables $z_{i}$:

$\displaystyle y_{i}=I\\{z_{i}>0\\},$ $\displaystyle\displaystyle
z_{i}=\log\lambda_{i}+\varepsilon_{i},\quad\varepsilon_{i}\sim
f_{\varepsilon}(\varepsilon_{i}),$ (A.22)

where $f_{\varepsilon}(\varepsilon)$ is the pdf of the cdf
$F_{\varepsilon}(\varepsilon)$.
$f_{\varepsilon}(\varepsilon)=\phi(\varepsilon)$ is equal to the standard
normal pdf for a probit model and equal to
$f_{\varepsilon}(\varepsilon)=e^{\varepsilon}/(1+e^{\varepsilon})^{2}$ for a
logit model.

Provided that the quantile function $F_{\varepsilon}^{-1}(p)$ is available in
closed form, it is easy to sample $z_{i}$ from the marginalized density
$z_{i}|\lambda_{i},{\mathbf{y}}$ in the latent variable model (A.22). The
latent variables $z_{i}$ are pairwise independent and the posterior
$z_{i}|y_{i},\lambda_{i}$ is given by

$\displaystyle p(z_{i}|y_{i},\lambda_{i})\propto
p(y_{i}|z_{i})f_{\varepsilon}(z_{i}-\log\lambda_{i}).$

The posterior of $z_{i}$ is $f_{\varepsilon}(z_{i}-\log\lambda_{i})$ truncated
to $(-\infty,0]$, if $y_{i}=0$, and truncated to $[0,\infty)$, if $y_{i}=1$,
hence, $z_{i}=\log\lambda_{i}+\varepsilon_{i}$, where $\varepsilon_{i}\sim
f_{\varepsilon}(\varepsilon_{i})I\\{\varepsilon_{i}>-\log\lambda_{i}\\}$, if
$y_{i}=1$, and $\varepsilon_{i}\sim
f_{\varepsilon}(\varepsilon_{i})I\\{\varepsilon_{i}<-\log\lambda_{i}\\}$, if
$y_{i}=0$.999To simulate $\varepsilon$ from a distribution
$F_{\varepsilon}(\varepsilon)$ truncated to $[a,b]$ we simulate a uniform
random number $U$ and define either
$\varepsilon=F_{\varepsilon}^{-1}(F_{\varepsilon}(a)+U(F_{\varepsilon}(b)-F_{\varepsilon}(a)))$
or
$\varepsilon=F_{\varepsilon}^{-1}(F_{\varepsilon}(b)-U(F_{\varepsilon}(b)-F_{\varepsilon}(a)))$.
Therefore:

$\displaystyle
z_{i}=\log\lambda_{i}+F_{\varepsilon}^{-1}(y_{i}+U_{i}(1-y_{i}-\pi_{i})),$
(A.23)

where $U_{i}\sim\mathcal{U}\left[0,1\right]$ and $\pi_{i}=\mbox{\rm
Pr}(y_{i}=1|\lambda_{i})=F_{\varepsilon}(\log\lambda_{i})$, where
$F_{\varepsilon}^{-1}(p)=\Phi^{-1}(p)$ for the probit model and
$F_{\varepsilon}^{-1}(p)=\log p-\log(1-p)$ for the logit model.

####  A.3.2 Adaptive Metropolis Hastings

Johndrow et al. (2019) suggest to use a simplified adaptive MH algorithm based
on the ideas outlined in haario2001adaptive. Corresponding to their sampling
scheme, we update the probit/logit/multinomial coefficients $\theta$ using a
time-inhomogenous proposal $\theta_{i}^{\star}$ for each component
$\theta_{i}$ of $\theta$ using the proposal kernel

$q_{k}(\theta_{ik},\theta_{i}^{\star})=\phi((\theta_{i}^{\star}-\theta_{ik})^{2}/\sigma_{ik}),\quad\sigma_{ik}=\frac{s}{k}\sum^{k-1}_{j=0}(\theta_{ij}-\bar{\theta}_{ik})^{2},$
(A.24)

where $\phi$ is the standard univariate Gaussian density and
$\bar{\theta}_{ik}$ is the average posterior draw across all previous Gibbs
iterations. That is, they use the variance of all previous Gibbs iterations to
construct on-the-fly proposal densities. For the first 100 iterations, we use
a proposal density with unit variance.

####  A.3.3 Sampling the scale parameter $\delta$ for binomial models

The likelihood contribution for each $y_{i}$ in the boosted Pólya-Gamma
mixture representation of a binomial model (57) reads:

$\displaystyle p(\tilde{w}_{i}|\delta,\log\tilde{\lambda}_{i})\propto$ (A.25)
$\displaystyle\frac{1}{\sqrt{\delta}}\exp\left\\{-\frac{\omega_{{w},i}\tilde{w}_{i}^{2}}{2\delta}+\frac{\tilde{w}_{i}\kappa_{{w},i}}{\sqrt{\delta}}-\frac{\omega_{{w},i}(\log\tilde{\lambda}_{i})^{2}}{2\delta}+\frac{\log\tilde{\lambda}_{i}}{\delta}\left(\tilde{w}_{i}\omega_{{w},i}-\sqrt{\delta}\kappa_{{w},i}\right)\right\\},\quad
y_{i}>0,$ $\displaystyle
p(\tilde{v}_{i}|\delta,\log\tilde{\lambda}_{i})\propto$
$\displaystyle\frac{1}{\sqrt{\delta}}\exp\left\\{-\frac{\omega_{{v},i}\tilde{v}_{i}^{2}}{2\delta}+\frac{\tilde{v}_{i}\kappa_{{v},i}}{\sqrt{\delta}}-\frac{\omega_{{v},i}(\log\tilde{\lambda}_{i})^{2}}{2\delta}+\frac{\log\tilde{\lambda}_{i}}{\delta}\left(\tilde{v}_{i}\omega_{{v},i}-\sqrt{\delta}\kappa_{{v},i}\right)\right\\},\quad
y_{i}<N_{i}.$

The likelihood contributions (A.25) have to be combined with the prior
$p(\delta|\gamma)$ and with the prior distributions of the parameters in
$\log\tilde{\lambda}_{i}$. We provide details for a binomial regression model.

For a binomial regression model,
$\log\tilde{\lambda}_{i}={\mathbf{x}}_{i}\tilde{\boldsymbol{\beta}}$, where
the parameter $\tilde{\boldsymbol{\beta}}$ is defined exactly as for a binary
regression model:

$\displaystyle\tilde{\beta}_{j}=\left\\{\begin{array}[]{ll}\sqrt{\delta}\beta_{j},&1\leq
j<d,\\\ \sqrt{\delta}\beta_{j}+\gamma,&j=d.\\\ \end{array}\right.$ (A.28)

As in the binary case, we assume the following prior for
$\tilde{\boldsymbol{\beta}}$ in the boosted model:

$\displaystyle\tilde{\boldsymbol{\beta}}|\delta\sim\mathcal{N}_{d}\left({\mathbf{0}},\delta\boldsymbol{A}_{0}\right),\qquad$
(A.29)

where $A_{0,dd}$ is possibly larger than $A_{0,jj}$, for $1\leq j<d$, which
implies the following joint prior of $\boldsymbol{\beta}|\gamma,\delta$ in the
identified model:

$\displaystyle\boldsymbol{\beta}|\gamma,\delta\sim\mathcal{N}_{d}\left(\boldsymbol{a}_{0},\boldsymbol{A}_{0}\right),\qquad\boldsymbol{a}_{0}=\left(\begin{array}[]{c}{\mathbf{0}}_{(d-1)\times
1}\\\ -\frac{\gamma}{\sqrt{\delta}}\end{array}\right).$ (A.32)

This dependence on $\delta$ and $\gamma$ has to be taken into account, when we
move from the identified to the unidentified model and sample
$\delta^{\star}\sim\mathcal{G}^{-1}\left(d_{0},D_{0}\right)$ and sample
$\gamma^{\star}|\delta^{\star},{\beta}_{d}$ from

$\displaystyle\gamma^{\star}|\delta^{\star},\beta_{d}\sim\mathcal{N}\left(\sqrt{\delta^{\star}}g_{1},\delta^{\star}G_{1}\right),$
$\displaystyle G_{1}=G_{0}A_{0,dd}/(G_{0}+A_{0,dd}),\quad
g_{1}=-\beta_{d}\times G_{0}/(G_{0}+A_{0,dd}).$

The boosted model reads

$\displaystyle\tilde{w}_{i}={\mathbf{x}}_{i}\tilde{\boldsymbol{\beta}}+\sqrt{\delta}\kappa_{{w},i}/\omega_{{w},i}+\tilde{\varepsilon}_{{w},i},\quad
y_{i}>0,$
$\displaystyle\tilde{v}_{i}={\mathbf{x}}_{i}\tilde{\boldsymbol{\beta}}+\sqrt{\delta}\kappa_{{v},i}/\omega_{{v},i}+\tilde{\varepsilon}_{{v},i},\quad
y_{i}<N_{i}.$

with $\tilde{\varepsilon}_{{w},i}=\sqrt{\delta}\varepsilon_{{w},i}$ and
$\tilde{\varepsilon}_{{v},i}=\sqrt{\delta}\varepsilon_{{v},i}$. The choice
equation is defined in (61).

To obtain the joint posterior
$p(\tilde{\boldsymbol{\beta}},\delta|\tilde{\boldsymbol{z}},\gamma)$, we
combine the likelihood contributions (A.25) with the conditional prior
$p(\delta|\gamma)$ given in (8) and the prior distribution
$p(\tilde{\boldsymbol{\beta}}|\delta)$ given in (A.29):

$\displaystyle
p(\tilde{\boldsymbol{\beta}},\delta|\tilde{\boldsymbol{z}},\gamma)\propto\prod_{i:y_{i}>0}p(\tilde{w}_{i}|\delta,\tilde{\boldsymbol{\beta}})\prod_{i:y_{i}<N_{i}}p(\tilde{v}_{i}|\delta,\tilde{\boldsymbol{\beta}})p(\tilde{\boldsymbol{\beta}}|\delta)p(\delta|\gamma).$
(A.34)

It can be shown that the joint posterior
$p(\tilde{\boldsymbol{\beta}},\delta|\tilde{\boldsymbol{z}},\gamma)$
factorizes in the following way

$p(\tilde{\boldsymbol{\beta}},\delta|\tilde{\boldsymbol{z}},\gamma)\propto
f_{N}(\tilde{\boldsymbol{\beta}};{\mathbf{b}}_{N}(\delta),\delta{\mathbf{B}}_{N})\left(\frac{1}{\delta}\right)^{d_{I}+1}\exp\left\\{-\frac{D_{I}}{\delta}\right\\}\exp\left\\{\frac{B_{I}}{\sqrt{\delta}}\right\\}$

with

$\displaystyle{\mathbf{b}}_{N}(\delta)={\mathbf{B}}_{N}\boldsymbol{m}_{N}(\delta),\,\,\boldsymbol{m}_{N}(\delta)=\sum_{i=1}^{N}m_{i}{\mathbf{x}}_{i}^{\top},\qquad{\mathbf{B}}_{N}^{-1}=\boldsymbol{A}_{0}^{-1}+\sum_{i=1}^{N}M_{i}{\mathbf{x}}_{i}^{\top}{\mathbf{x}}_{i},$
(A.35) $\displaystyle
m_{i}=I\\{y_{i}>0\\}\left(\tilde{w}_{i}\omega_{{w},i}-\sqrt{\delta}\kappa_{{w},i}\right)+I\\{y_{i}<N_{i}\\}\left(\tilde{v}_{i}\omega_{{v},i}-\sqrt{\delta}\kappa_{{v},i}\right),$
$\displaystyle
M_{i}=I\\{y_{i}>0\\}\omega_{{w},i}+I\\{y_{i}<N_{i}\\}\omega_{{v},i},$
$\displaystyle
d_{I}=d_{0}+1/2+\frac{1}{2}\sum_{i=1}^{N}\left(I\\{y_{i}>0\\}+I\\{y_{i}<N_{i}\\}\right),$
(A.36) $\displaystyle
D_{I}=D_{0}+\gamma^{2}/(2G_{0})+\frac{1}{2}\sum_{i=1}^{N}\left(I\\{y_{i}>0\\}\tilde{w}_{i}^{2}\omega_{{w},i}+I\\{y_{i}<N_{i}\\}\tilde{v}_{i}^{2}\omega_{{v},i}\right)-\frac{1}{2}\boldsymbol{m}_{a}^{\top}{\mathbf{B}}_{N}\boldsymbol{m}_{a},$
$\displaystyle
B_{I}=\sum_{i=1}^{N}\left(I\\{y_{i}>0\\}\tilde{w}_{i}\kappa_{{w},i}+I\\{y_{i}<N_{i}\\}\tilde{v}_{i}\kappa_{{v},i}\right)-\boldsymbol{m}_{a}^{\top}{\mathbf{B}}_{N}\boldsymbol{m}_{b},$

where
$f_{N}(\tilde{\boldsymbol{\beta}};{\mathbf{b}}_{N}(\delta),\delta{\mathbf{B}}_{N})$
is the density of a
$\mathcal{N}_{d}\left({\mathbf{b}}_{N}(\delta),\delta{\mathbf{B}}_{N}\right)$-distribution
and $\boldsymbol{m}_{a}$ and $\boldsymbol{m}_{b}$ are defined below.

The last terms in $D_{I}$ and $B_{I}$ result from completing the squares in
$\tilde{\boldsymbol{\beta}}$. Note that

$\displaystyle
f_{N}(\tilde{\boldsymbol{\beta}};{\mathbf{b}}_{N}(\delta),\delta{\mathbf{B}}_{N})=\frac{1}{\sqrt{2\pi\delta}^{d}\sqrt{|{\mathbf{B}}_{N}|}}\exp\left\\{-\frac{1}{2\delta}(\tilde{\boldsymbol{\beta}}-{\mathbf{b}}_{N}(\delta))^{\top}{\mathbf{B}}_{N}^{-1}(\tilde{\boldsymbol{\beta}}-{\mathbf{b}}_{N}(\delta))\right\\}$
$\displaystyle=\frac{1}{\sqrt{2\pi\delta}^{d}\sqrt{|{\mathbf{B}}_{N}|}}\exp\left\\{-\frac{1}{2\delta}\tilde{\boldsymbol{\beta}}^{\top}{\mathbf{B}}_{N}^{-1}\tilde{\boldsymbol{\beta}}+\frac{1}{\delta}\boldsymbol{m}_{N}(\delta)^{\top}\tilde{\boldsymbol{\beta}}\right\\}\exp\left\\{-\frac{1}{2\delta}\boldsymbol{m}_{N}(\delta)^{\top}{\mathbf{B}}_{N}\boldsymbol{m}_{N}(\delta)\right\\}.$

Hence, the factor
$\exp\left\\{\frac{1}{2\delta}\boldsymbol{m}_{N}(\delta)^{\top}{\mathbf{B}}_{N}\boldsymbol{m}_{N}(\delta)\right\\}$
has to be added to (A.34) to compensate for completing the squares in
$f_{N}(\tilde{\boldsymbol{\beta}};{\mathbf{b}}_{N}(\delta),\delta{\mathbf{B}}_{N})$.
Using that

$\displaystyle\boldsymbol{m}_{N}(\delta)=\boldsymbol{m}_{a}-\sqrt{\delta}\boldsymbol{m}_{b},$
(A.37)
$\displaystyle\boldsymbol{m}_{a}=\sum_{i=1}^{N}\left(I\\{y_{i}>0\\}\tilde{w}_{i}\omega_{{w},i}+I\\{y_{i}<N_{i}\\}\tilde{v}_{i}\omega_{{v},i}\right){\mathbf{x}}_{i}^{\top},$
$\displaystyle\boldsymbol{m}_{b}=\sum_{i=1}^{N}\left(I\\{y_{i}>0\\}\kappa_{{w},i}+I\\{y_{i}<N_{i}\\}\kappa_{{v},i}\right){\mathbf{x}}_{i}^{\top},$

we obtain:

$\displaystyle\exp\left\\{\frac{1}{2\delta}\boldsymbol{m}_{N}(\delta)^{\top}{\mathbf{B}}_{N}\boldsymbol{m}_{N}(\delta)\right\\}\propto\exp\left\\{\frac{1}{2\delta}\boldsymbol{m}_{a}^{\top}{\mathbf{B}}_{N}\boldsymbol{m}_{a}\right\\}\exp\left\\{-\frac{1}{\sqrt{\delta}}\boldsymbol{m}_{a}^{\top}{\mathbf{B}}_{N}\boldsymbol{m}_{b}\right\\}.$

The marginal density
$p(\delta|\gamma,\tilde{\boldsymbol{z}},{\mathbf{y}})\propto\left(\frac{1}{\delta}\right)^{d_{I}+1}e^{-D_{I}/\delta}e^{B_{I}/\sqrt{\delta}}$
is not available in closed form, unless all $N_{i}=1$. In this case
$\kappa_{{w},i}=\kappa_{{v},i}=0$ for all $i=1,\dots,N$ and hence $B_{I}=0$
and $p(\delta|\gamma,\tilde{\boldsymbol{z}},{\mathbf{y}})$ follows the
inverted gamma distribution $\mathcal{G}^{-1}\left(d_{I},D_{I}\right)$. It can
be shown that in any case $D_{I}>0$, provided that $D_{0}>0$.

Figure A.1: Approximating $p(\delta|\cdot)$ by the inverted gamma densities
$\delta\sim\mathcal{G}^{-1}\left(d_{I}^{\star},D_{I}^{\star}\right)$ (left-
hand side) and
$\sqrt{\delta}\sim\mathcal{G}^{-1}\left(b_{I}^{\star},B_{I}^{\star}\right)$
(right-hand side) for a case where $B_{I}>0$ and $B_{I}<0$ (bottom).

To sample from

$p(\delta|\gamma,\tilde{\boldsymbol{z}},{\mathbf{y}})\propto\left(\frac{1}{\delta}\right)^{d_{I}+1}e^{-D_{I}/\delta}e^{B_{I}/\sqrt{\delta}}$

in Algorithm 3 the following resampling technique is used: Choose an
“auxiliary prior” $\pi(\delta)$ for resampling such that the mode and the
curvature of $\pi(\delta)$ coincide with the mode $\delta_{M}$ and the
curvature $I_{p}$ of the posterior
$p(\delta|\gamma,\tilde{\boldsymbol{z}},{\mathbf{y}})$ which are given by:

$\displaystyle\delta_{M}=\frac{16D_{I}^{2}}{\left[B_{I}+\sqrt{B_{I}^{2}+16D_{I}(d_{I}+1)}\right]^{2}},\qquad
I_{p}=-\frac{\sqrt{B_{I}^{2}+16D_{I}(d_{I}+1)}}{4\cdot(\delta_{M})^{\frac{5}{2}}}.$

Resampling works as follows. $L$ draws $\delta^{(l)}\sim\pi(\delta)$,
$l=1,\ldots,L$, from the auxiliary prior are resampled using weights
proportional to the “auxiliary likelihood”
$\ell(\delta)=p(\delta|\gamma,\tilde{\boldsymbol{z}},{\mathbf{y}})/\pi(\delta)$,
given by:

$\displaystyle\log\ell(\delta)$ $\displaystyle\propto$ $\displaystyle\log
p(\delta|\gamma,\tilde{\boldsymbol{z}},{\mathbf{y}})-\log\pi(\delta)$
$\displaystyle\propto$
$\displaystyle-(d_{I}+1)\log\delta-\frac{D_{I}}{\delta}+\frac{B_{I}}{\sqrt{\delta}}-\log\pi(\delta).$

The desired draw from $p(\delta|\cdot)$ is given by $\delta^{(l^{\star})}$,
where $l^{\star}\sim\mbox{\rm MulNom}\left(1;w_{1},\ldots,w_{L}\right)$ and
the weights $w_{l}\propto\ell(\delta^{(l)})$ are normalized to 1.

The auxiliary likelihood $\ell(\delta)$ is expected to be rather flat over the
support of $\pi(\delta)$, see Figure A.1 for illustration. Hence,
$w_{1},\ldots,w_{L}$ is expected to be close to a uniform distribution and $L$
can be pretty small ($L=5$ or $L=10$ should be enough).

The factorization of the posterior suggests two distribution families as
auxiliary prior $\pi(\delta)$. First, the inverted Gamma prior
$\delta\sim\mathcal{G}^{-1}\left(d_{I}^{\star},D_{I}^{\star}\right)$ with mode
$\delta_{IG}$ and curvature $I_{IG}$ of the pdf given by:

$\displaystyle\displaystyle\delta_{IG}=\frac{D_{I}^{\star}}{d_{I}^{\star}+1},\qquad
I_{IG}=-\frac{(d_{I}^{\star}+1)^{3}}{(D_{I}^{\star})^{2}}.$

Matching the mode, i.e. $\delta_{IG}=\delta_{M}$, and the curvature, i.e.
$I_{IG}=I_{p}$, to the posterior, i.e.

$\displaystyle
I_{IG}=-\frac{(d_{I}^{\star}+1)^{3}}{(D_{I}^{\star})^{2}}=-\frac{d_{I}^{\star}+1}{(\delta_{IG})^{2}}=-\frac{d_{I}^{\star}+1}{(\delta_{M})^{2}}=I_{p},$
$\displaystyle\delta_{IG}=\frac{D_{I}^{\star}}{d_{I}^{\star}+1}=\delta_{M},$

yields following optimal choice for the parameters
$(d_{I}^{\star},D_{I}^{\star})$:

$\displaystyle d_{I}^{\star}=-I_{p}\cdot(\delta_{M})^{2}-1,\quad
D_{I}^{\star}=\delta_{M}(d_{I}^{\star}+1).$

The log likelihood ratio reads:

$\displaystyle\log\ell(\delta)\propto-(d_{I}-d_{I}^{\star})\log\delta-\frac{D_{I}-D_{I}^{\star}}{\delta}+\frac{B_{I}}{\sqrt{\delta}}.$

Second, provided that $B_{I}^{\star}>0$, the translated inverted Gamma prior
$\sqrt{\delta}\sim\mathcal{G}^{-1}\left(2b_{I}^{\star},B_{I}^{\star}\right)$
101010As follows from the law of transformation of densities, a random
variable $\delta$, where the transformed variable
$\sqrt{\delta}\sim\mathcal{G}^{-1}\left(2a,b\right)$, has density
$\displaystyle
p_{SQ}(\delta;2a,b)=\frac{b^{2a}}{2\Gamma(2a)}\left(\frac{1}{\delta}\right)^{a+1}e^{-\frac{b}{\sqrt{\delta}}}.$
(A.38) To sample from such a density, we sample
$X\sim\mathcal{G}^{-1}\left(2a,b\right)$ and take the square, i.e.
$\delta=X^{2}$. with mode $\delta_{SQ}$ and the curvature $I_{SQ}$ of the pdf
given by:

$\displaystyle\displaystyle\delta_{SQ}=\frac{(B_{I}^{\star})^{2}}{4(b_{I}^{\star}+1)^{2}},\qquad
I_{SQ}=-\frac{8(b_{I}^{\star}+1)^{5}}{(B_{I}^{\star})^{4}}.$

Matching the mode, i.e. $\delta_{SQ}=\delta_{M}$, and the curvature, i.e.
$I_{SQ}=I_{p}$, to the posterior, i.e.

$\displaystyle
I_{SQ}=-\frac{8(b_{I}^{\star}+1)^{5}}{(B_{I}^{\star})^{4}}=-\frac{b_{I}^{\star}+1}{2(\delta_{SQ})^{2}}=-\frac{b_{I}^{\star}+1}{2(\delta_{M})^{2}}=I_{p},$
$\displaystyle\delta_{SQ}=\frac{(B_{I}^{\star})^{2}}{4(b_{I}^{\star}+1)^{2}}=\delta_{M},$

yields following optimal choice for the parameters
$(b_{I}^{\star},b_{I}^{\star})$:

$\displaystyle b_{I}^{\star}=-2I_{p}\cdot(\delta_{M})^{2}-1,\quad
B_{I}^{\star}=2(b_{I}^{\star}+1)\sqrt{\delta_{M}}.$

The log likelihood ratio reads:

$\displaystyle\log\ell(\delta)\propto-(d_{I}-b_{I}^{\star})\log\delta-\frac{D_{I}}{\delta}+\frac{B_{I}+B_{I}^{\star}}{\sqrt{\delta}}.$

See again Figure A.1 for the optimal values for two examples, one with
$B_{I}>0$ and one with $B_{I}<0$.

###  A.4 More Details on the UPG sampler for binary state space models

As for binary regression models, we use the priors
$\gamma|\delta\sim\mathcal{N}\left(0,\delta
G_{0}\right),\delta\sim\mathcal{G}^{-1}\left(d_{0},D_{0}\right)$ and, specify
the following priors in the boosted SSM (67):

$\displaystyle\tilde{\beta}_{j0}|\tilde{\beta}_{j},\theta_{j},\delta\sim\mathcal{N}\left(\tilde{\beta}_{j},\delta
P_{0,{jj}}\theta_{j}\right),\quad\tilde{\beta}_{j}|\delta\sim\mathcal{N}\left(0,\delta
A_{0,jj}\right),\quad\theta_{j}\sim\mathcal{G}^{-1}\left(c_{0},C_{0}\right),$
(A.39)

where $A_{0,jj},j=1,\ldots,d$ and $c_{0}$ as well as $C_{0}$ are fixed
hyperparameters. Marginalizing with respect to $\tilde{\beta}_{j}$ yields
following prior for the initial value $\tilde{\beta}_{j0}$ in (67):

$\displaystyle\tilde{\beta}_{j0}|\delta,\theta_{j}\sim\mathcal{N}\left(0,\delta(A_{0,jj}+P_{0,{jj}}\theta_{j})\right).$
(A.40)

Matching back through transformation (68) yields the prior in the original SSM
(65):

$\displaystyle\beta_{j}\sim\mathcal{N}\left(0,A_{0,jj}\right),j<d,\quad\beta_{d}|\gamma,\delta\sim\mathcal{N}\left(-\frac{\gamma}{\sqrt{\delta}},A_{0,dd}\right).$
(A.41)

We observe the same prior dependence between $\gamma$, $\delta$ and
$\beta_{d}$ as for the simple binary regression model:

$\displaystyle\gamma|\delta,\beta_{d}\sim\mathcal{N}\left(\sqrt{\delta}g_{1},\delta
G_{1}\right),\quad G_{1}=G_{0}A_{0,dd}/(G_{0}+A_{0,dd}),\quad
g_{1}=-\beta_{d}\times G_{0}/(G_{0}+A_{0,dd}).$ (A.42)

Note that in Step (b-2a) of Algorithm 4, $\delta^{\mathrm{new}}$ is sampled
from the marginalized posterior
$\delta^{\mathrm{new}}|\gamma^{\mathrm{new}},\tilde{\boldsymbol{z}}$,
$\theta_{1},\ldots,\theta_{d}$ where the entire state process and the initial
values are integrated out. This is possible, since all process variances, the
observation variance, and the variance of the initial state vector are scaled
with $\delta$, see Appendix A.4.1.

Algorithm 4 The ultimate Pólya-Gamma sampler for binary state space models.

Choose starting values for $(\theta_{1},\ldots,\theta_{d})$, and
$\lambda_{1},\ldots,\lambda_{T}$ and repeat the following steps:

  * (a-1)

For each $t=1,\ldots,T$, sample
$z_{t}=\log\lambda_{t}+F_{\varepsilon}^{-1}(y_{t}+U_{t}(1-y_{t}-\pi_{t}))$ in
the SSM (65), where $U_{t}\sim\mathcal{U}\left[0,1\right]$ and
$\pi_{t}=F_{\varepsilon}(\log\lambda_{t})$. Sample
$\omega_{t}|z_{t},\lambda_{t},y_{t}\sim\mathcal{PG}\left(2,|z_{t}-\log\lambda_{t}|\right)$
for the logit SSM.

  * (a-2)

Sample $\delta^{\star}\sim\mathcal{G}^{-1}\left(d_{0},D_{0}\right)$ and
$\gamma^{\star}|\delta^{\star}$ from (A.42). Move to the boosted SSM (67)
using transformation (68) to define $\tilde{z}_{t}$, $t=1,\ldots,T$.

  * (b)

Sampling in the boosted conditionally Gaussian SSM (67):

    * (b-1)

Sample $\gamma^{\mathrm{new}}=\sqrt{G_{0}^{\star}}T_{2d_{0}}^{-1}(U\times
T_{2d_{0}}(O(\tilde{\boldsymbol{z}})/\sqrt{G_{0}^{\star}})+(1-U)\times
T_{2d_{0}}(L(\tilde{\boldsymbol{z}})/\sqrt{G_{0}^{\star}}))$, where
$U\sim\mathcal{U}\left[0,1\right]$, $T_{\nu}(\cdot)$ is the cdf of a $t_{\nu}$
distribution, $L(\tilde{\boldsymbol{z}})=\max_{t:y_{t}=0}\tilde{z}_{t}$ and
$O(\tilde{\boldsymbol{z}})=\min_{t:y_{t}=1}\tilde{z}_{t}$.

    * (b-2)

Sample $\delta^{\mathrm{new}}$ and
$\tilde{\boldsymbol{\beta}}_{0},\ldots,\tilde{\boldsymbol{\beta}}_{T}$ jointly
conditional on $\theta_{1},\ldots,\theta_{d}$ based on the marginal prior
$\tilde{\beta}_{j0}|\delta,\theta_{j}$ given in (A.40), using the methods
described in Appendix A.4.1 with the process variances being equal to
$\delta\theta_{j}$ and heteroscedastic observation errors
$\sigma_{t}^{2}=\delta/\omega_{t}$ (where $\omega_{t}$ is known):

      * (b-2a)

Sample $\delta^{\mathrm{new}}$ from
$\delta{|\gamma^{\mathrm{new}},\tilde{\boldsymbol{z}},\theta_{1},\ldots,\theta_{d}}\sim\mathcal{G}^{-1}\left(d_{T},D_{T}\right)$,
where $d_{T}$ and $D_{T}$ are given by:

$\displaystyle{d_{T}}={d_{0}+\frac{T+1}{2}},\quad{D_{T}}={D_{0}+\frac{(\gamma^{\mathrm{new}})^{2}}{2G_{0}}+\frac{1}{2}\sum_{t=1}^{T}\frac{(\tilde{z}_{t}-\hat{z}_{t|t-1})^{2}}{S_{t|t-1}}},$

and $\hat{z}_{t|t-1}$ and $S_{t|t-1}$ are obtained from a Kalman filter with
$\delta=1$, see (A.48).

      * (b-2b)

Sample
$\tilde{\boldsymbol{\beta}}_{0},\ldots,\tilde{\boldsymbol{\beta}}_{T}|{\delta^{\mathrm{new}},\tilde{\boldsymbol{z}},\theta_{1},\ldots,\theta_{d}}$
using forward-filtering-backward-sampling.

    * (b-3)

Update $(\tilde{\beta}_{j},\theta_{j})|\delta^{\mathrm{new}}$ under the prior
$\tilde{\beta}_{j0}|\tilde{\beta}_{j},\theta_{j},\delta^{\mathrm{new}}\sim\mathcal{N}\left(\tilde{\beta}_{j},\delta^{\mathrm{new}}P_{0,{jj}}\theta_{j}\right)$
given by (A.39):

      * (b-3a)

Sample
$\tilde{\beta}_{j}|\theta_{j},\tilde{\beta}_{j0},\delta^{\mathrm{new}},$ for
$j=1,\ldots,d$:

$\displaystyle\tilde{\beta}_{j}|\theta_{j},\delta^{\mathrm{new}},\tilde{\beta}_{j0}\sim\mathcal{N}\left(\frac{A_{0,jj}\tilde{\beta}_{j0}}{A_{0,jj}+\theta_{j}P_{0,{jj}}},\delta^{\mathrm{new}}\frac{A_{0,jj}\theta_{j}P_{0,{jj}}}{A_{0,jj}+\theta_{j}P_{0,{jj}}}\right).$

      * (b-3b)

Sample
$\theta_{j}|\delta^{\mathrm{new}},\tilde{\beta}_{j},\\{\tilde{\beta}_{jt}\\}\sim\mathcal{G}^{-1}\left(c_{0}+(T+1)/2,C_{j}\right)$
for $j=1,\ldots,d$ where

$\displaystyle
C_{j}=C_{0}+\frac{1}{2\delta^{\mathrm{new}}}\left(\frac{(\tilde{\beta}_{j0}-\tilde{\beta}_{j})^{2}}{P_{0,{jj}}}+\sum_{t=1}^{T}(\tilde{\beta}_{jt}-\tilde{\beta}_{j,t-1})^{2}\right).$

    * (b-4)

Move back to the identified model through the inverse of transformation (68),
i.e. $j=1,\ldots,d$,
$\beta_{{j}}=(\tilde{\beta}_{j}-I\\{j=d\\}\gamma^{\mathrm{new}})/\sqrt{\delta^{\mathrm{new}}}$
and
$\beta_{{jt}}=(\tilde{\beta}_{jt}-I\\{j=d\\}\gamma^{\mathrm{new}})/\sqrt{\delta^{\mathrm{new}}}$,
$t=0,\ldots,T$.

####  A.4.1 Linear Gaussian state space models with unknown scaling

A linear Gaussian state space model with unknown scaling is defined as the
following dynamic stochastic system for $t=1,\ldots,T$:

$\displaystyle{\boldsymbol{\beta}}_{t}={\mathbf{F}}_{t}{\boldsymbol{\beta}}_{t-1}+{\mathbf{w}}_{t},\qquad{\mathbf{w}}_{t}\sim\mathcal{N}_{d}\left({\mathbf{0}},\delta{\mathbf{Q}}_{t}\right),$
(A.43)
$\displaystyle{\mathbf{y}}_{t}={\mathbf{H}}_{t}{\boldsymbol{\beta}}_{t}+\boldsymbol{\epsilon}_{t},\qquad\boldsymbol{\epsilon}_{t}\sim\mathcal{N}_{r}\left({\mathbf{0}},\delta{\mathbf{R}}_{t}\right),$
(A.44)

where
${\boldsymbol{\beta}}_{0}\sim\mathcal{N}_{d}\left(\hat{\boldsymbol{\beta}}_{0|0},\delta{\mathbf{P}}_{0|0}\right)$.
The key variables in these formulations are the state variable
${\boldsymbol{\beta}}_{t}$ and the observation variable ${\mathbf{y}}_{t}$.
The state variable ${\boldsymbol{\beta}}_{t}$ is a latent $d$-dimensional
random vector which is observed only indirectly through the effect it has on
the distribution of the $r$-dimensional time series ${\mathbf{y}}_{t}$. Note
that the system (A.43) and (A.44) is a linear Gaussian state space model where
all (co)variances depend on the same unknown scaling factor $\delta$. Let
${\mathbf{\boldsymbol{\vartheta}}}$ be the vector containing all unknown
parameters in the system matrices ${\mathbf{F}}_{t}$, ${\mathbf{H}}_{t}$,
${\mathbf{Q}}_{t}$ and ${\mathbf{R}}_{t}$ apart from $\delta$.

In filtering, the goal is to estimate
${\boldsymbol{\beta}}_{0},\ldots,{\boldsymbol{\beta}}_{T}|{\mathbf{y}},\delta,{\mathbf{\boldsymbol{\vartheta}}}$
given data ${\mathbf{y}}=({\mathbf{y}}_{1},\ldots,{\mathbf{y}}_{T})$ for fixed
model parameters ${\mathbf{\boldsymbol{\vartheta}}}$. For a linear Gaussian
state space model with unknown scaling, the filter density
$p({\boldsymbol{\beta}}_{t}|\delta,{\mathbf{\boldsymbol{\vartheta}}},{\mathbf{y}}^{t})$
at time $t$, where
${\mathbf{y}}^{t}=({\mathbf{y}}_{1},\ldots,{\mathbf{y}}_{t})$, is equal to a
normal distribution, scaled with $\delta$, i.e.

$\displaystyle{\boldsymbol{\beta}}_{t}|\delta,{\mathbf{\boldsymbol{\vartheta}}},{\mathbf{y}}^{t}\sim\mathcal{N}_{d}\left(\hat{\boldsymbol{\beta}}_{t|t},\delta\,{\mathbf{P}}_{t|t}\right),$

where $\hat{\boldsymbol{\beta}}_{t|t}$ and ${\mathbf{P}}_{t|t}$ are determined
recursively for $t=1,\ldots,T$ from $\hat{\boldsymbol{\beta}}_{t-1|t-1},$
${\mathbf{P}}_{t-1|t-1}$ and ${\mathbf{y}}_{t}$. This results follows
immediately from applying the Kalman filter to the system (A.43) and (A.44).
From a computational point of view, this extended filter is a standard Kalman-
filter which is run with $\delta=1$, and where all covariances in the filter
density
${\boldsymbol{\beta}}_{t}|\delta,{\mathbf{\boldsymbol{\vartheta}}},{\mathbf{y}}^{t}$,
the propagated density
${\boldsymbol{\beta}}_{t}|\delta,{\mathbf{\boldsymbol{\vartheta}}},{\mathbf{y}}^{t-1}$
and the predictive density
$y_{t}|\delta,{\mathbf{\boldsymbol{\vartheta}}},{\mathbf{y}}^{t-1}$ are
multiplied by $\delta$ after running the filter.

We provide further details for TVP models for univariate time series $y_{t}$,
where $r=1$, ${\mathbf{F}}_{t}={{\mathbf{I}}}_{d}$,
${\mathbf{H}}_{t}={\mathbf{x}}_{t}$,
${\mathbf{Q}}_{t}={\mathbf{{\mathbf{Q}}}}=\mbox{\rm
Diag}\left(\theta_{1},\ldots,\theta_{d}\right)$,
${\mathbf{R}}_{t}=\sigma_{t}^{2}$,
$\hat{\boldsymbol{\beta}}_{0|0}=(\beta_{1},\ldots,\beta_{d})^{\prime}$ and
${\mathbf{P}}_{0|0}=\mbox{\rm
Diag}\left(\theta_{1}P_{0,{11}},\ldots,\theta_{d}P_{0,{dd}}\right)$:

  1. (a)

Propagation:
${\boldsymbol{\beta}}_{t}|\delta,{\mathbf{\boldsymbol{\vartheta}}},{\mathbf{y}}^{t-1}\sim\mathcal{N}_{d}\left(\hat{\boldsymbol{\beta}}_{t|t-1},\delta{\mathbf{P}}_{t|t-1}\right),$
where

$\displaystyle\hat{\boldsymbol{\beta}}_{t|t-1}=\hat{\boldsymbol{\beta}}_{t-1|t-1},\quad{{\mathbf{P}}_{t|t-1}}={\mathbf{P}}_{t-1|t-1}+{{\mathbf{{\mathbf{Q}}}}}.$

  2. (b)

Prediction:
$y_{t}|\delta,{\mathbf{\boldsymbol{\vartheta}}},{\mathbf{y}}^{t-1}\sim\mathcal{N}\left(\hat{y}_{t|t-1},\delta
S_{t|t-1}\right)$,

$\displaystyle\hat{y}_{t|t-1}={\mathbf{x}}_{t}\hat{\boldsymbol{\beta}}_{t|t-1},\quad
S_{t|t-1}={\mathbf{x}}_{t}{\mathbf{P}}_{t|t-1}{\mathbf{x}}_{t}^{\top}+{\sigma_{t}^{2}}.$
(A.45)

  3. (c)

Correction:
${\boldsymbol{\beta}}_{t}|\delta,{\mathbf{\boldsymbol{\vartheta}}},{\mathbf{y}}^{t}\sim\mathcal{N}_{d}\left(\hat{\boldsymbol{\beta}}_{t|t},\delta{\mathbf{P}}_{t|t}\right),$

$\displaystyle\hat{\boldsymbol{\beta}}_{t|t}=\hat{\boldsymbol{\beta}}_{t|t-1}+{{\mathbf{K}}_{t}}(y_{t}-\hat{y}_{t|t-1}),\quad{\mathbf{P}}_{t|t}=({\mathbf{I}}-{\mathbf{K}}_{t}{\mathbf{x}}_{t}){\mathbf{P}}_{t|t-1},$
$\displaystyle{{\mathbf{K}}_{t}}={{\mathbf{P}}_{t|t-1}}{\mathbf{x}}_{t}^{\top}S_{t|t-1}^{-1}.$

An advantage of writing the Kalman filter in this way is that the integrated
likelihood $p({\mathbf{y}}|\delta,{\mathbf{\boldsymbol{\vartheta}}})$ is
available in closed form. If $r=1$, then exploiting (A.45) we obtain:

$\displaystyle
p({\mathbf{y}}|\delta,{\mathbf{\boldsymbol{\vartheta}}})=\prod_{t=1}^{T}p(y_{t}|{\mathbf{y}}^{t-1},\delta,{\mathbf{\boldsymbol{\vartheta}}})\propto\left(\frac{1}{\delta}\right)^{\frac{T}{2}}\exp\left\\{-\frac{1}{2\delta}\sum_{t=1}^{T}\frac{(y_{t}-\hat{y}_{t|t-1})^{2}}{S_{t|t-1}}\right\\},$
(A.46)

where $\hat{y}_{t|t-1}$ and $S_{t|t-1}$ are the moments given in (A.45). This
allows to sample the scaling factor $\delta$ jointly with the state process
${\boldsymbol{\beta}}_{0},\ldots,{\boldsymbol{\beta}}_{T}$ conditional on
${\mathbf{\boldsymbol{\vartheta}}}$ and ${\mathbf{y}}$,

  * (a)

by sampling $\delta|{\mathbf{\boldsymbol{\vartheta}}},{\mathbf{y}}$ from the
inverted Gamma prior that results from combining the likelihood (A.46) with
the inverted gamma prior $\delta\sim\mathcal{G}^{-1}\left(d_{0},D_{0}\right)$,

$\displaystyle\delta|{\mathbf{\boldsymbol{\vartheta}}},{\mathbf{y}}\sim\mathcal{G}^{-1}\left(d_{T},D_{T}\right),$
(A.47) $\displaystyle d_{T}=d_{0}+\frac{T}{2},\quad
D_{T}=D_{0}+\frac{1}{2}\sum_{t=1}^{T}\frac{(y_{t}-\hat{y}_{t|t-1})^{2}}{S_{t|t-1}};$

  * (b)

by sampling
${\boldsymbol{\beta}}_{0},\ldots,{\boldsymbol{\beta}}_{T}|{\mathbf{y}},\delta,{\mathbf{\boldsymbol{\vartheta}}}$
conditional on $\delta$ using FFBS (Frühwirth-Schnatter 1994).

As mentioned above, the moments $\hat{y}_{t|t-1}$ and $S_{t|t-1}$ in Step (a)
are obtained from running a Kalman-filter with $\delta=1$. Note that we need
not rerun the Kalman-filter in Step (b), all covariances are just multiplied
by the scaling factor $\delta$ sampled in Step (a).

##### Application to boosting binary state space models.

This result is useful for the boosted state space model (67), because the
working parameter $\delta^{\mathrm{new}}$ can be sampled without conditioning
on the state process
$\tilde{\boldsymbol{\beta}}_{0},\ldots,\tilde{\boldsymbol{\beta}}_{T}$ and the
initial values $\tilde{\beta}_{1},\ldots,\tilde{\beta}_{d}$. In this case,
${\mathbf{\boldsymbol{\vartheta}}}=(\theta_{1},\ldots,\theta_{d})$ and the
Kalman filter described above is applied with
$\hat{\boldsymbol{\beta}}_{0|0}={\mathbf{0}}_{d}$ and
${\mathbf{P}}_{0|0}=\mbox{\rm
Diag}\left(A_{0,11}+\theta_{1}P_{0,{11}},\ldots,A_{0,dd}+\theta_{d}P_{0,{dd}}\right)$.
The Kalman filter reads:

$\displaystyle\tilde{\boldsymbol{\beta}}_{t}|\delta,{\mathbf{\boldsymbol{\vartheta}}},\tilde{\boldsymbol{z}}^{t-1}\sim\mathcal{N}_{d}\left(\hat{\boldsymbol{\beta}}_{t-1|t-1},\delta{\mathbf{P}}_{t|t-1}\right),\quad{\mathbf{P}}_{t|t-1}={\mathbf{P}}_{t-1|t-1}+\mbox{\rm
Diag}\left(\theta_{1},\ldots,\theta_{d}\right),$ (A.48)
$\displaystyle\tilde{z}_{t}|\delta,{\mathbf{\boldsymbol{\vartheta}}},\tilde{\boldsymbol{z}}^{t-1}\sim\mathcal{N}\left(\hat{z}_{t|t-1},\delta
S_{t|t-1}\right),\quad\hat{z}_{t|t-1}={\mathbf{x}}_{t}\hat{\boldsymbol{\beta}}_{t-1|t-1},\quad
S_{t|t-1}={\mathbf{x}}_{t}{\mathbf{P}}_{t|t-1}{\mathbf{x}}_{t}^{\top}+{\omega_{t}},$
$\displaystyle\tilde{\boldsymbol{\beta}}_{t}|\delta,{\mathbf{\boldsymbol{\vartheta}}},\tilde{\boldsymbol{z}}^{t}\sim\mathcal{N}_{d}\left(\hat{\boldsymbol{\beta}}_{t|t},\delta{\mathbf{P}}_{t|t}\right),$
$\displaystyle\hat{\boldsymbol{\beta}}_{t|t}=\hat{\boldsymbol{\beta}}_{t-1|t-1}+{\mathbf{K}}_{t}(\tilde{z}_{t}-\hat{z}_{t|t-1}),\quad{\mathbf{P}}_{t|t}=({\mathbf{I}}-{\mathbf{K}}_{t}{\mathbf{x}}_{t}){\mathbf{P}}_{t|t-1},,\quad{\mathbf{K}}_{t}={{\mathbf{P}}_{t|t-1}}{\mathbf{x}}_{t}^{\top}S_{t|t-1}^{-1}.$

Therefore, the likelihood (A.46) takes the following form:

$\displaystyle
p(\tilde{\boldsymbol{z}}|\delta,{\mathbf{\boldsymbol{\vartheta}}})=\prod_{t=1}^{T}p(\tilde{z}_{t}|\tilde{\boldsymbol{z}}^{t-1},\delta,{\mathbf{\boldsymbol{\vartheta}}})\propto\left(\frac{1}{\delta}\right)^{\frac{T}{2}}\exp\left\\{-\frac{1}{2\delta}\sum_{t=1}^{T}\frac{(\tilde{z}_{t}-\hat{z}_{t|t-1})^{2}}{S_{t|t-1}}\right\\}.$

###  A.5 Additional simulation results

This section provides some additional simulation results for probit, logit,
multinomial and binomial regression models. We investigate a number of
simulation scenarios. First, we vary the sample size $N$ with only one
instance of $y_{i}=1$, implying increasing imbalance with a larger sample
size. Second, we fix $N=1000$ and vary the intercept in the data generating
process, where large absolute intercepts imply greater imbalance than
intercepts that are close to zero (corresponding to the fully balanced case).
In addition, we vary the number of choices $\\{0,\ldots,m\\}$ in the
multinomial case and the number of successes $N_{i}$ in the binomial case.

The results for the probit model are discussed in detail in Section A.5.1.
Graphical summaries of the additional simulation results of the UPG models for
the binary logit model, the multinomial model for the binomial logit model are
provided in Section A.5.2. Put briefly, the additional simulation results
suggest that the proposed UPG approach can compete with the efficiency ”gold
standard” of the Metropolis-Hastings algorithm. At the same time, several
extensions of the UPG approach are easily implemented due to the conditional
Gaussian posterior of the regression coefficients.

(a) Varying sample sizes with one success.

(b) Varying intercepts with N = 1000.

Figure A.2: Inefficiency factors across 50 simulation runs for binary probit
models.

####  A.5.1 Probit results

For binary probit models, we compare the marginal data augmentation approach
from Imai & van Dyk (2005) (I&vD),111111We use the UPG algorithm with $\gamma$
deterministically fixed at $0$ which corresponds to the marginal data
augmentation algorithm in Imai & van Dyk (2005).

the UPG approach outlined in this article, the adaptive Metropolis-Hastings
step proposed in Johndrow et al. (2019) (AMH) as well as the original probit
data augmentation scheme outlined in Albert & Chib (1993) (A&C).

Figure A.3: Inefficiency factors varying sample size with one success in
binary logit models.

The results for the probit model are summarized graphically in Figure A.2 and
are provided numerically in Table A.1. It is obvious that the commonly used
standard data augmentation approach from A&C has severe efficiency problems
when increasing sparsity in the dependent variable and therefore
imbalancedness of the data. Increasing inefficiency can also be observed for
the I&vD algorithm as seen in Figure A.2 and Table A.1. On the other hand,
constantly high levels of efficiency are observed for UPG and the AMH
algorithm. In the probit case, the extreme intercept values of $-3$ and $3$
imply success counts of close to $0$ and close to $N$, respectively. In
general, AMH and UPG exhibit efficient estimation behavior over the full range
of intercepts. On the contrary, I&vD as well as A&C show increasing
inefficiency with increasing imbalance.

Figure A.4: Inefficiency factors for varying number of categories $K=m+1$ in
multinomial logit models with one success and $N=1000$.

####  A.5.2 Additional results

Additional results are provided in Figure A.3 for the binary logit model, in
Figure A.5 and Figure A.4 for the multinomial model and in Figure A.6 for the
binomial logit model. Selected numerical results for these models can be found
in Table A.2 to Table A.5.

Figure A.5: Inefficiency factors for varying sample size in multinomial logit
models with one success and three categories. Figure A.6: Varying number of
successes $\text{N}_{i}$ in binomial logit models with intercept
$\beta_{d}=-3$ and $N=500$.

####  A.5.3 Tabulated simulation results

Probit N = 10 N = 100 N = 1000 N = 10000 N = 20000 N = 50000 I&vD
$\text{IF}^{\text{B}}$ $1.8676$ $4.8926$ $21.1759$ $109.6955$ $182.3035$