# Conjugate priors for count and rounded data regression

Daniel R. Kowal Dobelman Family Assistant Professor, Department of Statistics,
Rice University<EMAIL_ADDRESS>
###### Abstract

Discrete data are abundant and often arise as counts or rounded data. However,
even for linear regression models, conjugate priors and closed-form posteriors
are typically unavailable, thereby necessitating approximations or Markov
chain Monte Carlo for posterior inference. For a broad class of count and
rounded data regression models, we introduce conjugate priors that enable
closed-form posterior inference. Key posterior and predictive functionals are
computable analytically or via direct Monte Carlo simulation. Crucially, the
predictive distributions are discrete to match the support of the data and can
be evaluated or simulated jointly across multiple covariate values. These
tools are broadly useful for linear regression, nonlinear models via basis
expansions, and model and variable selection. Multiple simulation studies
demonstrate significant advantages in computing, predictive modeling, and
selection relative to existing alternatives.

KEYWORDS: Bayesian statistics; discrete data; generalized linear models;
integer data; prediction.

##  1 Introduction

Discrete data commonly arise either as natural counts or via data coarsening
such as rounding. Bayesian regression models must account for this
discreteness in order to provide a coherent data-generating process. However,
discrete data often present challenging distributional features, including
zero-inflation, over/under-dispersion, boundedness or censoring, and heaping.
Customized Bayesian models that target these features are usually too complex
for analytical posterior and predictive analysis and can incur significant
computational burdens from the requisite approximation algorithms, such as
Markov chain Monte Carlo. The contribution of this paper is to provide
conjugate priors, closed-form posteriors, and efficient Monte Carlo posterior
and predictive simulation algorithms for a broad class of discrete data
regression models.

The impact of data coarsening is often nonignorable, both for coherency of the
data-generating process and for correct and efficient inference. Dempster and
Rubin, (1983) emphasized the surprisingly large impact of rounded data
adjustments in least squares regression, which notably do not decay with the
sample size. More broadly, Heitjan, (1989) reviewed inference and asymptotics
for grouped data, where continuous data are observable only up to some
coarsening such as an interval. Heitjan and Rubin, (1991) studied multiple
sources of data coarsening, including rounding, heaping, and, censoring, as
well as random coarsening operators.

The perspective on discrete data offered by grouping or rounding has also been
popularized as a modeling strategy for count and binary data. Canale and
Dunson, (2011) introduced rounded Gaussian mixture models for count-valued
Bayesian nonparametrics, while Canale and Dunson, (2013) similarly applied
rounded Gaussian processes for count regression modeling. Kowal and Canale,
(2020) developed extensions that incorporated both rounding and transformation
operators to adapt linear, additive, and BART regression models for count and
rounded data. Although these works demonstrate the utility of discretized
continuous data models for count data, none of these methods provides analytic
posterior or predictive distributions, and instead use Markov chain Monte
Carlo for inference. Among related non-Bayesian approaches, Siegfried and
Hothorn, (2020) and Kowal and Wu, (2021) developed count data regression with
an emphasis on learning a transformation for more adaptive discrete
distributional modeling. A similar strategy was deployed by Jia et al., (2021)
for count time series analysis.

Conjugate priors for discrete data regression have received renewed attention,
including for binary (Durante,, 2019) and multinomial (Fasano and Durante,,
2020) data. Indeed, the closed-form posteriors for probit linear regression
derived by Durante, (2019) are a special case of our results. However, there
are additional challenges in our setting. First, our results generalize from
binary data $y\in\\{0,1\\}$ to discrete data $y\in\mathbb{Z}$, including count
data, rounded data, and bounded or censored (discrete) data. These
generalizations require a different class of distributions for conjugate
modeling, and therefore distinct descriptions and computations of the relevant
posterior and predictive quantities. Second, our framework incorporates a
transformation, which is irrelevant for binary data but widely-used in count
data models. For instance, count data are often modeled (incoherently) using
continuous regression models after log or square-root transformations, while
copula models incorporate transformations via the marginal distributions (Pitt
et al.,, 2006). Third, we prioritize closed-form predictive distributions that
are discrete to match the support of the data and joint across multiple
covariate values. Naturally, a predictive distribution supported on (a subset
of) $\mathbb{Z}$ is a substantial generalization of a binary predictive
distribution, which is completely characterized by the success probability
$\mbox{pr}(\tilde{y}=1\mid y)$. Our results admit direct predictive
probability evaluations and Monte Carlo predictive sampling for discrete data
$y\in\mathbb{Z}$, including for sparse regression, nonlinear models, and
model-averaged scenarios.

Bayesian regression models for count data commonly build upon the Poisson
distribution. Except in special cases, the posterior and predictive
distributions do not have known closed forms and therefore require
approximations such as Markov chain Monte Carlo. Further, the Poisson
distribution is often inadequate for modeling, especially in the case of
over/under-dispersion, zero-inflation, or heaping. Generalizations such as the
Conway-Maxwell-Poisson distribution (Sellers and Shmueli,, 2010) or the
discrete Weibull distribution (Peluso et al.,, 2019) provide more flexible
count distribution models, yet are burdened by more intensive computations.
Polson et al., (2013) introduced a data augmentation strategy for Gibbs
sampling with Negative Binomial likelihoods, while Bradley et al., (2018)
proposed log-gamma distributions that are conditionally conjugate to the
Poisson distribution. These approaches employ Markov chain Monte Carlo for all
posterior and predictive inference and typically require further modifications
to account for zero-inflation, bounded or censored data, and heaping.

##  2 Conjugate priors for count and rounded data

###  2.1 Transformation and rounding models

For paired data $\\{(x_{i},y_{i})\\}_{i=1}^{n}$ with covariates
$x_{i}\in\mathbb{R}^{p}$ and integer-valued responses $y_{i}\in\mathbb{Z}$,
consider the following class of models:

$\displaystyle y_{i}$ $\displaystyle=h\circ g^{-1}(z_{i})$ (1) $\displaystyle
z_{i}$ $\displaystyle\sim\Pi_{\theta}(\cdot\mid x_{i})$ (2)

where $h:\mathbb{R}\to\mathbb{Z}$ is a _rounding operator_ ,
$g:\mathbb{R}\to\mathbb{R}$ is a (monotone increasing) _transformation_ , and
$\Pi_{\theta}(\cdot\mid x_{i})$ is a covariate-dependent _continuous
distribution_ parametrized by $\theta$. Informally, (1)–(2) is designed to
adapt a continuous data model $\Pi_{\theta}$ for coherent modeling of discrete
data $\\{y_{i}\\}_{i=1}^{n}$. We focus primarily on the Gaussian linear model
for the latent continuous data $z_{i}$,

$z_{i}=x_{i}^{\prime}\theta+\epsilon_{i},\quad\epsilon_{i}\sim
N(0,\sigma^{2})$ (3)

independently for $i=1,\ldots,n$, with modifications available for
heteroskedastic or dependent errors and nonlinear versions (see Section 2.5).

It is useful to decouple the rounding operator $h$, which defines the support
of $y$, and the transformation $g$, which is a smooth monotone function.
First, consider the rounding operator $h$, which is defined via a partition
$\\{\mathcal{A}_{j}\\}$ of $\mathbb{Z}$ such that $y_{i}=j$ whenever $z_{i}\in
g(\mathcal{A}_{j})$.

######  Example 1.

Suppose that we observe rounded data instead of continuous data. Let the
rounded digit be the first decimal place, which can be satisfied by pre-
multiplying the data by the appropriate power of 10 such that
$y\in\mathbb{Z}$. When $g(t)=t$ is the identity, $z_{i}$ represents the
unobserved continuous data version of $y_{i}$. In this case, $y_{i}=j$
whenever $z_{i}\in\mathcal{A}_{j}=[j-0.5,j+0.5)$, so $h(t)=\lfloor
t+0.5\rfloor$ and $\lfloor\cdot\rfloor$ is the floor function.

######  Example 2.

Suppose that $y\in\\{0,\ldots,y_{max}\\}$ is count-valued and bounded by some
known $y_{max}$. The rounding operator defined by $h(t)=0$ for $t<1$,
$h(t)=\lfloor t\rfloor$ for $t\in[1,y_{max}]$, and $h(t)=y_{max}$ for
$t>y_{max}$ ensures the correct support for the data-generating process.
Equivalently, the partition is $\mathcal{A}_{0}=(-\infty,1)$,
$\mathcal{A}_{j}=[j,j+1)$ for $j=1,\ldots,y_{max}-1$, and
$A_{y_{max}}=[y_{max},\infty)$. In addition, (1)–(2) coupled with this choice
of $h$ produces a valid likelihood for censored data (Kowal and Canale,,
2020). The case of unbounded counts simply requires $y_{max}=\infty$. When
$y_{max}=1$, our results reproduce Durante, (2019) for probit regression.

Next, consider the monotone increasing transformation $g$.

######  Example 3.

The (signed) Box-Cox family (Box and Cox,, 1964) provides a parametric option,
$g(t;\lambda)=\\{\mbox{sgn}(t)|t|^{\lambda}\\}/\lambda$, for $\lambda>0$ with
$g(t;\lambda=0)=\log(t)$ for $t>0$. In addition to the log transformation,
this family includes the (shifted and scaled) identity ($\lambda=1$) and
square-root $(\lambda=1/2)$ transformations. These transformations are common
for (incoherent) application of continuous models (2) to count data without
the rounding operator $h$; for example, the square-root transformation is the
variance-stabilizing transformation of the Poisson distribution. Kowal and
Canale, (2020) used this family for Bayesian modeling of (1)–(2).

######  Example 4.

Using copula-style models, an alternative strategy is to link the
transformation to the cumulative distribution function of $y$ (Liu et al.,,
2009; Siegfried and Hothorn,, 2020; Kowal and Wu,, 2021). Suppose
$\mathcal{A}_{j}=[a_{j},a_{j+1})$ as in the examples above. The cumulative
distribution functions of $y$ and $z$ are related as follows:

$F_{y}(j)=\mbox{pr}(y\leq
j)=\mbox{pr}\\{z<g(a_{j+1})\\}=F_{z}\\{g(a_{j+1})\\}.$ (4)

By model (2), the latent $z$ are Gaussian. Consider a working marginal
distribution $z\sim N(\mu_{z},\sigma_{z}^{2})$, so
$F_{z}(t)=\Phi\\{(t-\mu_{z})/\sigma_{z}\\}$ for $t\in\mathbb{R}$. Rearranging
(4) motivates the following estimator:

$\hat{g}_{0}(a_{j+1})=\bar{y}+\hat{s}_{y}\Phi^{-1}\\{\tilde{F}_{y}(j)\\}$ (5)

along with the appropriate limiting behavior, i.e.,
$\hat{g}_{0}(-\infty)=-\infty$ and $\hat{g}_{0}(\infty)=\infty$. Since the
transformation is identifiable only up to location ($\mu_{z}$) and scale
($\sigma_{z}$), the estimator (5) matches the first two moments of $z$ with
those of $y$: $\mu_{z}=\bar{y}$ is the sample mean and
$\sigma_{z}=\hat{s}_{y}$ is the sample standard deviation of
$\\{y_{i}\\}_{i=1}^{n}$. Similarly, (5) substitutes a point estimate for
$F_{y}$ using the nonparametric estimator
$\tilde{F}_{y}(j)=n(n+1)^{-1}\hat{F}_{y}(j)$ that rescales the empirical
cumulative distribution function
$\hat{F}_{y}(j)=n^{-1}\sum_{i=1}^{n}\mathbb{I}\\{y_{i}\leq j\\}$ to avoid
boundary issues. Lastly, Kowal and Wu, (2021) proposed a smooth interpolation
of $\hat{g}_{0}$ on the observed values,
$\hat{g}({y_{i}})=\hat{g}_{0}({y_{i}})$ for $i=1,\ldots,n$, where where
$\hat{g}$ is a monotonic cubic interpolating spline (Fritsch and Carlson,,
1980). The limits of $\hat{g}_{0}$ are also imposed.

The distinction between $\hat{g}_{0}$ and $\hat{g}$ highlights the importance
of the rounding operator: $\hat{g}_{0}$ is a step function with
$h\circ\hat{g}_{0}^{-1}=\hat{g}_{0}^{-1}$ and therefore does not need the
rounding operator. However, the data-generating process of (1)–(2) with
transformation $\hat{g}_{0}$ is supported only on the observed values:
$\mbox{pr}(y=j)=0$ whenever $g(a_{j})=g(a_{j+1})$, which occurs for
$\hat{g}_{0}$ when $j\not\in\\{y_{i}\\}_{i=1}^{n}$. Instead, the smoothness of
$\hat{g}$ appropriately expands the support of the data-generating process,
while the interpolation ensures that $\hat{g}$ preserves (5) at the observed
data values $\\{y_{i}\\}_{i=1}^{n}$.

######  Example 5.

The estimator (5) may instead use parametric distributions, such as Poisson or
Negative Binomial distribution functions. In this case, the role of the
transformation is to map the latent data model (2) such that the marginal
distribution of $y$ conforms to a pre-specified parametric family. The
accompanying parameters of $F_{y}$ must be estimated; a useful default is
method-of-moments estimation, which recycles $\bar{y}$ and $\hat{s}_{y}$. Jia
et al., (2021) adopted a related approach for count time series modeling,
while Pitt et al., (2006) similarly used parametric regression models for the
marginals in a multivariate Gaussian copula model.

######  Example 6.

Semiparametric copula models similarly apply parametric or nonparametric
estimates of $F_{y}$, but then define an extended rank (probit) likelihood
that only uses the ordering of the data (Hoff,, 2007; Feldman and Kowal,,
2021). In particular, the ordering of $\\{y_{i}\\}$ imposes an equivalent
ordering on $\\{z_{i}\\}$, which ensures that conditioning on the data is
equivalent to applying a constraint region on $z_{i}$. This result is closely
linked to the framework in the subsequent section.

###  2.2 Selection distributions

Under model (1)–(2), the likelihood for $\theta$ is

$\mbox{pr}(y_{i}=j\mid\theta)=\mbox{pr}\\{z_{i}\in
g(\mathcal{A}_{j})\mid\theta\\}=\Pi_{\theta}\\{g(\mathcal{A}_{j})\mid
x_{i}\\}$ (6)

with conditional independence across $i=1,\ldots,n$. Given a prior
$p(\theta)$, the posterior distribution is

$p(\theta\mid y)=\frac{p(\theta)\mbox{pr}\\{z\in
g(\mathcal{A}_{y})\mid\theta\\}}{\mbox{pr}\\{z\in
g(\mathcal{A}_{y})\\}}=\frac{p(\theta)\Pi_{\theta}\\{g(\mathcal{A}_{y})\mid
x\\}}{\int p(\theta)\Pi_{\theta}\\{g(\mathcal{A}_{y})\mid x\\}d\theta}$ (7)

where $\\{z\in g(\mathcal{A}_{y})\\}=\\{z_{1}\in
g(\mathcal{A}_{y_{1}}),\ldots,z_{n}\in g(\mathcal{A}_{y_{n}})\\}$ is
elementwise.

Previous attempts at Bayesian modeling of (1)–(2) have applied Gibbs sampling
algorithms for estimating functionals of (7) (Canale and Dunson,, 2011, 2013;
Kowal and Canale,, 2020). The common strategy is to introduce a data-
augmentation step in the Markov chain Monte Carlo algorithm and iteratively
sample from $[z\mid y,\theta]$ and $[\theta\mid z,y]=[\theta\mid z]$. However,
these Markov chain Monte Carlo algorithms fail to provide general insights
about the posterior or predictive distribution under (1)–(2) and require
customized diagnostics and lengthy simulation runs.

Here, we characterize the posterior distribution explicitly:

######  Theorem 1.

When $\mathcal{C}=g(\mathcal{A}_{y})$ is a measurable subset of of
$\mathbb{R}^{n}$, the posterior distribution $p(\theta\mid y)$ under model
(1)–(2) is a _selection distribution_.

The result is self-evident from (6)–(7): a selection distribution is defined
by random variables of the form $[\theta\mid z\in\mathcal{C}]$ and
parametrized in terms of the joint distribution for $(\theta,z)$ and the
measurable constraint region $\mathcal{C}$. Hence, the observation that
$p(\theta\mid y)=p\\{\theta\mid z\in g(\mathcal{A}_{y})\\}$ is sufficient.

General properties of selection distributions are provided by Arellano-Valle
et al., (2006). When the prior $p(\theta)$ is closed under a set of
transformations or closed under marginalization, these properties propagate to
the posterior $p(\theta\mid y)$. When $(\theta,z)$ have a joint multivariate
elliptically contoured distribution, the posterior density and distribution
functions can be expressed in terms of the location parameters, dispersion
parameters, and density generator function; see Arellano-Valle et al., (2006)
for the explicit forms.

For a more tractable analysis, suppose that the prior $p(\theta)$ and the
latent data model $\Pi_{\theta}(\cdot\mid x_{i})$ are both Gaussian.
Specifically, construct the joint distribution

$\begin{pmatrix}z\\\ \theta\end{pmatrix}\sim
N_{n+p}\left\\{\begin{pmatrix}\mu_{z}\\\
\mu_{\theta}\end{pmatrix},\begin{pmatrix}\Sigma_{z}&\Sigma_{z\theta}\\\
\Sigma_{z\theta}^{\prime}&\Sigma_{\theta}\end{pmatrix}\right\\}$ (8)

for $z\in\mathbb{R}^{n}$ and $\theta\in\mathbb{R}^{p}$, where the moments may
depend on the covariates $\\{x_{i}\\}_{i=1}^{n}$. Then the posterior may be
characterized as follows:

######  Theorem 2.

Under the model (1) and (8), the posterior distribution is _selection normal_
, denoted $[\theta\mid
y]\sim\mbox{SLCT-N}_{n,p}(\mu_{z},\mu_{\theta},\Sigma_{z},\Sigma_{\theta},\Sigma_{z\theta},\mathcal{C})$,
with constraint region $\mathcal{C}=g(\mathcal{A}_{y})$, density

$p(\theta\mid
y)=\phi_{p}(\theta;\mu_{\theta},\Sigma_{\theta})\frac{\bar{\Phi}_{n}\\{\mathcal{C};\Sigma_{z\theta}\Sigma_{\theta}^{-1}(\theta-\mu_{\theta})+\mu_{z},\Sigma_{z}-\Sigma_{z\theta}\Sigma_{\theta}^{-1}\Sigma_{z\theta}^{\prime}\\}}{\bar{\Phi}_{n}(\mathcal{C};\mu_{z},\Sigma_{z})},$
(9)

and moment generating function $M_{[\theta\mid
z\in\mathcal{C}]}(s)=\exp\big{(}s^{\prime}\mu_{\theta}+s^{\prime}\Sigma_{\theta}s/2\big{)}\bar{\Phi}_{n}(\mathcal{C};\Sigma_{z\theta}s+\mu_{z},\Sigma_{z})/\bar{\Phi}_{n}(\mathcal{C};\mu_{z},\Sigma_{z})$,
where $\phi_{p}(\cdot;\mu,\Sigma)$ denotes the Gaussian density function of a
Gaussian random variable with mean $\mu$ and covariance $\Sigma$ and
$\bar{\Phi}_{n}(\mathcal{C};\mu,\Sigma)=\int_{\mathcal{C}}\phi_{n}(x;\mu,\Sigma)\
dx$.

The posterior distribution is defined by the joint moments of $(\theta,z)$ in
(8) and the constraint region implied by (1). Crucially, there exists a more
constructive representation:

######  Theorem 3.

Let $V_{0}\sim N_{n}(0,\Sigma_{z})$ and $V_{1}\sim
N_{p}(0,\Sigma_{\theta}-\Sigma_{z\theta}^{\prime}\Sigma_{z}^{-1}\Sigma_{z\theta})$
be independent with $V_{0}^{(\mathcal{C}-\mu_{z})}=[V_{0}\mid
V_{0}\in\mathcal{C}-\mu_{z}]$ for $\mathcal{C}=g(\mathcal{A}_{y})$. Then
$[\theta\mid
y]=\mu_{\theta}+V_{1}+\Sigma_{z\theta}^{\prime}\Sigma_{z}^{-1}V_{0}^{(\mathcal{C}-\mu_{z})}$.

Most important, Theorem 3 enables direct Monte Carlo simulation from the
posterior, rather than iterative Markov chain Monte Carlo. Specific examples
are provided subsequently.

While simplifications are often available depending on the model form (see
Sections 2.3 \- 2.5), the primary computational bottleneck for evaluating the
density (9) is the integration of an $n$-dimensional multivariate normal
density, whereas the Monte Carlo sampler in Theorem 3 is limited by the draw
from the $n$-dimensional truncated normal. Our computing strategies leverage
recent algorithmic developments. In particular, Botev, (2017) introduced a
minimax exponential tilting strategy that accurately estimates the normalizing
constant of a high-dimensional multivariate normal distribution and provides
an efficient importance sampling algorithm for truncated multivariate normal
variables with an accept-reject algorithm that achieves a high acceptance
rate. These algorithms are efficient and accurate for moderate $n$, such as
$n\approx 500$, and are implemented in the TruncatedNormal package in R.

As a secondary consequence of Theorem 3, the posterior mean and posterior
variance of the regression coefficients are readily available:

$\displaystyle E(\theta\mid y)$
$\displaystyle=\mu_{\theta}+\Sigma_{z\theta}^{\prime}\Sigma_{z}^{-1}E(V_{0}\mid
V_{0}\in\mathcal{C}-\mu_{z})$ (10) $\displaystyle\mbox{var}(\theta\mid y)$
$\displaystyle=\Sigma_{\theta}+\Sigma_{z\theta}^{\prime}\Sigma_{z}^{-1}\\{\Sigma_{z}-\mbox{var}(V_{0}\mid
V_{0}\in\mathcal{C}-\mu_{z})\\}\Sigma_{z}^{-1}\Sigma_{z\theta}$ (11)

which depend on the moments of $V_{0}^{(\mathcal{C}-\mu_{z})}=[V_{0}\mid
V_{0}\in\mathcal{C}-\mu_{z}]$. These and other terms can be simplified
depending on the model (8).

###  2.3 Posterior distribution for the linear model

Consider the linear model (3) with a Gaussian prior $\theta\sim
N(\mu_{\theta},\Sigma_{\theta})$ and let $X=(x_{1},\ldots,x_{n})^{\prime}$
denote the $n\times p$ covariate matrix. The posterior distribution of the
regression coefficients is available as a direct consequence of Theorem 2:

######  Lemma 1.

For model (1) and (3) with $\theta\sim N(\mu_{\theta},\Sigma_{\theta})$, the
posterior distribution is $[\theta\mid
y]\sim\mbox{SLCT-N}_{n,p}(\mu_{z}=X\mu_{\theta},\mu_{\theta},\Sigma_{z}=X\Sigma_{\theta}X^{\prime}+\sigma^{2}I_{n},\Sigma_{\theta},\Sigma_{z\theta}=X\Sigma_{\theta},\mathcal{C}=g(\mathcal{A}_{y}))$.

Lemma 1 enables evaluation of the posterior density (Theorem 2), direct Monte
Carlo simulation from the posterior (Theorem 3), and computation of the
posterior moments (equations (10) and (11)) for the linear model.

To obtain more analytically tractable posterior quantities, consider the
$g$-prior (Zellner,, 1986) for $\theta$. Monte Carlo sampling from the
posterior is simplified as follows:

######  Lemma 2.

Suppose $\Sigma_{\theta}=\psi\sigma^{2}(X^{\prime}X)^{-1}$ for $\psi>0$. Under
model (1) and (3), a posterior draw $\theta^{*}\sim p(\theta\mid y)$ is
obtained as follows:

  1. 1.

Simulate $V_{0}^{*}\sim N_{n}(0,\sigma^{2}\\{\psi
X(X^{\prime}X)^{-1}X^{\prime}+I_{n}\\})$ truncated to
$g(\mathcal{A}_{y})-X\mu_{\theta}$.

  2. 2.

Simulate $V_{1}^{*}\sim
N_{p}(0,\sigma^{2}\psi(1+\psi)^{-1}(X^{\prime}X)^{-1})$.

  3. 3.

Set
$\theta^{*}=\mu_{\theta}+V_{1}^{*}+\psi(1+\psi)^{-1}(X^{\prime}X)^{-1}X^{\prime}V_{0}^{*}$.

Lemma 2 provides direct simulation from the posterior distribution of $\theta$
using standard regression functionals of $X$. Lemma 2 also simplifies the
posterior expectation of the regression coefficients:

$E(\theta\mid
y)=\mu_{\theta}+\frac{\psi}{1+\psi}(X^{\prime}X)^{-1}X^{\prime}\hat{z}$ (12)

where $\hat{z}=E\\{z\mid z\in
g(\mathcal{A}_{y})-X\mu_{\theta}\\}=E(V_{0}^{*})$ is the marginal expectation
of $z$ truncated to $g(\mathcal{A}_{y})-X\mu_{\theta}$. Using Lemma 2,
$E(\theta\mid y)$ is easily estimable by replacing $\hat{z}$ with the sample
mean of Monte Carlo draws of $V_{0}^{*}$. When the prior is centered at zero,
$\mu_{\theta}=0$ and $\hat{z}=E\\{z\mid z\in g(\mathcal{A}_{y}\\}=E(z\mid y)$
is the posterior expectation of the latent data $z$ unconditional on $\theta$.
By comparison, the zero-centered $g$-prior with a Gaussian likelihood for $y$
produces a posterior expectation of the same form, but with the observed data
$y$ in place of $\hat{z}$. This comparison is particularly insightful in the
context of the rounded data scenario from Example 1. In this case, $z_{i}$ is
the latent continuous version of the rounded data $y_{i}$, so
$\hat{z}_{i}=E(z_{i}\mid y)$ is a reasonable proxy for the unobservable
continuous data in the posterior expectation (12).

More broadly, Lemma 1 suggests a conjugate prior for model (1) and (3). In
particular, the Gaussian distribution is a special case of the selection
normal:
$\theta\sim\mbox{SLCT-N}_{1,p}(\mu_{z}=0,\mu_{\theta},\Sigma_{z}=1,\Sigma_{\theta},\Sigma_{z\theta}=0_{p}^{\prime},\mathcal{C}=\mathbb{R})$,
where the moments and constraints on $z$ are irrelevant as long as
$\Sigma_{z\theta}=0$. By generalizing these additional arguments, we obtain a
richer class of conjugate priors:

######  Lemma 3.

Under the prior
$\theta\sim\mbox{SLCT-N}_{n_{0},p}(\mu_{z_{0}},\mu_{\theta},\Sigma_{z_{0}},\Sigma_{\theta},\Sigma_{z_{0}\theta},\mathcal{C}_{0})$
and the model (1) and (3), the posterior is $[\theta\mid
y]\sim\mbox{SLCT-N}_{n_{0}+n,p}(\mu_{z_{1}},\mu_{\theta},\Sigma_{z_{1}},\Sigma_{\theta},\Sigma_{z_{1}\theta},\mathcal{C}_{1})$
with

$\mu_{z_{1}}=\begin{pmatrix}\mu_{z_{0}}\\\ X\mu_{\theta}\end{pmatrix},\
\Sigma_{z_{1}}=\begin{pmatrix}\Sigma_{z_{0}}&\Sigma_{z_{0}\theta}X^{\prime}\\\
X\Sigma_{z_{0}\theta}^{\prime}&X\Sigma_{\theta}X^{\prime}+\sigma^{2}I_{n}\end{pmatrix},\
\Sigma_{z_{1}\theta}=\begin{pmatrix}\Sigma_{z_{0}\theta}\\\
X\Sigma_{\theta}\end{pmatrix},\ \mathcal{C}_{1}=\mathcal{C}_{0}\times
g(\mathcal{A}_{y}).$

Perhaps the most useful setting of Lemma 3 is sequential updating based on
multiple datasets, say $y^{(1)}$ and $y^{(2)}$. When these datasets are
conditionally independent and follow the same model (1) and (3), the total
posterior decomposes as $p(\theta\mid y^{(1)},y^{(2)})\propto p(\theta\mid
y^{(1)})p(y^{(2)}\mid\theta)$. Under a Gaussian prior for $p(\theta)$, Lemma 1
implies that the partial posterior $p(\theta\mid y^{(1)})$ is selection
normal, which serves as the prior for the second dataset update. Hence, Lemma
3 provides the direct updates from the partial posterior $p(\theta\mid
y^{(1)})$ to the total posterior $p(\theta\mid y^{(1)},y^{(2)})$.

###  2.4 Prediction

A crucial feature of the model (1)–(2) is that the data-generating process and
therefore the predictive distribution is discrete and matches the support of
$y$. Consider the posterior predictive distribution at the $\tilde{n}\times p$
covariate matrix
$\tilde{X}=(\tilde{x}_{1},\ldots,\tilde{x}_{\tilde{n}})^{\prime}$. Because of
the link in (1), the joint predictive distribution of
$\tilde{y}=(\tilde{y}_{1},\ldots,\tilde{y}_{\tilde{n}})^{\prime}$ can be
expressed via the latent variables $\tilde{z}$:

$\mbox{pr}(\tilde{y}\mid y)=\mbox{pr}\\{\tilde{z}\in
g(\mathcal{A}_{\tilde{y}})\mid y\\}=\int\mbox{pr}\\{\tilde{z}\in
g(\mathcal{A}_{\tilde{y}})\mid\theta\\}p(\theta\mid y)d\theta.$ (13)

We consider multiple computing strategies for (functionals of) the posterior
predictive distribution. The simplest option is to leverage the Monte Carlo
sampling algorithms for $p(\theta\mid y)$ (e.g., Theorem 3) to obtain latent
predictive samples from $[\tilde{z}\mid y]$, and then map those to $\tilde{y}$
via (1):

######  Lemma 4.

A predictive draw $\tilde{y}^{*}\sim p(\tilde{y}\mid y)$ for model (1)–(2) is
generated as follows: sample $\theta^{*}\sim p(\theta\mid y)$; sample
$\tilde{z}^{*}\sim\Pi_{\theta=\theta^{*}}(\cdot\mid\tilde{X})$; set
$\tilde{y}^{*}=h\circ g^{-1}(\tilde{z}^{*})$.

The result is a direct consequence of the conditional independence assumptions
in (2) and the deterministic link in (1). The samples are discrete to match
the support of $y$ and capture the joint dependence among predictive variables
at multiple covariates $\tilde{X}$. For the linear model (3), the posterior
sampling step is available via Monte Carlo sampling (e.g., Lemma 2). As a
result, Lemma 4 produces Monte Carlo samples from the joint posterior
predictive distribution. The latent predictive step is simply
$\tilde{z}^{*}\sim
N_{\tilde{n}}(\tilde{X}\theta^{*},\sigma^{2}I_{\tilde{n}})$, so Lemma 4
requires minimal additional computational cost given $\theta^{*}\sim
p(\theta\mid y)$.

An alternative Monte Carlo sampler for the linear model circumvents the
posterior sampling and instead draws directly from the latent posterior
predictive distribution of $[\tilde{z}\mid y]$, which can be derived
explicitly:

######  Lemma 5.

For the linear model (3) with a Gaussian prior $\theta\sim
N(\mu_{\theta},\Sigma_{\theta})$, the latent predictive distribution at
$\tilde{X}$ is $[\tilde{z}\mid
y]\sim\mbox{SLCT-N}_{n,\tilde{n}}(X\mu_{\theta},\tilde{X}\mu_{\theta},X\Sigma_{\theta}X^{\prime}+\sigma^{2}I_{n},\tilde{X}\Sigma_{\theta}\tilde{X}^{\prime}+\sigma^{2}I_{\tilde{n}},X\Sigma_{\theta}\tilde{X}^{\prime},\mathcal{C}=g(\mathcal{A}_{y}))$.

Direct Monte Carlo simulation proceeds using Theorem 3 applied to the latent
predictive selection normal distribution, which we simplify for the $g$-prior:

######  Lemma 6.

Suppose $\Sigma_{\theta}=\psi\sigma^{2}(X^{\prime}X)^{-1}$ for $\psi>0$. Under
model (1) and (3), a joint predictive draw $\tilde{y}^{*}\sim p(\tilde{y}\mid
y)$ at the $\tilde{n}\times p$ covariate matrix $\tilde{X}$ is obtained as
follows:

  1. 1.

Simulate $V_{0}^{*}\sim N_{n}(0,\sigma^{2}\\{\psi
X(X^{\prime}X)^{-1}X^{\prime}+I_{n}\\})$ truncated to
$g(\mathcal{A}_{y})-X\mu_{\theta}$.

  2. 2.

Simulate $\tilde{V}_{1}^{*}\sim
N_{\tilde{n}}(0,\sigma^{2}\\{\psi(1+\psi)^{-1}\tilde{X}(X^{\prime}X)^{-1}\tilde{X}^{\prime}+I_{\tilde{n}}\\})$.

  3. 3.

Compute
$\tilde{z}^{*}=\tilde{X}\mu_{\theta}+\tilde{V}_{1}^{*}+\psi(1+\psi)^{-1}\tilde{X}(X^{\prime}X)^{-1}X^{\prime}V_{0}^{*}$.

  4. 4.

Set $\tilde{y}^{*}=h\circ g^{-1}(\tilde{z}^{*})$.

Lemma 6 provides direct Monte Carlo sampling from the joint posterior
predictive distribution using standard regression functionals of the covariate
matrix $X$ and the predictive covariate matrix $\tilde{X}$. By comparison, the
samplers in Lemma 2 and Lemma 6 share the same simulation step for
$V_{0}^{*}$, but the second step differs: Lemma 2 requires a draw from a
$p$-dimensional multivariate normal distribution, while Lemma 6 draws from a
$\tilde{n}$-dimensional multivariate normal distribution. As a result, Lemma 6
may provide computational savings relative to Lemma 4 when $\tilde{n}$ is
small, $p$ is large, and the posterior samples of $\theta$ are not needed.

Analytic computations for predictive inference that avoid Monte Carlo sampling
are also available. The predictive distribution is computable using the
marginal density of (1) and (3):

$\mbox{pr}(\tilde{y}\mid y)=\mbox{pr}(\tilde{y},y)/\mbox{pr}(y).$ (14)

The marginal distribution is the constant term in the denominator of (9),

$\mbox{pr}(y)=\bar{\Phi}_{n}\\{\mathcal{C}=g(\mathcal{A}_{y});\mu_{z},\Sigma_{z}\\},$
(15)

where $\mu_{z}$ and $\Sigma_{z}$ are determined by the model (2), e.g.,
$\mu_{z}=X\mu_{\theta}$ and
$\Sigma_{z}=X\Sigma_{\theta}X^{\prime}+\sigma^{2}I_{n}$ in the linear model
setting of Lemma 1. The joint distribution $\mbox{pr}(\tilde{y},y)$ proceeds
similarly via the latent data $(\tilde{z},z)$:
$\mbox{pr}(\tilde{y},y)=\bar{\Phi}_{n+\tilde{n}}\\{\mathcal{C}=g(\mathcal{A}_{\tilde{y}})\times
g(\mathcal{A}_{y});\mu_{\tilde{z}z},\Sigma_{\tilde{z}z}\\},$ where

$\mu_{\tilde{z}z}=\begin{pmatrix}\tilde{X}\mu_{\theta}\\\
X\mu_{\theta}\end{pmatrix},\quad\Sigma_{\tilde{z}z}=\begin{pmatrix}\tilde{X}\Sigma_{\theta}\tilde{X}^{\prime}+\sigma^{2}I_{\tilde{n}}&\tilde{X}\Sigma_{\theta}X^{\prime}\\\
X\Sigma_{\theta}\tilde{X}^{\prime}&X\Sigma_{\theta}X^{\prime}+\sigma^{2}I_{n}\end{pmatrix}$

for the linear model setting of Lemma 1. These direct probability computations
of $\mbox{pr}(\tilde{y}\mid y)$ are most useful when $\tilde{n}$ is small,
such as $\tilde{n}=1$.

Lastly, the predictive probabilities provide analytic point predictions, such
as $E(\tilde{y}\mid y)=\sum_{j\in{\rm supp}(y)}j\ \mbox{pr}(\tilde{y}=j\mid
y).$When $\mbox{supp}(y)=\mathbb{Z}$ or is large, these predictive probability
evaluations may be computationally prohibitive. In that case, the Monte Carlo
simulation is a more appealing option for computing predictive moments and
other functionals.

###  2.5 Nonlinear basis regression

The linear model (3) can accommodate nonlinear regressions via basis
expansions. Consider modeling $z_{i}=z(\tau_{i})$ as a smooth function of
$\tau$ with the following adaptation of (3):

$z_{i}=b^{\prime}(\tau_{i})\theta+\epsilon_{i},\quad\epsilon_{i}\sim
N(0,\sigma^{2})$ (16)

where $b^{\prime}(\tau)=(b_{1}(\tau),\ldots,b_{p}(\tau))$ is a vector of known
basis functions, such as splines, Fourier functions, or wavelets. Canale and
Dunson, (2011) applied a version of model (16) for count data analysis, but
required a Markov chain Monte Carlo algorithm for posterior inference.

Let $X=(b(\tau_{1}),\ldots,b(\tau_{n}))^{\prime}$ denote the basis matrix.
Basis expansions such as (16) are usually accompanied by a smoothness prior of
the form $\theta\sim N_{p}(0,\psi\sigma^{2}\Omega^{-})$, where $\psi>0$
controls the smoothness and $\Omega$ is a known positive semidefinite penalty
matrix. For common choices of $\Omega$, Scheipl et al., (2012) showed that $X$
and $\Omega$ can be reparametrized such that the distribution of $X\theta$ is
unchanged, but $\Omega=I_{p}$ is the identity and $X^{\prime}X$ is diagonal.
Hence, we proceed under the assumption that
$X^{\prime}X=\mbox{diag}\\{d_{j}\\}$ and $\theta\sim
N_{p}(0,\psi\sigma^{2}I_{p})$.

We focus on the predictive distribution $\mbox{pr}\\{\tilde{y}(\tau)\mid y\\}$
across multiple points $\tau$. Monte Carlo simulation proceeds directly
without requiring samples from the posterior:

######  Lemma 7.

Suppose $\theta\sim N_{p}(0,\psi\sigma^{2}I_{p})$ for $\psi>0$ and
$X^{\prime}X=\mbox{diag}\\{d_{j}\\}$. Under model (1) and (16), a joint
predictive draw $\tilde{y}^{*}\sim p(\tilde{y}\mid y)$ at the observation
points $\tau_{1^{\prime}},\ldots,\tau_{\tilde{n}^{\prime}}$ with
$\tilde{X}=(b(\tau_{1^{\prime}}),\ldots,b(\tau_{\tilde{n}^{\prime}}))^{\prime}$
is obtained as follows:

  1. 1.

Simulate $V_{0}^{*}\sim N_{n}(0,\sigma^{2}(\psi XX^{\prime}+I_{n}))$ truncated
to $g(\mathcal{A}_{y})$.

  2. 2.

Simulate $\tilde{V}_{1}^{*}\sim N_{\tilde{n}}(0,\sigma^{2}[\psi\tilde{X}\
\mbox{diag}\\{1/(1+\psi d_{j})\\}\ \tilde{X}^{\prime}+I_{\tilde{n}}])$.

  3. 3.

Compute $\tilde{z}^{*}=\tilde{V}_{1}^{*}+\psi\tilde{X}\ \mbox{diag}\\{1-\psi
d_{j}/(1+\psi d_{j})\\}\ \tilde{X}^{\prime}V_{0}^{*}$.

  4. 4.

Set $\tilde{y}^{*}=h\circ g^{-1}(\tilde{z}^{*})$.

Lemma 7 provides Monte Carlo simulation from the joint predictive distribution
of a discrete and nonlinear regression model. Most notably, the key terms are
computable without any matrix inversions. The most demanding step is again the
simulation of $V_{0}^{*}$ from the $n$-dimensional truncated normal.
Fortunately, this draw is shared among all choices of points
$\tau_{1^{\prime}},\ldots,\tau_{\tilde{n}^{\prime}}$.

More broadly, Lemma 7 is applicable for predictive simulation from any linear
model (1) and (3) with a prior of the form $\theta\sim
N_{p}(0,\psi\sigma^{2}\Omega^{-})$ for a known positive semidefinite matrix
$\Omega$. Although the reparametrization of $X\theta$ changes the
interpretation of $X$ and $\theta$, it does not change the data-generating
process, so the predictive draws remain valid.

###  2.6 Model and variable selection

Suppose we have candidate models $\mathcal{M}_{1},\ldots,\mathcal{M}_{K}$,
where each $\mathcal{M}_{k}$ is of the form (1)–(2). Each candidate model
differs in a fixed attribute: examples include the transformation $g$, the
hyperparameters in the model (such as $\psi$), and linear (3) against
nonlinear (16) models. This setting also applies to variable selection by
replacing each covariate matrix $X$ with $X_{\gamma}$, where
$\gamma\in\\{0,1\\}^{p}$ is a variable inclusion indicator and each $\gamma$
corresponds to a model $\mathcal{M}_{k}$. We focus on two goals: model
selection and prediction via model-averaging.

Given model probabilities $\mbox{pr}(\mathcal{M}_{k})$ for $k=1,\ldots,K$, the
posterior probability of each model is

$\mbox{pr}(\mathcal{M}_{k}\mid
y)=\frac{\mbox{pr}(y\mid\mathcal{M}_{k})\mbox{pr}(\mathcal{M}_{k})}{\sum_{h=1}^{K}\mbox{pr}(y\mid\mathcal{M}_{h})\mbox{pr}(\mathcal{M}_{h})}$
(17)

where $\mbox{pr}(y\mid\mathcal{M}_{k})$ denotes the marginal likelihood under
$\mathcal{M}_{k}$. When the prior on $\theta$ and the latent continuous data
model for $z$ are both Gaussian in the form of (8) for each $\mathcal{M}_{k}$,
we derive the posterior probabilities directly:

######  Lemma 8.

For each model $\\{\mathcal{M}_{k}\\}_{k=1}^{K}$ of the form (1)–(2) with
transformation $g^{(k)}$, suppose that $(\theta^{(k)},z^{(k)})$ are jointly
Gaussian with $z^{(k)}\sim N_{n}(\mu_{z}^{(k)},\Sigma_{z}^{(k)})$. Then the
posterior model probabilities are

$\mbox{pr}(\mathcal{M}_{k}\mid
y)=\frac{\bar{\Phi}_{n}\\{\mathcal{C}=g^{(k)}(\mathcal{A}_{y});\mu_{z}^{(k)},\Sigma_{z}^{(k)}\\}\mbox{pr}(\mathcal{M}_{k})}{\sum_{h=1}^{K}\bar{\Phi}_{n}\\{\mathcal{C}=g^{(h)}(\mathcal{A}_{y});\mu_{z}^{(h)},\Sigma_{z}^{(h)}\\}\mbox{pr}(\mathcal{M}_{h})}.$
(18)

Lemma 8 is naturally informative for model selection, for example by selecting
$\mathcal{M}_{k}$ for which $\mbox{pr}(\mathcal{M}_{k}\mid y)$ is largest.
More broadly, the posterior probabilities in (18) enable model uncertainty in
posterior and predictive inference. Focusing on the latter setting, consider
the goal of prediction that averages over the uncertainty about the models
$\\{\mathcal{M}_{k}\\}_{k=1}^{K}$. Using the predictive simulators, which now
depend on the specific model $\mathcal{M}_{k}$, we propose the following Monte
Carlo sampler: for $s=1,\ldots,S$, sample
$k^{s}\sim\mbox{Categorical}(\pi_{1},\ldots,\pi_{K})$ with
$\pi_{k}=\mbox{pr}(\mathcal{M}_{k}\mid y)$ from (18); then sample
$\tilde{y}^{s}\sim\mbox{pr}(\tilde{y}\mid y,\mathcal{M}=\mathcal{M}_{k^{s}})$
from Lemma 4. This sampler produces joint and independent draws of
$\\{(\mathcal{M}^{s},\tilde{y}^{s})\\}_{s=1}^{S}\sim\mbox{pr}(\mathcal{M},\tilde{y}\mid
y)$. By computing functionals of $\\{\tilde{y}^{s}\\}$ only, we effectively
marginalize over the model uncertainty within
$\\{\mathcal{M}_{k}\\}_{k=1}^{K}$ yet retain a discrete and joint predictive
distribution at $\tilde{X}$.

###  2.7 The sparse means problem for discrete data

The classical sparse means problem considers continuous data
$[z_{i}\mid\theta_{i}]\sim N(\theta_{i},\sigma^{2})$ with the goal of
identifying the nonzero means, $\\{i:\theta_{i}\neq 0\\}$ (Castillo and van
der Vaart,, 2012). Instead, suppose we observe discrete data
$y_{i}\in\mathbb{Z}$ under model (1)–(2). In the context of Example 1, this
scenario corresponds to the normal means model for rounded data. We introduce
sparsity in $\\{\theta_{i}\\}_{i=1}^{n}$ via the spike-and-slab model

$[\theta_{i}\mid\gamma_{i}]\sim\gamma_{i}N(0,\sigma^{2}\psi)+(1-\gamma_{i})\delta_{\\{0\\}}$
(19)

where $\gamma_{i}\in\\{0,1\\}$ is the inclusion indicator with
$\mbox{pr}(\gamma_{i}=1\mid\pi)=\pi$. Conditional on the inclusion indicators,
this model is a special case of (3) with $X=I_{n}$ and $n=p$, where the prior
for the active components $[\theta_{i}\mid\gamma_{i}=1]\sim
N(0,\sigma^{2}\psi)$ is a special case of the Gaussian priors considered
previously. This scenario is also useful for wavelet models and can be
extended to other spike-and-slab priors.

The sparse means model with prior (19) offers certain computational
simplifications, but also introduces new challenges. In particular, inference
on the inclusion indicators $[\gamma\mid y]$ proceeds via Lemma 8. Since the
marginal latent mean is zero and the latent covariance is diagonal,
$\Sigma_{z}^{(\gamma)}=\sigma^{2}\mbox{diag}\\{1+\psi\gamma_{i}\\}_{i=1}^{n}$,
the key term in (18) simplifies considerably:

$\bar{\Phi}_{n}\\{\mathcal{C}=g(\mathcal{A}_{y});0,\Sigma_{z}^{(\gamma)}\\}=\prod_{i=1}^{n}\int_{g(y_{i}-0.5)}^{g(y_{i}+0.5)}\phi_{1}\\{x_{i};0,\sigma^{2}(1+\psi\gamma_{i})\\}\
dx_{i}.$ (20)

In place of an $n$-dimensional integral of a multivariate normal density, we
obtain a product of $n$ univariate normal integrals, which is a massive
advantage for computing with large $n$.

Despite these gains for evaluating each $\mbox{pr}(\gamma\mid y)$, the search
over the model space requires consideration of $2^{n}$ models, which is
computationally prohibitive even for moderate $n$. Instead, we may employ a
stochastic search to sample from $[\gamma\mid y]$. Specifically, we apply a
Gibbs sampler that cycles through $[\gamma_{i}\mid y,\gamma_{-i}]$ for
$i=1,\ldots,n$, which observes that $\mbox{pr}(\gamma_{i}=1\mid
y,\gamma_{-i})=\omega_{i}/(1+\omega_{i})$ for the (full conditional) odds of
inclusion

$\omega_{i}=\frac{\mbox{pr}(\gamma_{i}=1\mid
y,\gamma_{-i})}{\mbox{pr}(\gamma_{i}=0\mid
y,\gamma_{-i})}=\frac{\pi}{1-\pi}\frac{\mbox{pr}(y\mid\gamma_{-i},\gamma_{i}=1)}{\mbox{pr}(y\mid\gamma_{-i},\gamma_{i}=0)}$
(21)

where the terms of the form $\mbox{pr}(y\mid\gamma)$ are given by (20). We may
expand this Gibbs sampler to incorporate a prior on the inclusion probability,
$\pi\sim\mbox{Beta}(a_{\pi},b_{\pi})$ (Scott and Berger,, 2010), which adds a
simple sampling step of the
form$[\pi\mid\gamma]\sim\mbox{Beta}(a_{\pi}+\\#\\{\gamma_{i}=1\\},b_{\pi}+\\#\\{\gamma_{i}=0\\})$
with an accompanying update for $\pi$ in (21) at each Gibbs sweep of
$\gamma_{1},\ldots,\gamma_{n},\pi$.

Inference on the regression coefficients $\theta$ is readily available by
observing that $[\theta,\gamma\mid y]=[\theta\mid y,\gamma][\gamma\mid y]$.
Given samples of $\gamma$ using the above procedure, we sample from
$[\theta\mid y,\gamma]$ using the following fast version of Lemma 2: for
$\\{i:\gamma_{i}=1\\}$, sample $V_{0,i}^{*}\sim N\\{0,\sigma^{2}(1+\psi)\\}$
truncated to $[g(y_{i}-0.5),g(y_{i}+0.5)]$; sample $V_{1,i}^{*}\sim
N\\{0,\sigma^{2}\psi(1+\psi)^{-1}\\}$; and set
$\theta_{i}^{*}=V_{1,i}^{*}+\psi(1+\psi)^{-1}V_{0,i}^{*}$. We set
$\theta_{i}^{*}=0$ for all $\\{i:\gamma_{i}=0\\}$. This sampler is
parallelizable across $i=,1\ldots,n$ yet produces draws from the joint
posterior of $\theta$. The draws of $\theta$ are not needed for sampling the
inclusion indicators $\gamma$, but can be used to obtain predictive samples by
application of Lemma 4.

##  3 Empirical results

###  3.1 Direct Monte Carlo and Gibbs sampling

To demonstrate the computational advantages of direct Monte Carlo sampling for
posterior inference, we compare the proposed sampler from Lemma 2 with the
Gibbs sampler from Kowal and Canale, (2020). The two samplers are applied to
the same linear regression model (3) with a $g$-prior (see Lemma 2). We fix
$\sigma$ at the maximum likelihood estimator (Kowal and Wu,, 2021) and set
$\psi=1000$. The transformation and rounding operators for (1) are given by
Example 2 with $y_{max}=\infty$ and Example 4, respectively.

We simulate negative binomial count data with independent Gaussian covariates
using the design from Kowal and Canale, (2020). The focus here is exclusively
on computational performance. For varying sample sizes $n\in\\{100,200,500\\}$
and number of covariates $p\in\\{10,50\\}$, we run each algorithm for 1000
samples (after a burn-in of 1000 for the Gibbs sampler).

For each simulation, we compute the median effective sample size across the
regression coefficients $\\{\theta_{j}\\}_{j=1}^{p}$ and report the result as
the percent of the total number of simulations (1000). In addition, we record
the raw computing time (using R on a MacBook Pro, 2.8 GHz Intel Core i7). The
results are presented in Figure 1. Most notably, the proposed Monte Carlo
sampler achieves maximal efficiency for all $n,p$ considered. By comparison,
the Gibbs sampler produces 40-70% efficiency for $p=10$, but the efficiency
declines as $p$ increases, including only 5-20% efficiency for $n=100,p=50$.
In addition, the raw computing times are comparable for $n\in\\{100,200\\}$
regardless of $p$, while the proposed Monte Carlo sampler incurs some
additional cost only for $n=500$. In general, the proposed Monte Carlo sampler
demonstrates substantial efficiency gains with increasing $p$, especially when
$n$ is moderate.

Figure 1: Median effective sample size across $\\{\theta_{j}\\}_{j=1}^{p}$
reported as the percent of the total number of simulations (left) and raw
computing time for 1000 samples (right). Results are computed for 25 synthetic
datasets (jittered points) and summarized using medians (non-vertical lines)
with $n\in\\{100,200,500\\}$ and $p\in\\{10,50\\}$.

###  3.2 Prediction for nonlinear regression with synthetic data

We evaluate the predictive distributions from count-valued nonlinear
regression models using simulated data. We consider two data-generating
processes: one with a Negative Binomial likelihood (Neg-Bin) and one that
resembles (1) and (16) (Mixture-CDF). In both cases, we generate a latent
smooth curve $f_{true}(\tau)=\sum_{j=1}^{4}f_{j}(\tau)\beta_{j}$ on a grid of
$n=250$ points, where $\\{f_{j}\\}$ are orthonormal polynomials of degree $j$
and $\\{\beta_{j}\\}$ are independent normal variables with mean zero and
decreasing variance $4-j+1$. For the Negative Binomial case, we set the mean
function to be $\exp\\{f_{true}(\tau)\\}$ and the dispersion to be $10$, which
produces moderately overdispersed counts. We bound these data at $y_{max}=30$
to prevent exceedingly large values. For the second case, we specify the
rounding operator as in Example 2 with $y_{max}=30$. The transformation is
defined akin to Examples 4 or 5: we smooth the distribution function of a
random variable that places $1/2$ mass on a $\mbox{Poisson}(10)$ distribution,
$1/4$ mass uniformly on the set of heaps $\\{5,10,\ldots,25\\}$, and $1/4$
mass uniformly on the set of boundaries $\\{0,30\\}$. The data-generating
model (16) is completed by setting $f_{true}$ to be the mean of the latent
data and $\sigma=0.7$. The simulations are repeated for 100 iterations in each
case.

Monte Carlo samples from the predictive distribution under model (1) and (16)
are generated according to Lemma 7 for the following transformations:
identity, square-root, and log (see Example 3); the empirical cumulative
distribution function transformation (see Example 4); and the Poisson- and
Negative Binomial-based parametric transformations (see Example 5). In each
case, we set $\psi=1000$ and estimate $\sigma$ using the maximum likelihod
estimator. In addition, we compute the predictive distribution from a model-
averaged version of all of these transformations (with equal prior weights)
using the sampler following Lemma 8. To streamline the results, only a subset
of the transformations are presented, yet all are included in the model
averaging. For benchmarking, we include Poisson and Negative Binomial spline
regression models implemented in stan_gamm4 in the rstanarm package in R using
the default specifications.

The predictive distributions are evaluated using ranked probability scores,
which is a proper scoring rule for discrete data predictive distributions
(Czado et al.,, 2009). In addition, we compute mean interval width and
empirical coverage (on out-of-sample testing points) for 90% pointwise
prediction intervals. The results are presented in Figure 2.

Figure 2: Ranked probability score (left) and mean prediction interval widths
with empirical coverage (right) for Negative Binomial data (top) and Mixture-
CDF data (bottom); proposed methods are denoted (*).

For the Negative Binomial setting, the predictive distributions all perform
similarly, with the proposed model-averaged predictions offering slight
improvements over the fixed transformations. The prediction intervals are all
conservative, with the surprising result that the Negative Binomial-based
parametric transformation within the proposed framework offers narrower
intervals than the Negative Binomial model, even though the latter model is a
direct match for the data-generating process. One plausible reason is the
computational challenges of the Negative Binomial sampler: it is difficult to
learn both the dispersion parameter and the nonlinear mean function, and the
output from stan_gamm4 regularly reports low effective sample sizes.

For the second case, the gains of the proposed strategy are now more apparent:
the ranked probability scores especially favor the nonparametric
transformation (Example 4) and the model-averaged predictions. The prediction
intervals confirm these advantages, with an additional warning regarding the
Poisson model: the intervals remain narrower as in the previous case, but the
empirical coverage is now far below the nominal 90%. Again the nonparametric
transformation and the model-averaged predictions perform similarly: the model
probabilities place nearly all of the weight on the nonparametric
transformation, so this result is consistent.

###  3.3 Sparse means estimation for synthetic rounded data

Lastly, we evaluate the selection capabilities of the sparse means model under
rounding. The goal is to assess (i) whether an omission of the rounding
operator impacts selection ability and (ii) whether the transformation is
helpful for this task. Synthetic rounded data are generated as $y_{i}=\lfloor
z_{i}\rfloor$ with $z_{i}=\theta_{i}+\epsilon_{i}$ and $\epsilon_{i}\sim
N(0,1)$, where $\theta_{i}=\mu\gamma_{i}$ and $\mu$ denotes the signal
strength with variable inclusion indicator $\gamma_{i}\in\\{0,1\\}$. We use a
small or moderate proportion of signals,
$\\#\\{\gamma_{i}=1\\}/n\in\\{0.1,0.5\\}$. Results are presented for $\mu=2$
and $n=200$; increasing $\mu\geq 3$ uniformly improved results and $n=500$
produced similar conclusions.

For all competing methods, we apply the same sparsity prior (19) accompanied
by $\pi\sim\mbox{Beta}(1,1)$ and $\psi^{1/2}\sim\mbox{Uniform}(0,n^{1/2})$. We
include two variations of the proposed approach, both with the rounding
operator from Example 1: one uses the identity transformation (as in the data-
generating process) and the other uses the empirical cumulative distribution
function transformation (see Example 4). In addition, we include a Gaussian
model that omits the rounding operation entirely. To maintain focus on the
modeling specifications, we apply identical Gibbs sampling strategies as
outlined in Section 2.7; the Gaussian case simply replaces (21) with the
relevant Gaussian likelihoods. The additional Gibbs sampling step for all
models is $[\psi^{-1}\mid
y,\theta,\gamma]\sim\mbox{Gamma}\\{\\#\\{\gamma_{i}=1\\}/2-1/2,\sum_{i:\gamma_{i}=1}\theta_{i}^{2}/(2\sigma^{2})\\}$
truncated to $[1/n,\infty)$. Lastly, we fix $\sigma$ at the the maximum
likelihod estimator from a kmeans model with two clusters.

The primary inferential target is $\mbox{pr}(\gamma_{i}=1\mid y)$, which we
use to select the nonzero means. The receiving operator characteristic curves
for this selection mechanism (averaged over 100 simulations) and the
corresponding area under the curve values (presented across 100 simulations)
are in Figure 3. Perhaps most surprising, the nonparametric transformation
from Example 4 provides the best selection capabilities by a wide margin.
Comparing the identity models with and without rounding, we find that the
rounding operator offers only minor advantages for selection. Nonetheless, the
best performer belongs to the proposed modeling class (1)–(2), which is a
coherent data-generating process for the discrete (rounded) data.

Figure 3: Receiver operating characteristic curve (left) and corresponding
area under the curve (right) for sparse means selection with rounded data;
proposed methods are denoted (*).

##  4 Discussion

Despite the generality of the linear model (1) and (3), the posterior and
predictive derivations treat the variance components as known. Our analyses
estimate $\sigma$ using maximum likelihood, although other estimators are
available. Alternatively, we may place a prior on $\sigma$ and marginalize
over this unknown parameter. When $\sigma^{2}$ is assigned an inverse-Gamma
prior, the joint distribution of $(\theta,z)$ as in (8), but unconditional on
$\sigma^{2}$, is a multivariate $t$-distribution. The posterior distribution
$p(\theta\mid y)$ is therefore a _selection $t$-distribution_, for which the
density and distribution functions are given in Arellano-Valle et al., (2006).

For more general regression scenarios with $\epsilon\sim
N_{n}(0,\Sigma_{\epsilon})$ or with hyperpriors on
$(\mu_{\theta},\Sigma_{\theta})$ such as for hierarchical modeling, the
posterior derivations remain relevant and provide a blocked Gibbs sampling
strategy. Let $\Sigma=(\Sigma_{\epsilon},\mu_{\theta},\Sigma_{\theta})$ denote
the parameters excluding $\theta$. A blocked Gibbs sampler proceeds
iteratively by drawing $[z,\theta\mid y,\Sigma]=[\theta\mid y,\Sigma][z\mid
y,\theta,\Sigma]$ and $[\Sigma\mid y,z,\theta]$. The first block uses the
selection normal posterior for $\theta$, while $[z\mid y,\theta,\Sigma]\sim
N_{n}(X\theta,\Sigma_{\epsilon})$ truncated to $g(\mathcal{A}_{y})$. When
$\Sigma_{\epsilon}$ is diagonal, this step only requires $n$ draws from
independent and univariate truncated normals, which does not add much
computational cost. The remaining block typically satisfies $[\Sigma\mid
y,z,\theta]=[\Sigma\mid z,\theta]$ and therefore may use existing sampling
steps for variance parameters in (continuous) Gaussian models. Importantly,
our selection normal results provide a joint sampling step for $[z,\theta\mid
y,\Sigma]$, which protects against Markov chain Monte Carlo sampling
inefficiencies due to posterior correlations among $z$ and $\theta$.

In concurrent research, King and Kowal, (2021) introduce a dynamic linear
model in (2) for count time series data, which extends the dynamic probit
results of Fasano et al., (2021). The conjugacy afforded by Lemma 3 enables
analytic recursions for filtering and smoothing (selection normal)
distributions, which are accompanied by an optimal particle filter for online
inference. Similar conjugacy results and computational advantages may be
available for spatial, functional, or other dependent integer-valued data.

## References

  * Arellano-Valle et al., (2006) Arellano-Valle, R. B., Branco, M. D., and Genton, M. G. (2006).  A unified view on skewed distributions arising from selections.  Canadian Journal of Statistics, 34(4):581–601.
  * Botev, (2017) Botev, Z. I. (2017).  The normal law under linear restrictions: simulation and estimation via minimax tilting.  Journal of the Royal Statistical Society. Series B: Statistical Methodology, 79(1):125–148.
  * Box and Cox, (1964) Box, G. E. P. and Cox, D. R. (1964).  An analysis of transformations.  Journal of the Royal Statistical Society: Series B (Methodological), 26(2):211–243.
  * Bradley et al., (2018) Bradley, J. R., Holan, S. H., and Wikle, C. K. (2018).  Computationally efficient multivariate spatio-temporal models for high-dimensional count-valued data (with discussion).  Bayesian Analysis, 13(1):253–310.
  * Canale and Dunson, (2011) Canale, A. and Dunson, D. B. (2011).  Bayesian kernel mixtures for counts.  Journal of the American Statistical Association, 106(496):1528–1539.
  * Canale and Dunson, (2013) Canale, A. and Dunson, D. B. (2013).  Nonparametric Bayes modelling of count processes.  Biometrika, 100(4):801–816.
  * Castillo and van der Vaart, (2012) Castillo, I. and van der Vaart, A. (2012).  Needles and straw in a haystack: Posterior concentration for possibly sparse sequences.  The Annals of Statistics, 40(4):2069–2101.
  * Czado et al., (2009) Czado, C., Gneiting, T., and Held, L. (2009).  Predictive model assessment for count data.  Biometrics, 65(4):1254–1261.
  * Dempster and Rubin, (1983) Dempster, A. P. and Rubin, D. B. (1983).  Rounding error in regression: The appropriateness of Sheppard’s corrections.  Journal of the Royal Statistical Society: Series B (Methodological), 45(1):51–59.
  * Durante, (2019) Durante, D. (2019).  Conjugate Bayes for probit regression via unified skew-normal distributions.  Biometrika, 106(4):765–779.
  * Fasano and Durante, (2020) Fasano, A. and Durante, D. (2020).  A class of conjugate priors for multinomial probit models which includes the multivariate normal one.  arXiv preprint arXiv:2007.06944.
  * Fasano et al., (2021) Fasano, A., Rebaudo, G., Durante, D., and Petrone, S. (2021).  A closed-form filter for binary time series.  Statistics and Computing, 31(4):1–20.
  * Feldman and Kowal, (2021) Feldman, J. and Kowal, D. (2021).  A Bayesian Framework for Generation of Fully Synthetic Mixed Datasets.  Annals of Applied Statistics (revision submitted).
  * Fritsch and Carlson, (1980) Fritsch, F. N. and Carlson, R. E. (1980).  Monotone piecewise cubic interpolation.  SIAM Journal on Numerical Analysis, 17(2):238–246.
  * Heitjan, (1989) Heitjan, D. F. (1989).  Inference from grouped continuous data: a review.  Statistical science, 4(2):164–179.
  * Heitjan and Rubin, (1991) Heitjan, D. F. and Rubin, D. B. (1991).  Ignorability and coarse data.  The annals of statistics, pages 2244–2253.
  * Hoff, (2007) Hoff, P. D. (2007).  Extending the rank likelihood for semiparametric copula estimation.  The Annals of Applied Statistics, 1(1):265–283.
  * Jia et al., (2021) Jia, Y., Kechagias, S., Livsey, J., Lund, R., and Pipiras, V. (2021).  Latent Gaussian Count Time Series.  Journal of the American Statistical Association, (just-accepted):1–28.
  * King and Kowal, (2021) King, B. and Kowal, D. R. (2021).  Warped Dynamic Linear Models for Time Series of Counts.  arXiv preprint arXiv:2110.14790.
  * Kowal and Canale, (2020) Kowal, D. R. and Canale, A. (2020).  Simultaneous Transformation and Rounding (STAR) Models for Integer-Valued Data.  Electronic Journal of Statistics, 14(1):1744–1772.
  * Kowal and Wu, (2021) Kowal, D. R. and Wu, B. (2021).  Semiparametric count data regression for self-reported mental health.  arXiv preprint arXiv:2106.09114.
  * Liu et al., (2009) Liu, H., Lafferty, J., and Wasserman, L. (2009).  The nonparanormal: Semiparametric estimation of high dimensional undirected graphs.  Journal of Machine Learning Research, 10(10):2295–2328.
  * Peluso et al., (2019) Peluso, A., Vinciotti, V., and Yu, K. (2019).  Discrete Weibull generalized additive model: an application to count fertility data.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 68(3):565–583.
  * Pitt et al., (2006) Pitt, M., Chan, D., and Kohn, R. (2006).  Efficient Bayesian inference for Gaussian copula regression models.  Biometrika, 93(3):537–554.
  * Polson et al., (2013) Polson, N. G., Scott, J. G., and Windle, J. (2013).  Bayesian inference for logistic models using Pólya-Gamma latent variables.  Journal of the American Statistical Association, 108(504):1339–1349.
  * Scheipl et al., (2012) Scheipl, F., Fahrmeir, L., and Kneib, T. (2012).  Spike-and-slab priors for function selection in structured additive regression models.  Journal of the American Statistical Association, 107(500):1518–1532.
  * Scott and Berger, (2010) Scott, J. G. and Berger, J. O. (2010).  Bayes and empirical-Bayes multiplicity adjustment in the variable-selection problem.  The Annals of Statistics, pages 2587–2619.
  * Sellers and Shmueli, (2010) Sellers, K. F. and Shmueli, G. (2010).  A flexible regression model for count data.  The Annals of Applied Statistics, pages 943–961.
  * Siegfried and Hothorn, (2020) Siegfried, S. and Hothorn, T. (2020).  Count transformation models.  Methods in Ecology and Evolution, 11(7):818–827.
  * Zellner, (1986) Zellner, A. (1986).  On assessing prior distributions and Bayesian regression analysis with g-prior distributions.  Bayesian inference and decision techniques.

## Supplementary materials

Supplementary materials include R code to reproduce the simulation studies
from Section 3 and an R package available at www.GitHub.com/drkowal/rSTAR.

## Appendix

###### Proof.

_(Theorem 2)_ Since Theorem 1 implies that $p(\theta\mid y)$ is a selection
distribution under model (1)–(2), the joint Gaussian assumption in (8) implies
that $p(\theta\mid y)$ is selection normal. The remaining properties follow
from Arellano-Valle et al., (2006). ∎

###### Proof.

_(Theorem 3)_ The result follows from the properties of the selection normal
(Arellano-Valle et al.,, 2006) applied to model (1) and (8) via Theorem 2. ∎

###### Proof.

_(Lemma 1)_ The result follows by direct computation as a consequence of
Theorem 2 under the stated model (3) and prior $\theta\sim
N(\mu_{\theta},\Sigma_{\theta})$. ∎

###### Proof.

_(Lemma 2)_ Using Theorem 3, it is sufficient to compute $\Sigma_{z}$ and
$\Sigma_{\theta}-\Sigma_{z\theta}^{\prime}\Sigma_{z}^{-1}\Sigma_{z\theta}$
under the linear model with the $g$-prior. By computing
$\Sigma_{z}=\sigma^{2}\\{I_{n}+\psi X(X^{\prime}X)^{-1}X^{\prime}\\}$ and
applying the Woodbury identity, it follows that
$\Sigma_{z}^{-1}=\sigma^{-2}\\{I_{n}-\psi(1+\psi)^{-1}X(X^{\prime}X)^{-1}X^{\prime}\\}$,
so
$\Sigma_{z\theta}^{\prime}\Sigma_{z}^{-1}=\psi(1+\psi)^{-1}(X^{\prime}X)^{-1}X^{\prime}$
and
$\Sigma_{z\theta}^{\prime}\Sigma_{z}^{-1}\Sigma_{z\theta}=\sigma^{2}\psi^{2}(1+\psi)^{-1}(X^{\prime}X)^{-1}$.
Observing that
$\Sigma_{\theta}-\Sigma_{z\theta}^{\prime}\Sigma_{z}^{-1}\Sigma_{z\theta}=\sigma^{2}\psi(1+\psi)^{-1}(X^{\prime}X)^{-1}$
shows the result. ∎

###### Proof.

_(Lemma 3)_ The selection normal prior is equivalently defined by $[\theta\mid
z_{0}\in\mathcal{C}_{0}]$ for $(z_{0}^{\prime},\theta^{\prime})^{\prime}$
jointly Gaussian with moments given in the prior. The posterior is constructed
similarly: $[\theta\mid y]=[\theta\mid z_{0}\in\mathcal{C}_{0},z\in
g(\mathcal{A}_{y})]=[\theta\mid z\in\mathcal{C}_{0}\times g(\mathcal{A}_{y})]$
where $z_{1}=(z_{0}^{\prime},z^{\prime})^{\prime}$. It remains to identify the
moments of
$(z_{1}^{\prime},\theta^{\prime})^{\prime}=(z_{0}^{\prime},z^{\prime},\theta^{\prime})^{\prime}$.
For each individual block of $z_{0}$, $z$, and $\theta$ and the pairs
($z_{0},\theta$) and ($z,\theta$), the moments are provided by the selection
normal prior or the selection normal posterior in Lemma 1. All that remains is
the cross-covariance
$\mbox{Cov}(z_{0},z)=\mbox{Cov}(z_{0},X\theta+\epsilon)=\Sigma_{z_{0}\theta}X^{\prime}$.
∎

###### Proof.

_(Lemma 5)_ As a preliminary, consider a linear combination of a selection
normal variable with a Gaussian innovation:

######  Lemma 9.

Suppose
$\theta\sim\mbox{SLCT-N}_{n,p}(\mu_{z},\mu_{\theta},\Sigma_{z},\Sigma_{\theta},\Sigma_{z\theta},\mathcal{C})$
and let $\theta^{*}=A\theta+a$ where $A$ is a fixed $q\times p$ matrix and
$a\sim N_{q}(\mu_{a},\Sigma_{a})$ is independent of $\theta$. Then
$\theta^{*}\sim\mbox{SLCT-N}_{n,q}(\mu_{z},\mu_{\theta^{*}}=A\mu_{\theta}+\mu_{a},\Sigma_{z},\Sigma_{\theta^{*}}=A\Sigma_{\theta}A^{\prime}+\Sigma_{a},\Sigma_{z\theta^{*}}=\Sigma_{z\theta}A^{\prime},\mathcal{C})$.

###### Proof.

_(Lemma 9)_ First consider the distribution of $A\theta$. Since the moments of
the joint distribution $(z,\theta)$ are given by assumption, it follows that
$(z,A\theta)$ is jointly Gaussian and $[A\theta\mid
z\in\mathcal{C}]\sim\mbox{SLCT-N}_{n,q}(\mu_{z},A\mu_{\theta},\Sigma_{z},A\Sigma_{\theta}A^{\prime},\Sigma_{z\theta}A^{\prime},\mathcal{C})$.
Using the moment generating functions and noting independence between
$A\theta$ and $a$, it follows that
$M_{\theta^{*}}(s)=M_{A\theta+a}(s)=M_{A\theta}(s)M_{a}(s)$ where
$M_{A\theta}$ is in Theorem 2 and
$M_{a}(s)=\exp(s^{\prime}\mu_{a}+s^{\prime}\Sigma_{a}s/2)$. The result of this
product is the moment generating function of the stated selection normal
distribution. ∎

Returning to Lemma 5, the latent predictive variable $[\tilde{z}\mid\theta]$
in (13) can be represented as $\tilde{z}=\tilde{X}\theta+\tilde{\epsilon}$,
where $\tilde{\epsilon}\sim N_{n}(0,\sigma^{2}I_{n})$ is independent of
$\\{\epsilon_{i}\\}_{i=1}^{n}$ from (3). The integrand in (13) requires
marginalization over $\theta\sim p(\theta\mid y)$, which is selection normal
according to Lemma 1. Hence we obtain the latent predictive distribution
$[\tilde{z}\mid y]$ via the marginalization from Lemma 9. ∎

###### Proof.

_(Lemma 6)_ First, $\Sigma_{z}=\sigma^{2}\\{I_{n}+\psi
X(X^{\prime}X)^{-1}X^{\prime}\\}$ is unchanged from Lemma 2. Similarly, using
the generic notation from (8), we have
$\Sigma_{z\theta}^{\prime}\Sigma_{z}^{-1}=\psi(1+\psi)^{-1}\tilde{X}(X^{\prime}X)^{-1}X^{\prime}$
and
$\Sigma_{z\theta}^{\prime}\Sigma_{z}^{-1}\Sigma_{z\theta}=\sigma^{2}\psi^{2}(1+\psi)^{-1}\tilde{X}(X^{\prime}X)^{-1}\tilde{X}^{\prime}$.
Computing
$\Sigma_{\theta}-\Sigma_{z\theta}^{\prime}\Sigma_{z}^{-1}\Sigma_{z\theta}=\sigma^{2}\\{\psi(1+\psi)^{-1}\tilde{X}(X^{\prime}X)^{-1}\tilde{X}^{\prime}+I_{\tilde{n}}\\}$
shows the result for $\tilde{z}^{*}$. Since the link in (1) is deterministic,
the result for $\tilde{y}^{*}$ follows. ∎

###### Proof.

_(Lemma 7)_ Computing the selection normal terms in the notation of (8), we
have $\Sigma_{z}=\sigma^{2}(\psi XX^{\prime}+I_{n})$,
$\Sigma_{\theta}=\sigma^{2}(\psi\tilde{X}\tilde{X}^{\prime}+I_{\tilde{n}})$,
and $\Sigma_{z\theta}=\psi\sigma^{2}X\tilde{X}^{\prime}$. Let
$D_{\psi}=(\psi^{-1}I_{p}+X^{\prime}X)^{-1}=\mbox{diag}\\{\psi/(1+\psi
d_{j})\\}$, so $\Sigma_{z}^{-1}=\sigma^{-2}(I_{n}-XD_{\psi}X^{\prime})$ by the
Woodbury identity. Next, we derive
$\Sigma_{z\theta}^{\prime}\Sigma_{z}^{-1}=\psi\tilde{X}\\{X^{\prime}-(X^{\prime}X)D_{\psi}X^{\prime}\\}$
and
$\Sigma_{z\theta}^{\prime}\Sigma_{z}^{-1}\Sigma_{z\theta}=\psi^{2}\sigma^{2}\tilde{X}\\{(X^{\prime}X)-(X^{\prime}X)D_{\psi}(X^{\prime}X)\\}\tilde{X}^{\prime}$,
so
$\Sigma_{\theta}-\Sigma_{z\theta}^{\prime}\Sigma_{z}^{-1}\Sigma_{z\theta}=\psi\sigma^{2}\tilde{X}[I_{p}-\psi\\{(X^{\prime}X)-(X^{\prime}X)D_{\psi}(X^{\prime}X)\\}]\tilde{X}^{\prime}+\sigma^{2}I_{\tilde{n}}$.
Since $X^{\prime}X=\mbox{diag}(d_{j})$, the term in the square brackets is
diagonal with elements $1-\psi\\{d_{j}-d_{j}^{2}\psi/(1+\psi
d_{j})\\}=1/(1+\psi d_{j})$. Applying Lemma 6 shows the result. ∎

###### Proof.

_(Lemma 8)_ The result follows from (15) applied to each model
$\mathcal{M}_{k}$. ∎


$\mu_{1}$ is a discretization of Brownian motion in the real case or a
weighted uniform measure on a particular finite generating set in the non-
Archimedean case. It satisfies
$\delta(\mu_{1})=\delta(\Gamma)=\frac{\dim_{\mathrm{Haus}}(\partial X)}{2}$.

Consider a discrete subgroup $\Delta$ of the Lie group $G$ and assume that
$\delta(\Delta)\leq\frac{\delta(\mu_{1})}{2}$. Let $\nu$ be any weak-$*$
accumulation point of the sequence of probability measures
$\nu_{n}=\frac{1}{n}\sum_{i=1}^{n}\mu_{\textrm{0}}^{*i}\ast\mathrm{D}_{\Delta}$.
The key property of the probability measure $\mu_{0}$ implies that the
resulting random subgroup $\nu$ is _discrete_ and $\mu_{0}$-stationary [FG23,
Theorem 1.6]. The last paragraph of the preceding discussion shows that the
random subgroup $\nu$ is at the same time $\mu_{1}$-stationary.

We remark that in the zero characteristic non-Archimedean case there is no
need to consider the measure $\mu_{0}$ at all and we may work with the measure
$\mu_{1}$ directly.

We know from Theorem 9.2 that $\nu$-almost every _non-trivial_ subgroup
$\Lambda$ satisfies $\delta(\Lambda)>\frac{\delta(\mu_{1})}{2}$. One the other
hand, the critical exponent is a conjugation invariant and a lower semi-
continuous function on the Chabauty space of discrete subgroups by Proposition
10.4. So $\delta(\Lambda)\leq\delta(\Delta)$ must hold $\nu$-almost surely. We
conclude that $\nu=\mathrm{D}_{\\{e\\}}$ where $\\{e\\}$ is the trivial
subgroup of the Lie group $G$. In the particular the trivial subgroup
$\\{e\\}$ lies is the Chabauty closure of the conjugacy class of the subgroup
$\Delta$ by the Portmanteau theorem. This is equivalent to saying that the
subgroup $\Delta$ is not confined, see Proposition 10.2 for details. ∎

### General Gromov hyperbolic groups

In the case of free groups we recover the result of [Fra20, Corollary 3.2].

######  Theorem 10.6.

Let $F_{k}$ be the free group of rank $k\in\mathbb{N}$ and $X$ be the Cayley
graph of $F_{k}$ with respect to the standard generating set. Let $H$ be a
subgroup of $F_{k}$ with $\delta(H)\leq\frac{\log(2k-1)}{2}$. Then $H$ is not
confined.

###### Proof.

The proof follows the same strategy as the proof of Theorem 10.5. This time we
use the uniform probability measure $\mu_{k}$ supported on the standard
symmetric generating set for the group $F_{k}$. It is shown in [MM07] that the
measure $\mu_{k}$ satisfies

(10.5) $\delta(\mu_{k})=\delta(\Gamma)=\dim_{\mathrm{Haus}}(\partial
X)=\log(2k-1).$

The proof is simpler in this case as there is no need to worry about
discreteness of the limiting stationary random subgroup (for the group $F_{k}$
itself is discrete). ∎

Finally we state a weaker analogue of Theorems 10.5 and 10.6 that applies more
generally to all Gromov hyperbolic metric spaces.

######  Theorem 10.7.

Let $X$ be a proper Gromov hyperbolic geodesic metric space with isometry
group $\mathrm{Is}(X)$. Assume that the group $\mathrm{Is}(X)$ admits a
discrete subgroup $\Gamma$ which acts properly and cocompactly and which
intersects the elliptic radical $\mathrm{E}(\mathrm{Is}(X))$ trivially. If
$\Delta$ is any subgroup of $\Gamma$ with
$\delta(\Delta)<\frac{\delta(\Gamma)}{2}$ then $\Delta$ is not confined.

The above result applies in particular in the case where the group
$\mathrm{Is}(X)$ itself is discrete.

###### Proof of Theorem 10.7.

Let $\mu_{k}$ be the uniform probability measure supported on the finite set

(10.6) $\\{\gamma\in\Gamma\>:\>k\leq\|\gamma\|<k+1\\}$

for every $k\in\mathbb{N}$. The sequence of probability measures $\mu_{k}$
satisfies $\delta(\mu_{k})\xrightarrow{k\to\infty}\delta(\Gamma)$ [GMM18,
Theorem 1.4]. In particular $\delta(\Delta)\leq\frac{\delta(\mu_{k_{0}})}{2}$
holds for some fixed sufficiently large $k_{0}\in\mathbb{N}$. Let $\nu$ be any
weak-$*$ accumulation point of the sequence of probability measures
$\nu_{n}=\frac{1}{n}\sum_{i=1}^{n}\mu_{k_{0}}^{*i}\ast\mathrm{D}_{\Delta}$.
Since the convolution powers $\mu_{k_{0}}^{*i}$ are all supported on the
discrete subgroup $\Gamma$ and as $\Delta\leq\Gamma$, the resulting
$\mu_{k_{0}}$-stationary random subgroup $\nu$ is almost surely contained in
the discrete subgroup $\Gamma$. From this point we may conclude exactly as in
the proof of Theorem 10.5. ∎

######  Remark 10.8.

To replace the assumption $\delta(\Delta)<\frac{\delta(\Gamma)}{2}$ in Theorem
10.7 by a non-strict inequality, we require a probability measure $\mu$ on the
group $\Gamma$ whose support generates $\Gamma$ as a semigroup with finite
$p$-th moment for some $p>1$ and $\delta(\mu)=\delta(\Gamma)$. When $X$ is a
$CAT(-1)$-space, a measure with finite _first_ moment satisfying the other
required conditions was constructed by Connell and Muchnik [CM07]. It is
unclear to the authors whether this measure has finite $p$-th moment for some
$p>1$.

## References

  * [ABB+20] Miklos Abert, Nicolas Bergeron, Ian Biringer, Tsachik Gelander, Nikolay Nikolov, Jean Raimbault, and Iddo Samet.  On the growth of ${L^{2}}$-invariants of locally symmetric spaces, II: Exotic invariant random subgroups in rank one.  International Mathematics Research Notices, 2020(9):2588–2625, 2020\.
  * [AGV14] Miklós Abért, Yair Glasner, and Bálint Virág.  Kesten’s theorem for invariant random subgroups.  Duke Mathematical Journal, 163(3):465–488, 2014.
  * [BF] Uri Bader and Alex Furman.  Some ergodic properties of metrics on hypebolic groups.
  * [BH13] Martin Bridson and André Haefliger.  Metric spaces of non-positive curvature, volume 319.  Springer Science & Business Media, 2013.
  * [BHM08] Sébastien Blachère, Peter Haïssinsky, and Pierre Mathieu.  Asymptotic entropy and Green speed for random walks on countable groups.  In Annals of Probability, volume 36, pages 1134–1152, 2008.
  * [BHM11] Sébastien Blachère, Peter Haïssinsky, and Pierre Mathieu.  Harmonic measures versus quasiconformal measures for hyperbolic groups.  In Annales scientifiques de l’Ecole normale supérieure, volume 44, pages 683–721, 2011.
  * [BL] Werner Ballmann and François Ledrappier.  Discretization of positive harmonic functions on Riemannian manifolds and Martin boundary.
  * [BN15] Lewis Bowen and Amos Nevo.  Von Neumann and Birkhoff ergodic theorems for negatively curved groups.  In Annales scientifiques de lecole normale superieure, volume 48, pages 1113–1147, 2015.
  * [Bor63] Armand Borel.  Compact Clifford-Klein forms of symmetric spaces.  Topology, 2(1-2):111–122, 1963.
  * [BQ16a] Yves Benoist and Jean-François Quint.  Central limit theorem on hyperbolic groups.  Izvestiya Rossiiskoi Akademii Nauk. Seriya Matematicheskaya, 80(1):5–26, 2016.
  * [BQ16b] Yves Benoist and Jean-François Quint.  Random walks on reductive groups.  Springer, 2016.
  * [BS06] Uri Bader and Yehuda Shalom.  Factor and normal subgroup theorems for lattices in products of groups.  Inventiones mathematicae, 163(2):415–454, 2006.
  * [BT17] Ian Biringer and Omer Tamuz.  Unimodularity of invariant random subgroups.  Transactions of the American Mathematical Society, 369(6):4043–4061, 2017.
  * [BTMT12] Petra Bonfert-Taylor, Katsuhiko Matsuzaki, and Edwars Taylor.  Large and small covers of a hyperbolic manifold.  J. Geometric Analysis, 22:455–470, 2012.
  * [Car01] Lisa Carbone.  On the classification of rank one groups over non-Archimedean local fields.  In https://sites.math.rutgers.edu/ carbonel/pdfs/Classification.pdf. Citeseer, 2001\.
  * [CdCMT15] Pierre-Emmanuel Caprace, Yves de Cornulier, Nicolas Monod, and Romain Tessera.  Amenable hyperbolic groups.  Journal of the european Mathematical Society, 17(11):2903–2947, 2015.
  * [CM07] Christopher Connell and Roman Muchnik.  Harmonicity of quasiconformal measures and Poisson boundaries of hyperbolic spaces.  Geometric and Functional Analysis, 17:707–769, 2007.
  * [CM09a] Pierre-Emmanuel Caprace and Nicolas Monod.  Isometry groups of non-positively curved spaces: discrete subgroups.  Journal of Topology, 2(4):701–746, 2009.
  * [CM09b] Pierre-Emmanuel Caprace and Nicolas Monod.  Isometry groups of non-positively curved spaces: structure theory.  Journal of topology, 2(4):661–700, 2009.
  * [Coo93] Michel Coornaert.  Mesures de Patterson-Sullivan sur le bord d’un espace hyperbolique au sens de Gromov.  Pacific Journal of Mathematics, 159(2):241–270, 1993.
  * [Cor90] Kevin Corlette.  Hausdorff dimensions of limit sets I.  Inventiones mathematicae, 102(1):521–541, 1990.
  * [DGLL15] Bruno Duchesne, Yair Glasner, Nir Lazarovich, and Jean Lécureux.  Geometric density for invariant random subgroups of groups acting on $\mathrm{CAT}(0)$ spaces.  Geometriae Dedicata, 175(1):249–256, 2015.
  * [DSU17] Tushar Das, David Simmons, and Mariusz Urbański.  Geometry and dynamics in Gromov hyperbolic metric spaces, volume 218.  American Mathematical Soc., 2017.
  * [Dus] Matthieu Dussaule.  Groupes d’isométries et marches aléatoires dans les espaces hyperboliques: Rapport de stage de M2, ENS Lyon.
  * [Els74] Jürgen Elstrodt.  Die resolvente zum eigenwertproblem der automorphen formen in der hyperbolischen ebene. Teil III.  Mathematische Annalen, 208:99–132, 1974.
  * [EO23] Sam Edwards and Hee Oh.  Temperedness of $L^{2}(\Gamma\backslash G)$ and positive eigenfunctions in higher rank, 2023.
  * [FG23] Mikolaj Fraczyk and Tsachik Gelander.  Infinite volume and infinite injectivity radius.  Annals of Mathematics, 197(1):389–421, 2023.
  * [Fra20] Mikolaj Fraczyk.  Kesten’s theorem for uniformly recurrent subgroups.  Ergodic Theory and Dynamical Systems, 40(10):2778–2787, 2020.
  * [FT22] Behrang Forghani and Giulio Tiozzo.  Shannon’s theorem for locally compact groups.  The Annals of Probability, 50(1):61–89, 2022.
  * [Fur71] Harry Furstenberg.  Random walks and discrete subgroups of Lie groups.  Advances in probability and related topics, 1:1–63, 1971.
  * [Fur02] Alex Furman.  Random walks on groups and random transformations.  In Handbook of dynamical systems, volume 1, pages 931–1014. Elsevier, 2002.
  * [Gel15] Tsachik Gelander.  A lecture on invariant random subgroups.  arXiv preprint arXiv:1503.08402, 2015.
  * [GL19] Ilya Gekhtman and Arie Levit.  Critical exponents of invariant random subgroups in negative curvature.  Geometric and Functional Analysis, 29(2):411–439, 2019.
  * [GLM22] Tsachik Gelander, Arie Levit, and Gregori A. Margulis.  Effective discreteness radius of stabilizers for stationary actions.  Michigan Mathematical Journal, 72(1):389–438, 2022.
  * [GMM18] Sébastien Gouëzel, Frédéric Mathéus, and François Maucourant.  Entropy and drift in word hyperbolic groups.  Inventiones mathematicae, 211(3):1201–1255, 2018.
  * [Gro87] Mikhael Gromov.  Hyperbolic groups.  In Essays in group theory, pages 75–263. Springer, 1987.
  * [GT20] Ilya Gekhtman and Giulio Tiozzo.  Entropy and drift for Gibbs measures on geometrically finite manifolds.  Transactions of the American Mathematical Society, 373:2949–2980, 2020.
  * [Gui80] Yves Guivarch.  Sur la loi des grands nombres et le rayon spectral d’une marche aleatoire.  Asterisque, 74:47–98, 1980.
  * [Jaw95] Wojciech Jaworski.  Strong approximate transitivity, polynomial growth, and spread out random walks on locally compact groups.  Pacific Journal of Mathematics, 170(2):517–533, 1995.
  * [Kai94] Vadim A Kaimanovich.  The Poisson boundary of hyperbolic groups.  Comptes Rendus de l’Academie des Sciences-Serie I-Mathematique, 318(1):59–64, 1994.
  * [Kai96] Vadim A Kaimanovich.  Boundaries of invariant Markov operators: the identification problem.  London Mathematical Society lecture note series, pages 127–176, 1996.
  * [Kak51] Shizuo Kakutani.  Random ergodic theorems and Markoff processes with a stable distribution.  Technical report, Yale University, New Haven United States, 1951.
  * [Kec12] Alexander Kechris.  Classical descriptive set theory, volume 156.  Springer Science & Business Media, 2012.
  * [Kle99] Bruce Kleiner.  The local structure of length spaces with curvature bounded above.  Mathematische Zeitschrift, 231:409–456, 1999.
  * [KM99] Anders Karlsson and Gregory A Margulis.  A multiplicative ergodic theorem and nonpositively curved spaces.  Communications in mathematical physics, 208(1):107–123, 1999.
  * [KV83] Vadim A Kaimanovich and Anatoly M Vershik.  Random walks on discrete groups: boundary and entropy.  The annals of probability, 11(3):457–490, 1983.
  * [Led90] Francois Ledrappier.  Harmonic measures and Bowen-Margulis measures.  Israel J. Math., 71(3):275–287, 1990.
  * [Led01] François Ledrappier.  Some asymptotic properties of random walks on free groups.  Topics in probability and Lie groups: boundary theory, 28(117-152):7, 2001.
  * [LS84] Terry Lyons and Dennis Sullivan.  Function theory, random paths and covering spaces.  Journal of differential geometry, 19(2):299–323, 1984.
  * [Lub94] Alex Lubotzky.  Discrete groups, expanding graphs and invariant measures, volume 125.  Springer Science & Business Media, 1994.
  * [Mar91] Gregori A Margulis.  Discrete subgroups of semisimple Lie groups, volume 17.  Springer Science & Business Media, 1991.
  * [McM99] Curtis McMullen.  Hausdorff dimension and conformal dynamics. I. Strong convergence of Kleinian groups.  Journal of Differential Geometry, 51(3):471–515, 1999.
  * [MM07] Jean Mairesse and Frédéric Mathéus.  Random walks on free products of cyclic groups.  Journal of the London Mathematical Society, 75(1):47–66, 2007.
  * [MT18] Joseph Maher and Giulio Tiozzo.  Random walks on weakly hyperbolic groups.  Journal für die reine und angewandte Mathematik (Crelles Journal), 2018(742):187–239, 2018.
  * [MYJ20] Katsuhiko Matsuzaki, Yasuhiro Yabuki, and Johannes Jaerisch.  Normalizer, divergence type, and Patterson measure for discrete groups of the Gromov hyperbolic space.  Groups, Geometry, and Dynamics, 14(2):369–411, 2020.
  * [Nic89] Peter Nicholls.  The ergodic theory of discrete groups, volume 143.  Cambridge University Press, 1989.
  * [Osi17] D Osin.  Invariant random subgroups of groups acting on hyperbolic spaces.  Proceedings of the American Mathematical Society, 145(8):3279–3288, 2017.
  * [Pat76] Samuel James Patterson.  The limit set of a Fuchsian group.  Acta Mathematica, 136:241–273, 1976.
  * [Ser09] Jean-Pierre Serre.  Lie algebras and Lie groups: 1964 lectures given at Harvard University.  Springer, 2009.
  * [Sul87] Dennis Sullivan.  Related aspects of positivity in Riemannian geometry.  Journal of differential geometry, 25(3):327–351, 1987.
  * [Tan17] Ryokichi Tanaka.  Hausdorff spectrum of harmonic measure.  Ergodic Theory and Dynamical Systems, 37(1):277–307, 2017.
  * [Tan19] Ryokichi Tanaka.  Dimension of harmonic measures in hyperbolic spaces,.  Ergodic Theory and Dynamical Systems, 39(2):474–499, 2019.
  * [Ver00] Anatoly Vershik.  Dynamic theory of growth in groups: Entropy, boundaries, examples.  Uspekhi Mat. Nauk, 55:59–128, 2000.
  * [Woe00] Wolfgang Woess.  Random walks on infinite graphs and groups, volume 138.  Cambridge Tracts in Mathematics, 2000.
  * [Zim13] Robert J Zimmer.  Ergodic theory and semisimple groups, volume 81.  Springer Science & Business Media, 2013.