# Gibbs posterior concentration rates under sub-exponential type losses

Nicholas Syring111Department of Statistics, Iowa State University;
<EMAIL_ADDRESS> and Ryan Martin222Department of Statistics, North
Carolina State University<EMAIL_ADDRESS>
###### Abstract

Bayesian posterior distributions are widely used for inference, but their
dependence on a statistical model creates some challenges. In particular,
there may be lots of nuisance parameters that require prior distributions and
posterior computations, plus a potentially serious risk of model
misspecification bias. Gibbs posterior distributions, on the other hand, offer
direct, principled, probabilistic inference on quantities of interest through
a loss function, not a model-based likelihood. Here we provide simple
sufficient conditions for establishing Gibbs posterior concentration rates
when the loss function is of a sub-exponential type. We apply these general
results in a range of practically relevant examples, including mean
regression, quantile regression, and sparse high-dimensional classification.
We also apply these techniques in an important problem in medical statistics,
namely, estimation of a personalized minimum clinically important difference.

_Keywords and phrases:_ classification; generalized Bayes; M-estimation;
minimum clinically important difference; model misspecification.

##  1 Introduction

A major selling point of the Bayesian framework is that it is normative: to
solve a new problem, one only needs a statistical model/likelihood, a prior
distribution for the parameters of that model, and the tools necessary for
computing the corresponding posterior distribution. While the Bayesians’
obligation to specify a prior attracts the most criticism, their need for a
model/likelihood has a number of potentially negative consequences too,
especially when the quantity of interest has meaning independent of the
statistical model, like a population quantile. On the one hand, even if the
posited model is “right,” it is rare that all the parameters of that model are
relevant. For example, if we are interested in a quantile of a distribution,
then the focus is only on a specific real-valued function of the full set of
parameters. In such a case, the non-trivial effort invested in dealing with
these nuisance parameters, e.g., specifying reasonable prior distributions and
designing efficient computational algorithms, is effectively wasted. On the
other hand, an even greater concern is that, in the far more likely case where
the posited model is “wrong,” that model misspecification can negatively
impact one’s conclusions about the quantity of interest. For example, both the
gamma and Pareto models have a 0.8 quantile, and neither impose any
constraints on the value that quantile can take, but the quality of inferences
drawn about that quantile may vary significantly depending on which of these
two models is chosen.

The non-negligible dependence on the posited statistical model puts a burden
on the data analyst, and those reluctant to take that risk tend to opt for a
non-Bayesian approach. After all, if one can get a solution without specifying
a statistical model, then it is impossible to incur model misspecification
bias. But in taking such an approach, they give up the normative advantage of
Bayesian analysis. Is it possible to get the best of both worlds? That is, can
one construct a posterior distribution for the quantity of interest directly,
incorporating any prior information that might be available about the quantity
of interest, without specifying a statistical model and incurring the
associated model misspecification risks, and without the need for
marginalization over nuisance parameters? Fortunately, the answer is Yes, and
this is the present paper’s focus.

The so-called _Gibbs posterior distribution_ is the proper prior-to-posterior
update when data and the interest parameter are linked by a loss function
rather than a likelihood. For a derivation of the Gibbs posterior
distribution, see Catoni, (2004), Zhang, (2006), and Bissiri et al., (2016),
who show the Gibbs posterior solves a meaningful risk minimization problem,
and satisfies a relevant coherence property. Intuitively, the Gibbs and
Bayesian posterior distributions coincide when the loss function linking data
and parameter is a (negative) log-likelihood. In that case the properties of
the Gibbs posterior can be inferred from the literature on Bayesian
asymptotics in both the well-specified and misspecified contexts. For cases in
which the link between data and interest parameter is not a likelihood, the
large-sample behavior of the Gibbs posterior is less clear and elucidating
this behavior under some simple and fairly general conditions is our goal
here.

As a practical example, medical investigators are often interested to
determine if a treatment whose effect has been judged to be statistically
significant is also clinically significant in the sense that the patients
actually feel better post-treatment. Therefore, they are interested in
estimating and quantifying uncertainty about the effect size cutoff beyond
which patients feel better; this is called the minimum clinically important
difference or MCID for short, e.g., Jaescheke et al., (1989). The data
consists of $(x,y)$ pairs, where $x$ is a continuous diagnostic measure, e.g.,
the patient’s estimated effect size, and $y$ is a binary indicator of whether
the patient felt better post-treatment or not. Estimation of the MCID boils
down to a classification problem and we have found that the standard Bayesian
approaches, e.g., assuming a logistic regression model and working with the
marginal posterior of the MCID, often do not perform well, due to its
sensitivity to model misspecification bias. Alternatively, one could consider
a more sophisticated version of binary regression, e.g., one where the link
function is modeled nonparametrically (Choudhuri et al.,, 2007), but this
seems like overkill when the quantity of interest is a scalar. We found that a
Gibbs posterior distribution, as described above, provided a very reasonable
and robust solution to the MCID problem; see Syring and Martin, (2017). There
are even more challenging versions of the MCID problem, ones that seek a
“personalized” or subject-specific cutoff that depends on a set of additional
covariates attached to subjects. In these cases, the MCID is high- or even
infinite-dimensional, and our previous Gibbs posterior analysis is not
equipped to handle such a complicated situation. A primary motivation for the
present work was to develop a Gibbs posterior distribution theory that could
handle these new MCID applications; see Section 6.

In the following sections we lay out and apply conditions under which a Gibbs
posterior distribution concentrates, asymptotically, on a neighborhood of the
true value of the inferential target as the sample sizes increases. Our focus
is not on the most general set of sufficient conditions for concentration;
rather, we aim for conditions that are both widely applicable and easily
verified. To that end we consider loss functions that satisfy an inequality
similar to the moment-generating function bound for sub-exponential random
variables (e.g., Boucheron et al.,, 2012). We can apply this condition in a
variety of problems from mean and quantile regression to classification, and
in both fixed- and high-dimensional settings. As an added advantage, our
conditions lead to short, straightforward proofs of concentration.

Section 2 provides some background information and formally defines the Gibbs
posterior distribution. In Section 3, we state our theoretical objectives and
present our main results, namely, sets of sufficient conditions under which
the Gibbs posterior achieves a specified asymptotic concentration rate. A
unique attribute of the Gibbs posterior distribution is its dependence on a
tuning parameter called the learning rate, and our results in this section
cover a constant learning rate, a vanishing sequence of learning rates, and
even random learning rate sequences. Section 4 further discusses verifying our
conditions and extends our conditions and main results to handle certain
unbounded loss functions. Section 5 applies our general theorems to establish
Gibbs posterior concentration rates in a number of practically relevant
examples, including quantile regression, estimation of the area under the
receiver operating characteristic curve, nonparametric curve estimation, and
high-dimensional sparse classification. Section 6 focuses on our motivating
example, namely, estimation of the MCID function, where we establish
concentration rates for our proposed Gibbs posterior distribution. Some
concluding remarks are given in Section 7, and the proofs of our technical
results are postponed to the appendices.

##  2 Background on Gibbs posteriors

###  2.1 Notation and definitions

Consider a measurable space $(\mathbb{U},\mathcal{U})$, with $\mathcal{U}$ a
$\sigma$-algebra of subsets of $\mathbb{U}$, on which a probability measure
$P$ is defined. A random element $U\sim P$ need not be a scalar, and many of
the applications we have in mind involve $U=(X,Y)$ or $U=(X,Y,Z)$, where $Y$
denotes a “response” variable and $X$ or $(X,Z)$ denotes a “predictor”
variable, and $P$ encodes the dependence between the entries in $U$. Then the
real-world phenomenon under investigation is determined by $P$ and our goal is
to make inference on a relevant feature of $P$, which we define as a given
functional $\theta=\theta(P)$, taking values in $\Theta$. Expressing $\theta$
as a functional of $P$ reveals that $\theta$ typically is not a parameter of
some statistical model, it is an inherent feature of the phenomenon under
investigation, something that exists and has real-world meaning beyond the
statistical model. Note that $\theta$ could be finite-, high-, or even
infinite-dimensional, and examples of all three types are investigated here.

The specific way $\theta$ relates to $P$ guides our posterior construction.
Suppose there is a loss function, $\ell_{\theta}(u)$, that measures how
closely a generic value of $\theta$ agrees with a data point $u$. (As is
customary, “$\theta$” will denote both the quantity of interest and a generic
value of that quantity; when we need to distinguish the true from a generic
value, we will write “$\theta^{\star}$.”) For example, if $u=(x,y)$ is a
predictor–response pair, and $\theta$ is a function, then the loss might be

$\ell_{\theta}(u)=|y-\theta(x)|\quad\text{or}\quad\ell_{\theta}(u)=1\\{y\neq\theta(x)\\},$
(1)

depending on whether $y$ is continuous or discrete/binary, where $1(A)$
denotes the indicator function at the event $A$. Another common situation is
when one specifies a statistical model, say,
$\\{P_{\theta}:\theta\in\Theta\\}$, indexed by a parameter $\theta$, and sets
$\ell_{\theta}(u)=-\log p_{\theta}(u)$, where $p_{\theta}$ is the density of
$P_{\theta}$ with respect to some fixed dominating measure. In all of these
cases, the idea is that a loss is incurred when there is a certain discrepancy
between $\theta$ and the data point $u$. Then our inferential target is the
value of $\theta$ that minimizes the risk or average loss/discrepancy.

######  Definition 1.

Consider a real-valued loss function $\ell_{\theta}(u)$ defined on
$\mathbb{U}\times\Theta$, and define the risk function
$R(\theta)=P\ell_{\theta}$, the expected loss with respect to $P$; throughout,
we write $Pf$ to denote expected value of $f(U)$ with respect to $U\sim P$.
Then the inferential target is

$\theta^{\star}\in\arg\min_{\theta\in\Theta}R(\theta).$ (2)

Our focus here will be on case where the risk function minimizer,
$\theta^{\star}$, is unique, so that the “$\in$” in (2) is an equality. Given
that estimation/inference is our goal, it is only natural to assume that the
inferential target is identifiable. But identifiability is not a trivial
assumption so, fortunately, is not absolutely necessary for our theory.
Indeed, the main results in Section 3 remain valid even if the risk minimizer
is not unique, and we make a few brief comments about this extension in the
discussion following Theorem 1.

The risk function itself is unavailable—because it depends on the unknown
$P$—and, therefore, so is $\theta^{\star}$. However, suppose that we have an
independent and identically distributed (iid) sample
$U^{n}=(U_{1},\ldots,U_{n})$ of size $n$, with each $U_{i}$ having marginal
distribution $P$ on $\mathbb{U}$. The iid assumption is not crucial, but it
makes the notation and discussion easier; an extension to independent but not
identically distributed (inid) cases is discussed in the context of an
important example in Section 5.4. In general, we have data $U^{n}$ taking
values in the measurable space $(\mathbb{U}^{n},\mathcal{U}^{n})$, with joint
distribution denoted by $P^{n}$. From here, we can naturally replace the
unknown risk function in Definition 1 with an empirical version and proceed
accordingly.

######  Definition 2.

For a loss function $\ell_{\theta}$ as described above, define the empirical
risk as

$R_{n}(\theta)=\mathbb{P}_{n}\ell_{\theta}=n^{-1}\sum_{i=1}^{n}\ell_{\theta}(U_{i}),$
(3)

where $\mathbb{P}_{n}=n^{-1}\sum_{i=1}^{n}\delta_{U_{i}}$, with $\delta_{u}$
the Dirac point-mass measure at $u$, is the empirical distribution based on
the sample $U^{n}$.

Naturally, if the inferential target is the risk minimizer, then it makes
sense to estimate that quantity based on data $U^{n}$ by minimizing the
empirical risk, i.e.,

$\hat{\theta}_{n}\in\arg\min_{\theta\in\Theta}R_{n}(\theta).$ (4)

This is precisely the M-estimator based on an objective function determined by
the loss $\ell_{\theta}$; when $R_{n}$ is differentiable, the root of
$\dot{R}_{n}$, the derivative of $R_{n}$, is a Z-estimator and
“$\dot{R}_{n}(\theta)=0$” is often called an estimating equation (Godambe,,
1991; van der Vaart,, 1998). Since $\theta\mapsto\ell_{\theta}$ need not be
smooth or convex, and $R_{n}$ is an average over a finite set of data, it is
possible that its minimizer is not unique, even if $\theta^{\star}$ is. These
computational challenges are, in fact, part of the motivation for considering
Gibbs posteriors in the first place, as we discuss below.

There is a rich literature on the asymptotic distribution properties of
M-estimators, which can be used for developing hypothesis tests and confidence
intervals (Maronna et al.,, 2006; Huber and Ronchetti,, 2009). As an
alternative, one might consider a Bayesian approach to quantify uncertainty,
but there is an immediate obstacle, namely, no statistical model/likelihood
connecting the data to the quantity of interest. If we did have a statistical
model, with a density function $p_{\theta}$, then the most natural loss is
$\ell_{\theta}(u)=-\log p_{\theta}(u)$ and the likelihood is
$\exp\\{-nR_{n}(\theta)\\}$. It is, therefore, tempting to follow that same
strategy for general losses, resulting in a sort of generalized posterior
distribution for $\theta$.

######  Definition 3.

Given a loss function $\ell_{\theta}$ and the corresponding empirical risk
$R_{n}$ in Definition 2, define the Gibbs posterior distribution as

$\Pi_{n}(d\theta)\propto
e^{-\omega\,nR_{n}(\theta)}\,\Pi(d\theta),\quad\theta\in\Theta,$ (5)

where $\Pi$ is a prior distribution on $\Theta$ and $\omega>0$ is a so-called
learning rate parameter (Holmes and Walker,, 2017; Syring and Martin,, 2019;
Grünwald,, 2012; van Erven et al.,, 2015), discussed more below. For
simplicity, the dependence of $\Pi_{n}$ on $\omega$ will generally be omitted
from the notation, but see Section 3.3.

The name is derived from its connection to the so-called Gibbs measures from
statistical physics. In that setting, our empirical risk is called the energy
function and the learning rate is the inverse temperature. We have and will
continue to assume that the right-hand side of (5) is integrable in $\theta$,
so that the proportionality constant is well-defined and equal to the
reciprocal of that integral. Integrability holds whenever the loss function is
bounded from below, like for those in (1), but this could fail in some cases
where the loss is not bounded away from $-\infty$, as can happen when
$\ell_{\theta}(u)$ is a negative log-density. In such cases, extra conditions
on the prior distribution would be required to ensure the Gibbs posterior is
well-defined.

An immediate advantage of this approach, compared to the M-estimation strategy
described above, is that the user is able to incorporate available prior
information about $\theta$ directly into the analysis. This is especially
important in cases where the quantity of interest has a real-world meaning or
interpretation, as opposed to being just a model parameter, where having
genuine prior information is the norm rather than the exception. Additionally,
even though there is no likelihood, the same computational methods, such as
Markov chain Monte Carlo (e.g., Chernozhukov and Hong,, 2003) and variational
approximations (e.g., Alquier et al.,, 2016), common in Bayesian analysis, can
be employed to numerically approximate the Gibbs posterior.

We have opted here to define the Gibbs posterior directly as an object to be
used and studied, but there is a more formal, more principled way in which
Gibbs posteriors emerge. In the PAC-Bayes literature, the goal is to construct
a randomized estimator—basically, a posterior distribution—that concentrates
in regions of $\Theta$ where the risk, $R(\theta)$, or its empirical version,
$R_{n}(\theta)$, is small (e.g., Valiant,, 1984; McAllester,, 1999; Alquier,,
2008; Guedj,, 2019). That is, the Gibbs posterior distribution can be viewed
as a solution to an optimization problem rather than a solution to the
updating-prior-beliefs problem. More formally, for a given prior $\Pi$ on
$\Theta$, suppose the goal is to find

$\inf_{\mu}\Bigl{\\{}\int R_{n}(\theta)\,\mu(d\theta)+(\omega
n)^{-1}K(\mu,\Pi)\Bigr{\\}},$ (6)

where the infimum is over all probability measures $\mu$ that are absolutely
continuous with respect to $\Pi$, and $K$ denotes the Kullback–Leibler
divergence. Then it can be shown that the unique solution is $\Pi_{n}$, the
Gibbs posterior defined in (5); see, Zhang, (2006) and Bissiri et al., (2016).
Therefore, the Gibbs posterior distribution is the measure minimizing a
penalized risk, averaged with respect to a given prior, $\Pi$.

###  2.2 Learning rate

Readers familiar with M-estimation may not recognize the learning rate,
$\omega$. This does not appear in the M-estimation context because all that
influences the optimization problem—and the corresponding asymptotic
distribution theory—is the shape of the loss/risk function, not its magnitude
or scale. On the other hand, the learning rate is an essential component of
the Gibbs posterior distribution in (5) since the distribution depends on both
the shape and scale of the loss function.

The data analyst must choose a learning rate, but there is no obviously
optimal way to do so. A Bayesian would be tempted to treat the learning rate
as an unknown parameter and assign it a prior distribution. However, since the
data enters the model only through the loss function, and since the loss
minimizer is insensitive to the learning rate, it is not clear that a prior-
to-posterior update will “learn” the “right” value of the learning rate.
Rather than viewing the learning rate as a model parameter we view it as a
tuning parameter similar to the role of the penalty parameter in ridge
regression. In penalized likelihood estimation, the penalty parameter controls
the strength of the likelihood relative to the penalty, and the learning rate
plays the same role for the Gibbs posterior by weighting the influence of data
versus the prior distribution. So, even though it is not a model parameter, it
is reasonable to use data-driven methods to select the learning rate similar
to how cross-validation is used to select penalty parameters. Data-driven
strategies for tuning the learning rate are available; see Holmes and Walker,
(2017), Syring and Martin, (2019), Lyddon et al., (2019), and Wu and Martin,
(2020, 2021).

In this paper we focus on how the learning rate affects posterior
concentration. In most cases, our proofs require the learning rate is no
larger than a constant upper bound. The upper bound depends on features of
$P$, which are generally unknown, so, in practice, the learning rate can be
taken to be a slowly vanishing sequence, which has a negligible effect on the
concentration rate. In more challenging examples, we require the learning rate
to vanish sufficiently quickly.

###  2.3 Relation to other generalized posterior distributions

We refer to any posterior distribution besides the well-specified Bayesian
posterior distribution as a generalized posterior distribution, with the most
familiar example being a misspecified Bayesian posterior. Readers may be
tempted to interpret Gibbs posterior distributions as such; however, these two
kinds of distributions are actually very different. Generally there is no
genuine likelihood function that links the data to the risk minimizer
$\theta$. Consequently, any model-based approach to inference on $\theta$ must
be indirect in the following sense. The data analyst will do their best to
identify a model that “fits” the data, but the empirical risk function will
rarely be among those models under consideration because it will not equal any
familiar negative log-likelihood. Therefore, the selected model is indexed by
some other parameter, say, $\gamma$, and, after specifying a prior for
$\gamma$ and then updating to a posterior, the posterior distribution for
$\theta=\theta(P_{\gamma})$ is obtained through marginalization. Ideally, the
posterior for $\gamma$ will concentrate around a point $\gamma^{\dagger}$ such
that the Kullback–Leibler divergence of $P_{\gamma}$ from $P$ is minimized
(e.g., Kleijn and van der Vaart,, 2006; De Blasi and Walker,, 2013;
Ramamoorthi et al.,, 2015). But there is no guarantee that
$\theta(P_{\gamma^{\dagger}})=\theta^{\star}$, hence the risk of
misspecification bias, not to mention the potential dangers associated with
marginalization (e.g., Fraser,, 2011; Martin,, 2019). The point is that
Bayesian posteriors are only misspecified accidentally, so any direct
posterior distribution for a risk minimizer based on an intentionally
misspecified model cannot claim to be Bayesian—it is a Gibbs posterior. Beyond
the semantics, there are practical consequences to mistaking Gibbs posteriors
as Bayesian: for example, if the context were Bayesian, then one would either
ignore the learning rate entirely or be tempted to treat the learning rate in
a Bayesian way, which does not work.

Treating whatever might be interpreted as “model misspecification” in a Gibbs
formulation as intentional provides some additional flexibility. Even in cases
where the quantity of interest is not characterized by an optimization problem
as in Definition 1, we might consider “reverse engineering” a suitable
loss/risk function such that the given quantity of interest solves the
corresponding optimization problem; see, e.g., Syring and Martin, (2020) and
Wang and Martin, (2020). And, even if a statistical model is within reach, one
may choose not to use it, e.g., to avoid dealing with nuisance parameters or
because the likelihood is computationally intractable. The point is that Gibbs
posteriors can be employed in a rather broad range of problems as an
alternative to the Bayesian formulation that insists on a statistical model.

A number of authors have studied generalized posterior distributions formed
using a “power likelihood”—essentially a Bayesian posterior in which the
likelihood is raised to a fractional power $\eta$ between $0$ and $1$; see,
e.g., Martin et al., (2017) and Grünwald and van Ommen, (2017) among others.
These $\eta$-generalized posteriors tend to be robust to misspecification of
the probability model, and data-driven methods to optimize $\eta$ are
developed in Grünwald and van Ommen, (2017). These posteriors coincide with
Gibbs posteriors based on a log-loss and with the learning rate $\omega$
corresponding to $\eta$.

Another variety of generalized posterior distribution comes from moment-based
estimation in which $P$ is not fully specified but is assumed to obey a number
of moment conditions of the form $Pf_{j}=0$ for functions $f_{1},\ldots
f_{J}$; and see Kim, (2002), Chernozhukov and Hong, (2003), and Chib et al.,
(2018). Kim, (2002) uses the moment conditions to construct a so-called
limited-information likelihood to use in place of a fully-specified likelihood
in a Bayesian formulation. Similarly, Chernozhukov and Hong, (2003) construct
a posterior by taking a (pseudo) log-likelihood equal to a quadratic form
determined by the set of moment conditions. Chib et al., (2018) studies
Bayesian exponentially-tilted empirical likelihood posterior distributions,
also defined by moment conditions. In some cases the Gibbs posterior
distribution coincides with the above moment-based methods. For instance, in
the special case that the risk $R$ is differentiable at $\theta^{\star}$ with
derivative $\dot{R}(\theta^{\star})$, risk-minimization corresponds to the
single moment condition $\dot{R}(\theta^{\star})=0$.

##  3 Asymptotic concentration rates

###  3.1 Objective

Important questions in the Bayesian literature concern the asymptotic
concentration properties of their posterior distributions. Roughly, if data
are generated from a distribution for which the quantity of interest takes
value $\theta^{\star}$, then, as the sample size increases, the posterior
distribution ought to concentrate its mass around that same $\theta^{\star}$.
As Shalizi, (2009) puts it, such concentration results offer a “frequentist
license for Bayesian practice.” Our goal here is to develop concentration rate
results for Gibbs posteriors that parallel those in the Bayesian literature.
As we will show, optimal concentration rates are possible with Gibbs
posteriors, so the robustness achieved by not specifying a statistical model
has no cost in terms of (asymptotic) efficiency.

Towards a precise statement of our objectives, for a fixed
$\theta^{\star}\in\Theta$, let $d(\theta;\theta^{\star})$ denote a divergence
measure on $\Theta$ in the sense that $d(\theta;\theta^{\star})\geq 0$ for all
$\theta$, with equality if and only if $\theta=\theta^{\star}$. It could, of
course, be a proper metric on $\Theta$ but we will not insist here on symmetry
or the triangle inequality. The divergence measure could also depend on the
sample size $n$ or other deterministic features of the problem at hand,
especially in the independent but not iid setting; see Section 5.4.

Our objective is to provide conditions under which the Gibbs posterior will
concentrate asymptotically, at a certain rate, around $\theta^{\star}$
relative to the divergence measure $d$. Throughout this paper,
$(\varepsilon_{n})$ denotes a deterministic sequence of positive numbers with
$\varepsilon_{n}\to 0$, which will be referred to as the Gibbs posterior
concentration rate.

######  Definition 4.

The Gibbs posterior $\Pi_{n}$ in (5) asymptotically concentrates around
$\theta^{\star}$ at rate (at least) $\varepsilon_{n}$, with respect to $d$, if

$P^{n}\Pi_{n}(\\{\theta:d(\theta;\theta^{\star})>M_{n}\varepsilon_{n}\\})\rightarrow
0,\quad\text{as $n\to\infty$,}$ (7)

where $M_{n}>0$ is either a (deterministic) sequence satisfying
$M_{n}\to\infty$ arbitrarily slowly or is a sufficiently large constant,
$M_{n}\equiv M$.

In the PAC-Bayes literature, the Gibbs posterior distribution is interpreted
as a “randomized estimator,” a generator of random $\theta$ values that tend
to make the risk difference small. For iid data and with risk divergence
$d(\theta;\theta^{\star})=\\{R(\theta)-R(\theta^{\star})\\}^{1/2}$, a
concentration result like that in Definition 4 makes this strategy clear since
the $\Pi_{n}$-probability of the event
$\\{\theta:R(\theta)-R(\theta^{\star})\leq\varepsilon_{n}^{2}\\}$ would be
asymptotically close to 1.

If the Gibbs posterior concentrates around $\theta^{\star}$ in the sense of
Definition 4, then any reasonable estimator derived from that distribution,
such as the mean, should inherit the $\varepsilon_{n}$ rate at
$\theta^{\star}$ relative to the divergence measure $d$. This can be made
formal under certain conditions on $d$; see, e.g., Corollary 1 in Barron et
al., (1999) and the discussion following the proof of Theorem 2.5 in Ghosal et
al., (2000).

Besides concentration rates, in certain cases it is possible to establish
distributional approximations to Gibbs posteriors, i.e., Bernstein-von Mises
theorems. Results for finite-dimensional problems and with sufficiently smooth
loss functions can be found in, e.g., Bhattacharya and Martin, (2020) and
Chernozhukov and Hong, (2003).

###  3.2 Conditions

Here we discuss a general strategy for proving Gibbs posterior concentration
and the kinds of sufficient conditions needed for the strategy to be
successful. To start, set
$A_{n}=\\{\theta:d(\theta,\theta^{\star})>M_{n}\varepsilon_{n}\\}\subset\Theta$.
Then our first step towards proving concentration is to express
$\Pi_{n}(A_{n})$ as the ratio

$\Pi_{n}(A_{n})=\frac{N_{n}(A_{n})}{D_{n}}=\frac{\int_{A_{n}}\exp[-\omega
n\\{R_{n}(\theta)-R_{n}(\theta^{\star})\\}]\,\Pi(d\theta)}{\int_{\Theta}\exp[-\omega
n\\{R_{n}(\theta)-R_{n}(\theta^{\star})\\}]\,\Pi(d\theta)}.$ (8)

The goal is to suitably upper and lower bound $N_{n}(A_{n})$ and $D_{n}$,
respectively, in such a way that the ratios of these bounds vanish.

Two sufficient conditions are discussed below, the first primarily dealing
with the loss function and aiding in bounding $N_{n}(A_{n})$ and the second
primarily concerning the prior distribution and aiding in bounding $D_{n}$.
Both conditions concern the excess loss
$\ell_{\theta}(U)-\ell_{\theta^{\star}}(U)$ and its mean and variance:

$m(\theta,\theta^{\star}):=P(\ell_{\theta}-\ell_{\theta^{\star}})\quad\text{and}\quad
v(\theta,\theta^{\star}):=P(\ell_{\theta}-\ell_{\theta^{\star}})^{2}-m(\theta,\theta^{\star})^{2}.$

####  3.2.1 Sub-exponential type losses

Our method for proving concentration of the Gibbs posterior distribution
requires a vanishing upper bound on the expectation of the numerator term
$N_{n}(A_{n})$ in (8). Since the integrand, $\exp[-\omega
n\\{R_{n}(\theta)-R_{n}(\theta^{\star})\\}]$, in $N_{n}(A_{n})$ is non-
negative, Fubini’s theorem says it suffices to bound the expected integrand.
Further, independence of $U^{n}$ implies

$P^{n}e^{-\omega
n\\{R_{n}(\theta)-R_{n}(\theta^{\star})\\}}=\\{Pe^{-\omega(\ell_{\theta}-\ell_{\theta^{\star}})}\\}^{n},$

which reveals that the key to bounding $P^{n}N_{n}(A_{n})$ is to bound
$Pe^{-\omega(\ell_{\theta}-\ell_{\theta^{\star}})}$, the expected
exponentiated excess loss. In order for the bound on $P^{n}N_{n}(A_{n})$ to
vanish it must be that $Pe^{-\omega(\ell_{\theta}-\ell_{\theta^{\star}})}<1$,
but this is not enough to identify the concentration rate in (7). Rather, the
speed at which $Pe^{-\omega(\ell_{\theta}-\ell_{\theta^{\star}})}$ vanishes
must be a function of $d(\theta,\theta^{\star})$, and we take this
relationship as our key condition for Gibbs posterior concentration. When the
condition holds we say the loss function is of sub-exponential type.

######  Condition 1.

There exists an interval $(0,\bar{\omega})$ and constants $K,\,r>0$, such that
for all $\omega\in(0,\bar{\omega})$ and for all sufficiently small $\delta>0$,
for $\theta\in\Theta$

$d(\theta;\theta^{\star})>\delta\implies
Pe^{-\omega(\ell_{\theta}-\ell_{\theta^{\star}})}<e^{-K\omega\delta^{r}}$ (9)

Following the above arguments, Condition 1 and the definition of $A_{n}$ imply
the finite-sample bound

$P^{n}N_{n}(A_{n})\leq e^{-K\omega nM_{n}^{r}\varepsilon_{n}^{r}}.$ (10)

The constant $r>0$ that appears in Condition 1 and our various bounds on the
prior distribution (discussed below) can take on different values depending on
the context. The most familiar is $r=2$, but some “non-regular” problems
require different values of $r$.

Condition 1 is closely related to the developments in Grünwald and Mehta,
(2020). In particular, they focus on a so-called annealed expectation which,
for a real-valued function $f$ of the random element $U\sim P$ and a constant
$\omega>0$, is defined as

$P^{\text{\sc ann}(\omega)}f=-\omega^{-1}\log Pe^{-\omega f}.$

With this definition, it is easy to verify that

$Pe^{-\omega(\ell_{\theta}-\ell_{\theta^{\star}})}=\exp\\{-\omega P^{\text{\sc
ann}(\omega)}(\ell_{\theta}-\ell_{\theta^{\star}})\\},$

so a suitable upper bound on the left-hand side, as we require in our
Condition 1, is equivalent to a corresponding lower bound on $P^{\text{\sc
ann}(\omega)}(\ell_{\theta}-\ell_{\theta^{\star}})$. The “strong central
condition” in Grünwald and Mehta, (2020) states that

$P^{\text{\sc ann}(\omega)}(\ell_{\theta}-\ell_{\theta^{\star}})\geq
0,\quad\text{for all $\theta$},$

which implies $Pe^{-\omega(\ell_{\theta}-\ell_{\theta^{\star}})}\leq 1$ and,
in turn, that the expectation of $N_{n}(A_{n})$ is bounded as $n\to\infty$.
The other conditions in Grünwald and Mehta, (2020) aim to lower-bound the
annealed expected loss difference by a suitable function of the risk
difference. For example, their Lemma 13 shows that a “witness condition”
implies

$P^{\text{\sc ann}(\omega)}(\ell_{\theta}-\ell_{\theta^{\star}})\gtrsim
R(\theta)-R(\theta^{\star}),$

so, if $R(\theta)-R(\theta^{\star})$ is of the order
$d(\theta,\theta^{\star})^{r}$ for some $r>0$, then we recover Condition 1.
Therefore, in a certain sense, our Condition 1 is exactly what is needed to
suitably control the numerator of the Gibbs posterior, and the relevant
conditions developed by Grünwald and Mehta, (2020) and others are designed to
imply our Condition 1.

Section 3.4 provides more details on verifying Condition 1, and the strategies
presented there all involve lower and upper bounds on
$m(\theta,\theta^{\star})$ and $v(\theta,\theta^{\star})$ in terms of
$d(\theta,\theta^{\star})$. For now, consider the special case of a well-
specified statistical model, i.e., a collection of densities
$\\{p_{\theta}:\theta\in\Theta\\}$ containing the “true” $\theta^{\star}$
satisfying $P=P_{\theta^{\star}}$. In this setting Condition 1 follows from
properties of the Rényi divergence measure; and, see van Erven and Harremoës,
(2014). If $\omega\in(0,1)$, then the expected exponentiated loss difference
is

$P_{\theta^{\star}}e^{-\omega(\ell_{\theta}-\ell_{\theta^{\star}})}=\int
p_{\theta}^{\omega}p_{\theta^{\star}}^{1-\omega}\,d\mu.$ (11)

In particular, van Erven and Harremoës, (2014) show

$P_{\theta^{\star}}e^{-\omega(\ell_{\theta}-\ell_{\theta^{\star}})}\leq\exp\bigl{\\{}-\tfrac{\omega(1-\omega)}{2}d_{\text{\sc
tv}}^{2}(p_{\theta},p_{\theta^{\star}})\bigr{\\}},$ (12)

where $d_{\text{\sc tv}}$ denotes the total variation distance on densities.
If separation between $\theta$ and $\theta^{\star}$ with respect to the
divergence measure $d$ implies separation in terms of total variation on
densities, then Condition 1 holds.

Besides Condition 1 other strategies for bounding $N_{n}(A_{n})$ in (8) have
been employed in the literature. At least for finite-dimensional problems like
those in Syring and Martin, (2017) and Bhattacharya and Martin, (2020),
empirical process techniques have proved useful. Generally, these require a
uniform law of large numbers, i.e., $\sup_{\theta}|R_{n}(\theta)-R(\theta)|\to
0$ in probability as $n\to\infty$, which can be challenging to verify in non-
parametric problems. Chernozhukov and Hong, (2003) require the stronger
condition that
$\sup_{d(\theta,\theta^{\star})>\delta}\\{R_{n}(\theta)-R_{n}(\theta^{\star})\\}\rightarrow
0$ in $P^{n}$-probability. When it holds they immediately recover an in-
probability bound on $N_{n}(A_{n})$ irrespective of $\omega$, but do not
obtain a finite-sample bound on $P^{n}\Pi_{n}(A_{n})$. Their analysis is
limited to finite-dimensional parameters and empirical risk functions
$R_{n}(\theta)$ that are smooth in a neighborhood of $\theta^{\star}$. This
latter condition excludes important examples like the misclassification-error
loss function used in classification (see Section 5.5), but facilitates
proving distributional convergence. The challenges in applying empirical
process theory to non-parametric problems motivated us to consider alternative
conditions, and led us to applying a version of Condition 1 in the context of
a non-parametric image segmentation problem (Syring and Martin,, 2020). Our
goal in this paper is to generalize and extend the approach used in that
special case.

There are situations in which Condition 1 does not hold, and in Section 4 we
address some of these important cases with suitable generalizations and
extensions of Condition 1. As an example, note that both Condition 1 and the
witness condition used in Grünwald and Mehta, (2020) are closely related to a
Bernstein-type condition relating the first two moments of the excess loss:
$P(\ell_{\theta}-\ell_{\theta^{\star}})^{2}\leq
c\\{R(\theta)-R(\theta^{\star})\\}^{\alpha}$ for constants $(c,\,\alpha)>0$.
For the loss function used in quantile estimation, the Bernstein condition is
generally not satisfied and, consequently, neither the witness condition nor
Condition 1 can be verified. Therefore, a new/modified approach is needed to
handle quantiles, and we present one such approach in Section 4.1 with
applications in Sections 5.1 and 5.6

####  3.2.2 Prior distribution

Generally, the prior distribution must place enough mass on certain “risk-
function” neighborhoods $G_{n}:=\\{\theta:m(\theta,\theta^{\star})\vee
v(\theta,\theta^{\star})\leq\varepsilon_{n}^{r}\\}$. This is analogous to the
requirement in the Bayesian literature that the prior place sufficient mass on
Kullback–Leibler neighborhoods of $\theta^{\star}$; see, e.g., Shen and
Wasserman, (2001, Theorem 2) and Ghosal et al., (2000, Equation 2.4). Some
version of the following prior bound is needed

$\log\Pi(G_{n})\gtrsim-n\varepsilon_{n}^{r},$ (13)

for $r$ as in Condition 1. Lemma 1 in Appendix A.1 along with (13) provides a
lower probability bound on the denominator term $D_{n}$ defined in (8). That
is,

$P^{n}\bigl{(}D_{n}\leq\tfrac{1}{2}\Pi(G_{n})e^{-2\omega
n\varepsilon_{n}^{r}}\bigr{)}\leq
P\left(D_{n}\leq\tfrac{1}{2}e^{-(C+2\omega)n\varepsilon_{n}^{r}}\right)\leq(n\varepsilon_{n}^{r})^{-1},$
(14)

where $b_{n}=\frac{1}{2}\Pi(G_{n})e^{-2\omega n\varepsilon_{n}^{r}}$, which,
according to (13), is of the order $e^{-C\omega n\varepsilon_{n}^{r}}$, for an
appropriate constant $C>0$. Both the form of the risk-function neighborhoods
and the precise lower bound in (13) depend on the concentration rate and the
learning rate, so the results in Section 3.3 all require their own version of
the above prior bound.

Notably, Grünwald and Mehta, (2020) only require bounds on the prior mass of
the simpler neighborhoods
$\\{\theta:m(\theta,\theta^{\star})\leq\varepsilon_{n}^{r}\\}$ without the
second moment term, $v$. Under their condition we can derive a lower bound on
$P^{n}D_{n}$ similar to Lemma 1 in Martin et al., (2013). However, our proofs
require an in-probability lower bound on $D_{n}$, which requires the stronger
prior concentration bound like that in (13). While there are important
examples where the lower bounds on $m(\theta,\theta^{\star})$ and
$v(\theta,\theta^{\star})$ are of different orders, none of the applications
we consider here have upper bounds of different orders. Therefore, the
stronger prior concentration condition in (13) does not cost us in terms of
the concentration rates we derive for the examples in Section 5. Moreover, as
discussed following the statement of Theorem 1 below, our finite-sample bounds
on the Gibbs posterior concentration are better compared to those in Grünwald
and Mehta, (2020), a direct consequence of our method of proof that requires
stronger prior concentration requirements in terms of “$(m\vee
v)$-neighborhoods” rather than Grünwald and Mehta’s weaker requirements in
terms of “$m$-neighborhoods.”

###  3.3 Main results

In this section we present general results on Gibbs posterior concentration.
The detailed proofs can be found in Appendix A. Our first result establishes
Gibbs posterior concentration under Condition 1 and a local prior positivity
condition for all sufficiently small, constant learning rates.

######  Theorem 1.

Let $\varepsilon_{n}$ be a vanishing sequence satisfying
$n\varepsilon_{n}^{r}\to\infty$ for a constant $r>0$. Suppose the prior
satisfies

$\log\Pi(\\{\theta:m(\theta,\theta^{\star})\vee
v(\theta,\theta^{\star})\leq\varepsilon_{n}^{r}\\})\geq-
C_{1}n\varepsilon_{n}^{r},$ (15)

for $C_{1}>0$ and for divergence measure $d$, the same $r>0$ as above, and
learning rate $\omega\in(0,\bar{\omega})$ for some $\bar{\omega}>0$. If the
loss function satisfies Condition 1, then the Gibbs posterior distribution in
(5) has asymptotic concentration rate $\varepsilon_{n}$, with $M_{n}\equiv M$
a large constant as in Definition 4.

For a brief sketch of the proof bound the posterior probability of $A_{n}$ by

$\displaystyle\Pi_{n}(A_{n})$
$\displaystyle\leq\frac{N_{n}(A_{n})}{D_{n}}\,1(D_{n}>b_{n})+1(D_{n}\leq
b_{n})$ $\displaystyle\leq b_{n}^{-1}N_{n}(A_{n})+1(D_{n}\leq b_{n}),$

where $b_{n}$ is as above, given explicitly in the proof Lemma 1. Taking
expectation of both sides, and applying (14) and the consequence (10) of
Condition 1, we get

$P^{n}\Pi_{n}(A_{n})\leq b_{n}^{-1}e^{-K\omega
nM_{n}^{r}\varepsilon_{n}^{r}}+(n\varepsilon_{n}^{r})^{-1}.$

Then the right-hand side is itself upper-bounded by a term of the order

$\max\\{e^{-K^{\prime}\omega
nM_{n}^{r}\varepsilon_{n}^{r}},(n\varepsilon_{n}^{r})^{-1}\\},$

which is generally the of the order $(n\varepsilon_{n}^{r})^{-1}$.

In contrast, Grünwald and Mehta, (2020) begin with Markov’s inequality,
bounding

$\Pi_{n}(A_{n})\leq(M_{n}\varepsilon_{n})^{-1}\Pi_{n}m(\cdot,\theta^{\star}),$

where they (almost) always take
$d(\theta,\theta^{\star})=m(\theta,\theta^{\star})$. They bound the posterior
expected risk on the right hand side using a PAC-Bayes inequality, their
(weaker) prior concentration condition, and other conditions like their
witness condition, depending on the context. In many cases, they achieve the
bound $\Pi_{n}m(\cdot,\theta^{\star})\lesssim\varepsilon_{n}$ on the posterior
expectation of the $m$ function so that $P^{n}\Pi_{n}(A_{n})\leq M_{n}^{-1}$,
which vanishes arbitrarily slowly.

Our stronger control over prior concentration relative to the condition in
Grünwald and Mehta, (2020) appears to be necessary for the in-probability
bound on $D_{n}$. It is this bound that enables our different method of proof
compared to
Grünwald and Mehta, (2020), and, ultimately, our improved finite sample bound.

For the case where the inferential target, the risk minimizer in (2), is not
unique, certain modifications of the above argument can be made, similar to
those in
Kleijn and van der Vaart, (2006), Theorem 2.4. Roughly, to our Theorem 1
above, we would add the requirements that (15) and Condition 1 hold uniformly
in $\theta^{\star}\in\Theta^{\star}$, where $\Theta^{\star}$ is the set of
risk minimizers. Then virtually the same proof strategy could be used to show
that

$\mathsf{E}\,\Pi_{n}(\\{\theta:d(\theta,\Theta^{\star})\gtrsim\varepsilon_{n}\\})\to
0,\quad n\to\infty,$

where
$d(\theta,\Theta^{\star})=\inf_{\theta^{\star}\in\Theta^{\star}}d(\theta,\theta^{\star})$
is the distance of $\theta$ from the set $\Theta^{\star}$.

Theorem 1 is quite flexible and can be applied in a range of settings; see
Section 5. However, one case in which it cannot be directly applied is when
$n\varepsilon_{n}^{r}$ is bounded. For example, in sufficiently smooth finite-
dimensional problems, we have $r=2$ and the target rate is
$\varepsilon_{n}=n^{-1/2}$. The difficulty is caused by the prior bound in
(15), since it is impossible—at least with a fixed prior—to assign mass
bounded away from 0 to a shrinking neighborhood of $\theta^{\star}$. One
option is to add a logarithmic factor to the rate, i.e., take
$\varepsilon_{n}=(\log n)^{k}n^{-1/2}$, so that $e^{-Cn\varepsilon_{n}^{2}}$
is a power of $n^{-1/2}$. Alternatively, we obtain the following result by
refining the proof of Theorem 1 and avoid slowing down the rate by a
logarithmic factor. Whenever our Theorem 2 holds, we achieve faster rates than
in Grünwald and Mehta, (2020) since their rate always includes a logarithmic
factor.

######  Theorem 2.

Consider a finite-dimensional $\theta$, taking values in
$\Theta\subseteq\mathbb{R}^{q}$ for some $q\geq 1$. Suppose that the target
rate $\varepsilon_{n}$ is such that $n\varepsilon_{n}^{r}$ is bounded for some
constant $r>0$. If the prior $\Pi$ satisfies

$\Pi(\\{\theta:m(\theta,\theta^{\star})\vee
v(\theta,\theta^{\star})\leq\varepsilon_{n}^{2}\\})\gtrsim\varepsilon_{n}^{q},$
(16)

and if Condition 1 holds, then the Gibbs posterior distribution in (5), with
any learning rate $\omega\in(0,\bar{\omega})$, has asymptotic concentration
rate $\varepsilon_{n}$ at $\theta^{\star}$ with respect to any divergence
$d(\theta,\theta^{\star})$ satisfying $\|\theta-\theta^{\star}\|\lesssim
d(\theta,\theta^{\star})\lesssim\|\theta-\theta^{\star}\|$ and for any
diverging, positive sequence $M_{n}$ in Definition 4.

The learning rate is critical to the Gibbs posterior’s performance, both
asymptotically and in finite samples, but it is challenging to determine the
upper bound $\bar{\omega}$ in examples. For instance, suppose the excess loss
$\ell_{\theta}-\ell_{\theta}^{\star}$ is normally distributed with variance
$\sigma^{2}(\theta)$ and obeys the Bernstein condition $\sigma^{2}(\theta)\leq
c\\{R(\theta)-R(\theta^{\star})\\}$ with $\alpha=1$. Normality is unimportant;
it is only used here to simplify this illustration. The Bernstein condition,
however, is important for verifying Condition 1 in most examples. The form of
the moment-generating function of the excess loss implies Condition 1 holds
for $d(\theta,\theta^{\star})=\\{R(\theta)-R(\theta^{\star})\\}^{1/2}$ with
$r=2$ if $\omega<2c^{-1}$. In practice, we can usually determine the exponent
$\alpha$ in the Bernstein condition, but we can rarely quantify the factor $c$
beyond simply showing the Bernstein condition holds for some $c>0$. Therefore,
we need the learning rate to be sufficiently small, but the meaning of
“sufficiently small” depends on unknowns involving $P$. One way to satisfy
$\omega<\bar{\omega}$ is with a positive, vanishing sequence
$\omega=\omega_{n}$. For any such $\omega_{n}$, even a sequence that vanishes
arbitrarily slowly, we have $\omega_{n}\in(0,\bar{\omega})$ for all
sufficiently large $n$. And, since $\omega_{n}$ vanishes slowly it has no
effect on the Gibbs posterior concentration rate; and see Sections 5.2 and
5.4. All we require to accommodate a vanishing learning rate is a slightly
stronger $\omega_{n}$-dependent version of the prior concentration bound in
(15).

######  Theorem 3.

Let $\varepsilon_{n}$ be a vanishing sequence and $\omega_{n}$ be a learning
rate sequence satisfying $n\omega_{n}\varepsilon_{n}^{r}\to\infty$ for a
constant $r>0$. Consider a Gibbs posterior distribution
$\Pi_{n}=\Pi_{n}^{\omega_{n}}$ in (5) based on this sequence of learning
rates. If the prior satisfies

$\log\Pi(\\{\theta:m(\theta,\theta^{\star})\vee
v(\theta,\theta^{\star})\leq\varepsilon_{n}^{r}\\})\geq-
Cn\omega_{n}\varepsilon_{n}^{r},$ (17)

and if Condition 1 holds for $\omega_{n}$, then the Gibbs posterior
distribution in (5), with learning rate sequence $\omega_{n}$, has
concentration rate $\varepsilon_{n}$ at $\theta^{\star}$ for a sufficiently
large constant $M>0$ in Definition 4.

The proof of Theorem 3 is almost identical to that of Theorem 1, hence
omitted. But to see the basic idea, we mention two key observations. First,
since Condition 1 holds for $\omega_{n}$ for all sufficiently large $n$, and
since $\omega_{n}$ is deterministic, by the same argument producing the bound
in (10), we get

$P^{n}N_{n}(A_{n})\leq
e^{-Kn\omega_{n}M_{n}^{r}\varepsilon_{n}^{r}}\quad\text{for all {\em
sufficiently large} $n$}.$

Second, the same argument producing the bound in (14) shows that

$P^{n}\bigl{(}D_{n}\leq\tfrac{1}{2}\Pi(G_{n})e^{-2n\omega_{n}\varepsilon_{n}^{r}}\bigr{)}\leq
P\left(D_{n}\leq\tfrac{1}{2}e^{-(C+2)n\omega_{n}\varepsilon_{n}^{r}}\right)\leq(n\omega_{n}\varepsilon_{n}^{r})^{-1}.$

Then the only difference between this situation and that in Theorem 1 is that,
here, the numerator bound only holds for “sufficiently large $n$,” instead of
for all $n$. This makes Theorem 3 is slightly weaker than Theorem 1 since
there are no finite-sample bounds.

When $\omega_{n}\equiv\omega$ the constant learning rate is absorbed by $C$
and there is no difference between the prior bounds in (15) and (17). But, the
prior probability assigned to the $(m\vee v)$-neighborhood of $\theta^{\star}$
does not depend on $\omega_{n}$, so if it satisfies (15), then the only way it
could also satisfy (17) is if $\varepsilon_{n}$ is bigger than it would have
been without a vanishing learning rate. In other words, Theorem 1 requires
$n\varepsilon_{n}^{r}\rightarrow\infty$ whereas Theorem 3 requires
$n\omega_{n}\varepsilon_{n}^{r}\rightarrow\infty$, which implies that for a
given vanishing $\omega_{n}$ Theorem 1 potentially can accommodate a faster
rate $\varepsilon_{n}$. Therefore, we see that a vanishing learning rate can
slow down the Gibbs posterior concentration rate if it does not vanish
arbitrarily slowly. There are applications that require the learning rate to
vanish at a particular $n-$dependent rate, and these tend to be those where
adjustments like in Section 4 are needed; see Sections 5.3.2 and 5.5.2.

Theorems 1–3 assume a deterministic learning rate, but there are instances in
which the learning rate may be data-dependent; see Section 2.2. One motivation
for a random learning rate comes from the potential for a vanishing learning
rate to slow concentration. As an alternative to a vanishing learning rate we
might consider using the data to determine a $\hat{\omega}_{n}$ that will
satisfy $\hat{\omega}_{n}<\bar{\omega}$ at least with high probability. For an
example involving a data-dependent learning rate see Section 5.2.

######  Theorem 4.

Fix a positive deterministic learning rate sequence $\omega_{n}$ such that the
conditions of Theorem 3 hold and as a result $\Pi_{n}^{\omega_{n}}$ has
asymptotic concentration rate $\varepsilon_{n}$. Consider a random learning
rate sequence $\hat{\omega}_{n}$ satisfying

$P^{n}(\omega_{n}/2<\hat{\omega}_{n}<\omega_{n})\rightarrow 1,\quad
n\to\infty.$ (18)

Then $\Pi_{n}^{\hat{\omega}_{n}}$, the Gibbs posterior distribution in (5)
scaled by the random learning rate sequence $\hat{\omega}_{n}$, also has
concentration rate $\varepsilon_{n}$ at $\theta^{\star}$ for a sufficiently
large constant $M>0$ in Definition 4.

###  3.4 Checking Condition 1

Of course, Condition 1 is only useful if it can be checked in practically
relevant examples. For the examples we consider in Section 5 in which
Condition 1 can be used, the loss differences
$\ell_{\theta}-\ell_{\theta^{\star}}$ are either bounded or sub-
exponential/sub-Gaussian.

####  3.4.1 Bounded excess losses

For bounded loss differences, i.e., where
$\ell_{\theta}(u)-\ell_{\theta^{\star}}(u)<C$ for all $(\theta,u)$, Lemma 7.26
in Lafftery et al., (2010) gives the following bound:

$Pe^{-\omega(\ell_{\theta}-\ell_{\theta^{\star}})}\leq\exp\bigl{[}-\omega
m(\theta,\theta^{\star})+\omega^{2}v(\theta,\theta^{\star})\bigl{\\{}\tfrac{\exp(C\omega)-1-C\omega}{C^{2}\omega^{2}}\bigr{\\}}\bigr{]}.$
(19)

Whether Condition 1 holds depends on the choice of $\omega$ and the
relationship between $m(\theta,\theta^{\star})$ and $v(\theta,\theta^{\star})$
as defined by the Bernstein condition:

$v(\theta,\theta^{\star})\leq
C_{1}m(\theta,\theta^{\star})^{\alpha},\quad\alpha\in(0,1],\quad C_{1}>0.$
(20)

Denote the bracketed expression in the exponent of (19) by
$C(\omega)=\tfrac{\exp(C\omega)-1-C\omega}{C^{2}\omega^{2}}$. When $\alpha=1$,
the $m$ and $v$ functions are of the same order and, if
$d(\theta,\theta^{\star})=m(\theta,\theta^{\star})$, and $\omega
C(\omega)\leq(2C_{1})^{-1}$ then Condition 1 holds with $K=1/2$. Since
$C(\omega)\approx C\omega$ for small $\omega$, it suffices to take the
learning rate a sufficiently small constant.

On the other hand, when $v$ is larger than $m$, i.e., the Bernstein condition
holds with $\alpha$ strictly less than 1, the only way to satisfy Condition 1
is to let $\omega$ depend on $\alpha$ and $m(\theta,\theta^{\star})$. Suppose
$m(\theta,\theta^{\star})\geq d(\theta,\theta^{\star})>\varepsilon_{n}$. For
all small enough $\varepsilon_{n}$, we can set
$\omega_{n}=(4C)^{-1/2}\varepsilon_{n}^{(1-\alpha)/2}$ and Condition 1 is
again satisfied with $K=1/2$.

The above strategies are implemented in the classification examples in Section
5.5 where we also discuss connections to the Tsybakov margin condition
(Tsybakov,, 2004) and the Massart noise condition (Massart and Nedelec,,
2006). For the former example, see Grünwald and Mehta, (2020), Theorem 22, and
the discussion that follows, where the learning rate is defined to be a
vanishing sequence depending on the concentration rate and the Bernstein
exponent $\alpha$.

####  3.4.2 Sub-exponential excess losses

Unbounded but light-tailed loss differences
$\ell_{\theta}-\ell_{\theta^{\star}}$ may also satisfy Condition 1. Both sub-
exponential and sub-Gaussian random variables (Boucheron et al.,, 2012, Sec.
2.3-4) admit an upper bound on their moment-generating functions. When the
loss difference $\ell_{\theta}-\ell_{\theta^{\star}}$ is sub-Gaussian

$Pe^{-\omega(\ell_{\theta}-\ell_{\theta^{\star}})}\leq\exp\bigl{\\{}\tfrac{\omega^{2}}{2}\sigma^{2}(\theta,\theta^{\star})-\omega
m(\theta,\theta^{\star})\bigr{\\}},$ (21)

for all $\omega$ and where the so-called variance proxy
$\sigma^{2}(\theta,\theta^{\star})$ may depend on $\theta$ and
$\theta^{\star}$. If $\ell_{\theta}-\ell_{\theta^{\star}}$ is sub-exponential,
then the above bound holds for all $\omega\leq(2b)^{-1}$ for some parameter
$b<\infty$ indexing the tail behavior of $P$.

If $\sigma^{2}(\theta,\theta^{\star})$ is upper-bounded by
$Lm(\theta,\theta^{\star})$ for a constant $L>0$, then the above bound can be
rewritten as

$Pe^{-\omega(\ell_{\theta}-\ell_{\theta^{\star}})}\leq\exp\bigl{\\{}-\omega
m(\theta,\theta^{\star})(1-\tfrac{\omega L}{2})\bigr{\\}}.$ (22)

Then Condition 1 holds for $\omega<2L^{-1}$ if the loss difference is sub-
Gaussian, and for $\omega<\min\\{2L^{-1},(2b)^{-1}\\}$ if the loss difference
is sub-exponential.

In practice it may be awkward to assume
$\ell_{\theta}(U)-\ell_{\theta^{\star}}(U)$ is sub-exponential, but in certain
problems it is sufficient to make such an assumption about features of $U$,
which may be more reasonable. See Section 5.4 for an application of this idea
to a fixed-design regression problem where the fact the response variable is
sub-Gaussian implies the excess loss is itself sub-Gaussian.

##  4 Extensions

###  4.1 Locally sub-exponential type loss functions

In some cases the moment generating function bound in Condition 1 can be
verified in a neighborhood of $\theta^{\star}$ but not for all
$\theta\in\Theta$. For example, consider a loss function such that
$\theta\mapsto\ell_{\theta}(u)$ is Lipschitz with respect to a metric
$\|\cdot\|$—not necessarily the divergence $d$—with uniformly bounded
Lipschitz constant $L=L(u)$, such that for some $\delta>0$

$\displaystyle\|\theta-\theta^{\star}\|<\delta$ $\displaystyle\implies
v(\theta,\theta^{\star})\lesssim m(\theta,\theta^{\star})$
$\displaystyle\|\theta-\theta^{\star}\|>\delta$ $\displaystyle\implies
v(\theta,\theta^{\star})\lesssim m(\theta,\theta^{\star})^{2}.$

In other words, the Bernstein condition is satisfied for different values of
$\alpha$ depending on $\theta$. This is the case in quantile estimation and
quantile regression problems; see Sections 5.1 and 5.6 below. The
characterization of these cases we give below is apparently new in the
literature. For example, Grünwald and Mehta, (2020) do not consider problems
in which the exponent in the Bernstein condition depends on $\theta$, and, as
a consequence, they can only verify their concentration rates for Lipschitz
loss functions when the parameter space is bounded, where the Bernstein
condition (20) holds trivially with $\alpha=0$; see Grünwald and Mehta,
(2020), Example 10.

For each $n$, partition $\Theta$ as
$\Theta=\Theta_{n}\cup\Theta_{n}^{\text{\sc c}}$, where
$\Theta_{n}=\\{\theta:\|\theta\|\leq\Delta_{n}\\}$. As the notation suggests,
$\Theta_{n}$ may depend on $n$; e.g., $\Delta_{n}$ may be an increasing
sequence so that $\Theta_{n}$ is a sieve, an increasing sequence of subsets of
$\Theta$. The important point is that, for sufficiently large $n$,
$\theta\in\Theta_{n}$ implies that $\|\theta-\theta^{\star}\|$ is bounded by a
multiple of $\Delta_{n}$. This, together with the Lipschitz property of the
loss and (19), implies

$Pe^{-\omega_{n}(\ell_{\theta}-\ell_{\theta^{\star}})}\leq\exp\bigl{\\{}C_{n}\omega_{n}^{2}v(\theta,\theta^{\star})-\omega_{n}m(\theta,\theta^{\star})\bigr{\\}},$

where $C_{n}=O(1+\Delta_{n}\omega_{n})$. Suppose that
$\omega_{n}\Delta_{n}=O(1)$, so we can replace $C_{n}$ by a constant $C$; this
is not a strong supposition since both $\omega_{n}$ and $\Delta_{n}$ are
chosen at the data analyst’s discretion. If there exist constants $k,K_{n}>0$
such that

$m(\theta,\theta^{\star})\geq
k\|\theta-\theta^{\star}\|^{2}\quad\text{and}\quad
v(\theta,\theta^{\star})\leq K_{n}\|\theta-\theta^{\star}\|^{2}\quad\text{for
all $\theta\in\Theta_{n}$},$ (23)

where $K_{n}$ may depend on $\Delta_{n}$, then the above upper bound
simplifies to

$\exp\\{-\omega_{n}\|\theta-\theta^{\star}\|^{2}(k-CK_{n}\omega_{n})\\}.$

If $\omega_{n}<\min\\{k(2CK_{n})^{-1},\Delta_{n}^{-1}\\}$, then Condition 1
holds for $\theta\in\Theta_{n}$. Since we generally need $\Delta_{n}$ to be an
increasing sequence, this requires we take a vanishing learning rate like in
Theorem 3. Condition 2 below describes this situation in greater generality.

######  Condition 2.

There exists an interval $(0,\bar{\omega})$, a constant $r>0$, and a sequence
$K_{n}>0$ such that for all $\omega\in(0,\bar{\omega})$ and for all
sufficiently small $\delta>0$,

$\theta\in\Theta_{n}\;\text{ and }\;d(\theta;\theta^{\star})>\delta\implies
Pe^{-\omega(\ell_{\theta}-\ell_{\theta^{\star}})}<e^{-K_{n}^{r}\omega\delta^{r}}$
(24)

Besides the requirement that $\theta$ belongs to the subset $\Theta_{n}$ of
the full parameter space, Conditions 1 and 2 differ in the exponential bound.
Where Condition 1 has a constant $K$ in the bound, the dependence on $n$ in
Condition 2 causes this value to depend on the sample size; hence, the
sequence $K_{n}$ appears in Condition 2 place of the constant $K$.

Condition 2 can be used exactly as Condition 1 to upper bound
$P^{n}N_{n}(A_{n}\cap\Theta_{n})$ by $e^{-\omega
nK_{n}^{r}M_{n}^{r}\varepsilon_{n}^{r}}$. Separately, we need to control the
Gibbs posterior probability of the set $A_{n}\cap\Theta_{n}^{c}$ in order to
prove concentration; see Theorem 5.

######  Theorem 5.

Let $\Theta_{n}$ be a sequence of subsets for which the loss function
satisfies Condition 2 for a sequence $K_{n}>0$, a constant $r>0$, and a
learning rate $\omega_{n}\in(0,\overline{\omega})$ for all sufficiently large
$n$. Let $\varepsilon_{n}$ be a vanishing sequence satisfying
$n\omega_{n}K_{n}^{r}\varepsilon_{n}^{r}\rightarrow\infty$. Suppose the prior
satisfies

$\log\Pi(\\{\theta:m(\theta,\theta^{\star})\vee
v(\theta,\theta^{\star})\leq(K_{n}\varepsilon_{n})^{r}\\})\gtrsim-
Cn\omega_{n}K_{n}^{r}\varepsilon_{n}^{r},$ (25)

for some $C>0$ and the same $K_{n}$, $r$, and $\omega_{n}$ as above. Then the
Gibbs posterior distribution in (5) satisfies

$\limsup_{n\to\infty}P^{n}\Pi_{n}(A_{n})\leq\limsup_{n\to\infty}P^{n}\Pi_{n}(A_{n}\cap\Theta_{n}^{\text{\sc
c}}).$

Consequently, if

$P^{n}\Pi_{n}(\Theta_{n}^{\text{\sc c}})\to 0\quad\text{as $n\to\infty$},$
(26)

then the Gibbs posterior has asymptotic concentration rate $\varepsilon_{n}$
for all large enough constants $M>0$ in Definition 4.

The condition in (26) may not be easy to check, but this challenge is not
unique to Gibbs posteriors—there are similar challenges faced in the analysis
of misspecified Bayesian posteriors. Of course, it can always be arranged by
restricting the support of the prior distribution to $\Theta_{n}$, which is
what Kleijn and van der Vaart, (2006) suggest in their Theorem 2.3 in the
Bayesian case and how we handle (26) for our infinite-dimensional Gibbs
posterior example in Section 5.6 and Appendix C.8. However, as (26) suggests,
this is not entirely necessary. Indeed, Kleijn and van der Vaart, (2006)
describe a sort of trade-off between model complexity and prior support, and
offer a more complicated form (Equation 2.13) of their sufficient condition,
which they do not explore in that paper. The fact is, without a well-specified
likelihood, checking a condition like our (26) or Equation (2.13) in Kleijn
and van der Vaart, (2006) is a challenge, at least in infinite-dimensional
problems. For finite-dimensional models, it may be possible to verify (26)
directly using properties of the loss functions. See Appendix B for an
intermediate result on convex loss functions which we apply to the example in
Section 5.1.

###  4.2 Surrogate sub-exponential type losses

When loss differences are unbounded and heavy-tailed, i.e., neither sub-
Gaussian nor sub-exponential, the moment-generating function in Condition 1
does not exist and the condition cannot be satisfied. A reasonable strategy in
such problems is to replace $\ell_{\theta}$ by the bounded surrogate loss
$\ell^{n}_{\theta}(u)=\ell_{\theta}(u)\wedge C_{n}$, where $C_{n}$ is a
diverging sequence. Following the strategy laid out in Section 3.4 to verify
Condition 1, we have

$Pe^{-\omega_{n}(\ell_{\theta}^{n}-\ell_{\theta^{\star}}^{n})}\leq\exp\bigl{[}-\omega_{n}m_{n}(\theta,\theta^{\star})+\omega_{n}^{2}v_{n}(\theta,\theta^{\star})\bigl{\\{}\tfrac{\exp(L_{n}\omega_{n})-1-L_{n}\omega_{n}}{L_{n}^{2}\omega_{n}^{2}}\bigr{\\}}\bigr{]},$

where
$m_{n}(\theta,\theta^{\star}):=P(\ell_{\theta}^{n}-\ell_{\theta^{\star}}^{n})$
and
$v_{n}(\theta,\theta^{\star}):=P(\ell_{\theta}^{n}-\ell_{\theta^{\star}}^{n})^{2}-m_{n}(\theta,\theta^{\star})^{2}$
are the mean and variance of the surrogate excess loss, and where we suppose
$\sup_{u}|\ell_{\theta}^{n}(u)-\ell_{\theta_{n}^{\star}}^{n}(u)|<L_{n}$, for
some diverging sequence $L_{n}$ depending on $C_{n}$. Suppose we take
$\omega_{n}=L_{n}^{-1}$ so that the bracketed expression in the exponent is
bounded by a constant. And suppose that, for a given divergence measure
$d(\theta,\theta^{\star})$, and for $C_{n}$ diverging sufficiently fast, we
have

$d(\theta,\theta^{\star})>\varepsilon_{n}\implies
m_{n}(\theta,\theta^{\star})\geq
c_{1}d(\theta,\theta^{\star})^{r}\quad\text{and}\quad
v_{n}(\theta,\theta^{\star})\leq c_{2}d(\theta,\theta^{\star})^{r}.$

This is not a strong supposition and is implied if
$m(\theta,\theta^{\star})=O(d(\theta,\theta^{\star})^{r})$, the Bernstein
condition holds with $\alpha=1$, and if $m_{n}(\theta,\theta^{\star})\approx
m(\theta,\theta^{\star})$ and $v_{n}(\theta,\theta^{\star})\approx
v(\theta,\theta^{\star})$ for large $n$. Then, the fact $\omega_{n}\rightarrow
0$ implies that

$P^{-\omega_{n}(\ell_{\theta}^{n}-\ell_{\theta^{\star}}^{n})}\leq
e^{-c_{1}L_{n}^{-1}\varepsilon_{n}^{r}/2},\quad\text{for all large $n$},$

and the following extended version of Condition 1 for surrogate, capped loss
functions holds with $K=1/2$.

######  Condition 3.

Let $\ell_{\theta}^{n}(u)=\ell_{\theta}(u)\wedge C_{n}$, for a sequence of
constants $C_{n}\to\infty$. There exists a sequence $\bar{\omega}_{n}>0$, and
constants $K,\,r>0$, such that for all sequences
$0<\omega_{n}\leq\bar{\omega}_{n}$ and for all sufficiently small $\delta>0$,
for $\theta\in\Theta$

$d(\theta;\theta^{\star})>\delta\implies
Pe^{-\omega_{n}(\ell_{\theta}^{n}-\ell^{n}_{\theta^{\star}})}<e^{-K\omega_{n}\delta^{r}}.$
(27)

######  Theorem 6.

Let $\varepsilon_{n}$ be a vanishing sequence, $C_{n}>0$ be a diverging
sequence, and $\omega_{n}$ be a learning rate sequence satisfying
$n\omega_{n}\varepsilon_{n}^{r}\to\infty$ for a constant $r>0$. Consider a
Gibbs posterior distribution $\Pi_{n}=\Pi_{n}^{\omega_{n}}$ in (5) based on
this sequence of learning rates. If the prior satisfies

$\log\Pi(\\{\theta:m_{n}(\theta,\theta^{\star})\vee
v_{n}(\theta,\theta^{\star})\leq\varepsilon_{n}^{r}\\})\geq-
C_{1}n\omega_{n}\varepsilon_{n}^{r},$ (28)

and if Condition 3 holds for the sequences $C_{n}$ and $\omega_{n}$, then the
Gibbs posterior distribution in (5), with learning rate sequence $\omega_{n}$,
has concentration rate $\varepsilon_{n}$ at $\theta^{\star}$ for a
sufficiently large constant $M>0$ in Definition 4.

Condition 3 and (28) are simply Condition 1 and (17) applied to the surrogate
loss function $\ell_{\theta}^{n}$ instead of the loss function $\ell_{\theta}$
that has risk minimized by $\theta^{\star}$. The proof of Theorem 6 can be
obtained simply by substituting $\ell_{\theta}^{n}$ for $\ell_{\theta}$ in the
proof of Theorem 3 and applying the corresponding conditions. Therefore, we
omit the proof of Theorem 6.

In Section 5.3, we apply Theorem 6 in the context of regression under squared-
error loss when the distribution of the response is potentially heavy-tailed,
i.e., with polynomial instead of exponential tails where the basic Condition 1
would typically hold.

##  5 Examples

This section presents several illustrations of the general theory presented in
Sections 3 and 4. The examples include uses of vanishing and even random
learning rate sequences. And, the strategies laid out in Sections 3.4, 4.1,
and 4.2 are put to use here to verify our sufficient conditions for Gibbs
posterior concentration. All proofs of results in this section can be found in
Appendix C.

###  5.1 Quantile regression

Consider making inferences on the $\tau^{\text{th}}$ conditional quantile of a
response $Y$ given a predictor $X=x$. We model this quantile, denoted
$Q_{Y|X=x}(\tau)$, as a linear combination of functions of $x$, that is,
$Q_{Y|X=x}(\tau)=\theta^{\top}f(x)$, for a fixed, finite dictionary of
functions $f(x)=(f_{1}(x),\dots,f_{J}(x))^{\top}$ and where
$\theta=(\theta_{1},\dots,\theta_{J})^{\top}$ is a coefficient vector with
$\theta\in\Theta$. Here we assume the model is well-specified so the true
conditional quantile is ${\theta^{\star}}^{\top}f(x)$ for some
$\theta^{\star}\in\Theta$. The standard _check loss_ for quantile estimation
is

$\ell_{\theta}(u)=(y-\theta^{\top}f(x))(\tau-1\\{y<\theta^{\top}f(x)\\}).$
(29)

See Appendix C.1 for a proof that $\theta^{\star}$ minimizes $R(\theta)$. It
can be shown that $\theta\mapsto\ell_{\theta}(u)$ is $L$-Lipschitz, with
$L<1$, and convex, so the strategy in Section 4.1 is helpful here for
verifying Condition 2 and Lemma 2 in Appendix B can be used to verify (26).

Inference on quantiles is a challenging problem from a Bayesian perspective
because the quantile is well-defined irrespective of any particular
likelihood. Sriram et al., (2013) interprets the check loss as the negative
log-density of an asymmetric Laplace distribution and constructs a
corresponding pseudo-posterior using this likelihood, but their posterior is
effectively a Gibbs posterior as in Definition 3.

With a few mild assumptions about the underlying distribution $P$, our general
result in Theorem 2 can be used to establish Gibbs posterior concentration at
rate $n^{-1/2}$.

######  Assumption 1.

  1. 1.

The marginal distribution of $X$ is such that $Pff^{\top}$ exists and is
positive definite;

  2. 2.

the conditional distribution of $Y$, given $X=x$, has at least one finite
moment and admits a continuous density $p_{x}$ such that
$p_{x}(\theta^{\star\top}f)$ is bounded away from zero for $P$-almost all $x$;
and

  3. 3.

the prior $\Pi$ has a density bounded away from 0 in a neighborhood of
$\theta^{\star}$.

######  Proposition 1.

Under Assumption 1, if the learning rate is sufficiently small, then the Gibbs
posterior concentrates at $\theta^{\star}$ with rate
$\varepsilon_{n}=n^{-1/2}$ with respect to
$d(\theta,\theta^{\star})=\|\theta-\theta^{\star}\|$.

###  5.2 Area under receiver operator characteristic curve

The receiver operator characteristic (ROC) curve and corresponding area under
the curve (AUC) are diagnostic tools often used to judge the effectiveness of
a binary classifier. Suppose a binary classifier produces a score $U$
characterizing the likelihood an individual belongs to Group $1$ versus Group
$0$. We can estimate an individual’s group by $1(U>t)$ where different values
of the cutoff score $t$ may provide more or less accurate estimates. Suppose
$U_{0}$ and $U_{1}$ are independent scores corresponding to individuals from
Group $0$ and Group $1$, respectively. The specificity and sensitivity of the
test of $H_{0}:\text{individual }i\text{ belongs to Group 0}$ that rejects
when $U>t$ are defined by $\text{spec}(t)=P(U_{0}<t)$ and
$\text{sens}(t)=P(U_{1}>t)$. When the type 1 and 2 errors of the test are
equally costly the optimal cutoff is the value of $t$ maximizing
$1-\text{spec}(t)+\text{sens}(t)$, or, in other words, the test maximizing the
sum of power and one minus the type 1 error probability. The ROC is the
parametric curve $(1-\text{spec}(t),\text{sens}(t))$ in $[0,1]^{2}$ which
provides a graphical summary of the tradeoff between Type 1 and Type 2 errors
for different choices of the cutoff. The AUC, equal to $P(U_{1}>U_{0})$, gives
an overall numerical summary of the quality of the binary classifier,
independent of the choice of threshold.

Our goal is to make posterior inferences on the AUC, but the usual Bayesian
approach immediately runs into the kinds of problems we see in the examples in
Sections 5.1 and later in 6. The parameter of interest is one-dimensional, but
it depends on a completely unknown joint distribution $P$. Within a Bayesian
framework, the options are to fix a parametric model for this joint
distribution and risk model misspecification or work with a complicated
nonparametric model. Wang and Martin, (2020) constructed a Gibbs posterior for
the AUC that avoids both of these issues.

Suppose $U_{0,1},\ldots,U_{0,m}$ and $U_{1,1},\ldots,U_{1,n}$ denote random
samples of size $m$ and $n$, respectively, of binary classifier scores for
individuals belonging to Groups 0 and 1, and denote $\theta=P(U_{1}>U_{0})$.
Wang and Martin, (2020) consider the loss function

$\ell_{\theta}(u_{0},u_{1})=\\{\theta-1(u_{1}>u_{0})\\}^{2},\quad\theta\in[0,1],$
(30)

for which the risk satisfies $R(\theta)=(\theta-\theta^{\star})^{2}$. If we
interpret $m=m_{n}$ as a function of $n$, then it makes sense to write the
empirical risk function as

$R_{n}(\theta)=\frac{1}{mn}\sum_{i=1}^{m}\sum_{j=1}^{n}\\{\theta-1(U_{1,i}>U_{0,j})\\}^{2}.$
(31)

Note the minimizer of the empirical risk is equal to

$\hat{\theta}_{n}=\frac{1}{mn}\sum_{i=1}^{m}\sum_{j=1}^{n}1(U_{1,i}>U_{0,j}).$
(32)

Wang and Martin, (2020) prove concentration of the Gibbs posterior at rate
$n^{-1/2}$ under the following assumption.

######  Assumption 2.

  1. 1.

The sample sizes $(m,n)$ satisfy $m(m+n)^{-1}\rightarrow\lambda\in(0,1)$.

  2. 2.

The prior distribution has a density function $\pi$ that is bounded away from
zero in a neighborhood of $\theta^{\star}$.

Wang and Martin, (2020) note that their concentration result holds for fixed
learning rates and deterministic learning rates that vanish more slowly that
$\min(m,n)^{-1}$. As discussed in Syring and Martin, (2019) one motivation for
choosing a particular learning rate is to calibrate Gibbs posterior credible
intervals to attain a nominal coverage probability, at least approximately.
With this goal in mind, Wang and Martin, (2020) suggest the following random
learning rate. Define the covariances

$\displaystyle\tau_{10}$
$\displaystyle=\text{Cov}\\{1(U_{1,1}>U_{0,1}),\,1(U_{1,1}>U_{0,2})\\}$
$\displaystyle\tau_{01}$
$\displaystyle=\text{Cov}\\{1(U_{1,1}>U_{0,1}),\,1(U_{1,2}>U_{0,1})\\}.$ (33)

Wang and Martin, (2020) note the asymptotic covariance of $\hat{\theta}_{n}$
is given by

$\frac{1}{m+n}\Bigl{(}\frac{\tau_{10}}{\lambda}+\frac{\tau_{01}}{1-\lambda}\Bigr{)},$
(34)

and that the Gibbs posterior variance can be made to match this, at least
asymptotically, by using the random learning rate

$\hat{\omega}_{n}=\frac{m+n}{2mn}\Bigl{(}\frac{\hat{\tau}_{10}}{\lambda}+\frac{\hat{\tau}_{01}}{1-\lambda}\Bigr{)}^{-1},$
(35)

where $\hat{\tau}_{10}$ and $\hat{\tau}_{01}$ are the corresponding empirical
covariances:

$\displaystyle\hat{\tau}_{10}$
$\displaystyle=\frac{2}{mn(n-1)}\sum_{i=1}^{m}\sum_{j\neq
j^{\prime}}1(U_{1,i}>U_{0,j})1(U_{1,i}>U_{0,j^{\prime}})-\hat{\theta}_{n}^{2},$
$\displaystyle\hat{\tau}_{01}$
$\displaystyle=\frac{2}{nm(m-1)}\sum_{j=1}^{n}\sum_{i\neq
i^{\prime}}1(U_{1,i}>U_{0,j})1(U_{1,i^{\prime}}>U_{0,j})-\hat{\theta}_{n}^{2}.$
(36)

The hope is that the Gibbs posterior with the learning rate $\hat{\omega}_{n}$
has asymptotically calibrated credible intervals. It turns out that the
concentration result in Wang and Martin, (2020) along with Theorem 4 imply the
Gibbs posterior with a slightly adjusted version of learning rate
$\hat{\omega}_{n}$ also concentrates at rate $n^{-1/2}$. The adjustment to the
learning rate has the effect of slightly widening Gibbs posterior credible
intervals, so their asymptotic calibration is not adversely affected.

######  Proposition 2.

Suppose Assumption 2 holds and let $a_{n}$ denote any diverging sequence.
Then, the Gibbs posterior with learning rate $a_{n}\hat{\omega}_{n}$
concentrates at rate $n^{-1/2}$ with respect to
$d(\theta,\theta^{\star})=|\theta-\theta^{\star}|$.

###  5.3 Regression with bounded predictions

Consider predicting a response $y\in\mathbb{R}$ using a linear function
$x^{\top}\theta$ by minimizing the sum of squared-error losses
$\ell_{\theta}(u)=(y-x^{\top}\theta)^{2}$, with $u=(x,y)\in\mathbb{R}^{J+1}$,
over a compact parameter space $\theta\in\Theta\subseteq\mathbb{R}^{J}$.
Suppose the covariate-response variable pairs $(X_{i},Y_{i})$ are iid with $X$
taking values in a compact subset $\mathcal{X}\subset\mathbb{R}^{J}$.
Compactness of the parameter and covariate spaces implies the predictions are
bounded, i.e., $|\theta^{\top}X|<B$ for some constant $B>0$. A more flexible,
non-parametric regression with unbounded predictions is discussed in Section
5.4 below. For this example, we focus on how the tail behavior of the response
variable affects posterior concentration; see Assumptions 3a and 3b below.

####  5.3.1 Thin-tailed response

It is reasonable to use bounded functions to predict thin-tailed responses
since extreme observations are unlikely. If the response is sub-exponential
then so is the excess loss, and by the argument outlined in Section 3.4 we can
verify Condition 1 for
$d(\theta,\theta^{\star})=\|\theta-\theta^{\star}\|_{2}$ and with $r=2$. The
Gibbs posterior distribution concentrates at rate $n^{-1/2}$ as a consequence
of Theorem 2.

######  Assumption 3a.

  1. 1.

The response $Y$, given $x$, is sub-exponential with parameters
$(\sigma^{2},\,b)$ for all $x$;

  2. 2.

the parameter space and covariate space are both compact subsets of
$\mathbb{R}^{J}$;

  3. 3.

the marginal distribution of $X$ is such that $PXX^{\top}$ exists and is
positive definite with eigenvalues bounded away from $0$; and

  4. 4.

the prior $\Pi$ has a density bounded away from $0$ in a neighborhood of
$\theta^{\star}$.

######  Proposition 3a.

If Assumption 3a holds, and the learning rate $\omega$ is a sufficiently small
constant, then the Gibbs posterior concentrates at rate
$\varepsilon_{n}=n^{-1/2}$ with respect to
$d(\theta,\theta^{\star})=\|\theta-\theta^{\star}\|_{2}$.

####  5.3.2 Heavy-tailed response

As discussed in Section 3.4, Condition 1 can be expected to hold when the
response is light-tailed, but not when it is heavy-tailed. Although it may not
be reasonable to use bounded functions to predict heavy-tailed responses,
restricting to bounded predictions simplifies the analysis and allows us to
focus specifically on how to handle heavy-tailed responses. In this example we
find that in order to accommodate heavy-tailed responses, we must modify the
loss function and, in turn, accept a slower concentration rate.

In particular, we follow the strategy of Section 4.2 and define a surrogate
loss that caps the squared-error loss by a number $t_{n}$ increasing in $n$,
that is, $\ell_{\theta}^{n}(u)=\ell_{\theta}(u)\wedge t_{n}$, where
$\ell_{\theta}(u)=(y-x^{\top}\theta)^{2}$ as above. Two key observations drive
our concentration result in Proposition 3b below. The first is that the tail
behavior of the response, given below in Assumption 3b and indexed by the
constant $s^{\star}>0$, implies the surrogate loss satisfies

$\displaystyle
m_{n}(\theta,\theta^{\star})\gtrsim\|\theta-\theta^{\star}\|_{2}^{2}-ct_{n}^{-(s^{\star}-3)/2},\,\text{and}$
$\displaystyle
v_{n}(\theta,\theta^{\star})\lesssim\|\theta-\theta^{\star}\|_{2}^{2}.$

If we choose $t_{n}$ and $\varepsilon_{n}$ such that
$\varepsilon_{n}^{2}t_{n}^{(s^{\star}-3)/2}\rightarrow\infty$ then the bounds
on $m_{n}$ and $v_{n}$ are of the same order for all $\theta$ satisfying
$\|\theta-\theta^{\star}\|_{2}\geq\varepsilon_{n}$. The second key observation
is that the cap of $t_{n}$ on the loss function implies the excess loss is, at
most, of the order $t_{n}^{1/2}$, which means the learning rate should be set
equal to $t_{n}^{-1/2}$. Finally, if we set $\varepsilon_{n}$ to satisfy
$n\omega_{n}\varepsilon_{n}^{2}\rightarrow\infty$, then Condition 3 holds with
$r=2$.

######  Assumption 3b.

The marginal distribution of $Y$ admits a density $f$ with polynomial tails,
i.e., there exists a constant $s^{\star}>3$ such that, for all sufficiently
large constants $C>0$,

$|y|>C\implies f(y)<|y|^{-s^{\star}}.$

Generally, one would be interested in the tail behavior of the conditional
distribution of $Y$, given $X$. However, since here the $X$\- and
$\theta$-spaces are compact, it is enough to consider the tails of the
marginal distribution of $Y$. Condition 3b would hold, for example, if the
marginal/conditional distribution of $Y$ were Student-t with more than two
degrees of freedom.

######  Proposition 3b.

Let $t_{n}=n^{2/(s-2)}$, $\omega_{n}=t_{n}^{-1/2}$, and $\varepsilon_{n}=(\log
n)^{1/2}n^{-(s-3)/(2s-4)}$ for any $3<s\leq s^{\star}$. If Assumption 3a(2–4).
and Assumption 3b hold, then the Gibbs posterior with learning rate
$\omega_{n}$ concentrates at rate $\varepsilon_{n}$ with respect to
$d(\theta,\theta^{\star})=\|\theta-\theta^{\star}\|_{2}$.

The minimum loss cap of $n^{2/(s^{\star}-2)}$, defined in Proposition 3b,
increases in $n$ but decreases in $s^{\star}$. This may seem counterintuitive
because a large value of $s^{\star}$ means the response $Y$ behaves more like
a sub-exponential random variable, for which no loss cap is necessary. And, in
other words, an uncapped loss is equivalent to a capped loss where the cap is
infinite. So we might expect the loss cap to be higher for large $s^{\star}$,
not lower. Instead, it may be helpful to interpret the loss cap $t_{n}$ in
terms of what level of cap is tolerable, rather than what level of cap is
necessary for concentration. When $s^{\star}$ is large the tails are light so
that even a small cap on the loss is tolerable, as it is very unlikely to
observe a loss larger than the cap. On the other hand, when $s^{\star}$ is
small it is likely to observe a large loss, so a small cap would substantially
alter the risk function. Therefore, as $s^{\star}$ decreases and tail behavior
becomes heavier larger caps are required so that the minimizer of the
surrogate risk is close to $\theta^{\star}$.

The surrogate loss strategy has two advantages. First, it enables us to derive
concentration rates via Condition 1 in problems for which that condition does
not hold naturally. Second, it relieves the data analyst from needing to make
a strong assumption about $P$. Since Proposition 3b holds for every
$s^{\star}\geq s$ the data analyst may choose to use the surrogate loss rather
than the squared-error loss as the default loss for mean regression. Capping
the loss can be viewed as a robust strategy that makes a trade off, requiring
a weaker condition on the response variable in exchange for a slower rate of
concentration.

###  5.4 Mean regression curve

Section 3 presented the Gibbs posterior concentration rate results in the
context of iid data, which was only to keep the presentation as simple as
possible. It is possible to extend the Gibbs formulation and the rate results
to other scenarios, in particular, independent but not identically distributed
(inid) data. Here we present the relevant details in the context of a specific
example.

Let $Y_{1},\ldots,Y_{n}$ be independent, where the marginal distribution of
$Y_{i}$ depends on a fixed covariate $x_{i}\in[0,1]$ through the mean, i.e.,
the expected value of $Y_{i}$ is $\theta^{\star}(x_{i})$, $i=1,\ldots,n$. For
simplicity, set $x_{i}=i/n$, corresponding to an equally-spaced design. Then
the goal is estimation of the mean function
$\theta^{\star}:[0,1]\to\mathbb{R}$, which resides in a specified function
class $\Theta$ defined below.

A natural starting point is to define an empirical risk based on squared error
loss, i.e.,

$R_{n}(\theta)=\frac{1}{n}\sum_{i=1}^{n}\\{Y_{i}-\theta(x_{i})\\}^{2}.$ (37)

However, any function $\theta$ that passes through the observations would be
an empirical risk minimizer, so some additional structure is needed to make
the solution to the empirical risk minimization problem meaningful. Towards
this, as is customary in the literature, we parametrize the mean function as a
linear combination of a fixed set of basis functions,
$f(x)=(f_{1}(x),\ldots,f_{J}(x))^{\top}$. That is, we consider only functions
$\theta=\theta_{\beta}$, where

$\theta_{\beta}(x)=\beta^{\top}f(x),\quad\beta\in\mathbb{R}^{J}.$ (38)

Note that we do not assume that $\theta^{\star}$ is of the specified form;
more specifically, we do not assume existence of a vector $\beta^{\star}$ such
that $\theta^{\star}=\theta_{\beta^{\star}}$. The idea is that the structure
imposed via the basis functions will force certain smoothness, etc., so that
minimization of the risk over this restricted class of functions would
identify a suitable estimate.

This structure changes the focus of our investigation from the mean function
$\theta$ to the $J$-vector of coefficients $\beta$. We now proceed by first
constructing a Gibbs posterior for $\beta$ and then obtain the corresponding
Gibbs posterior for $\theta$ by pushing the former through the mapping
$\beta\mapsto\theta_{\beta}$. In particular, define the empirical risk
function in terms of $\beta$:

$r_{n}(\beta)=R_{n}(\theta_{\beta})=\frac{1}{n}\sum_{i=1}^{n}\\{Y_{i}-\theta_{\beta}(x_{i})\\}^{2}=\tfrac{1}{n}(Y-F_{n}\beta)^{\top}(Y-F_{n}\beta),$
(39)

where $\beta\in\mathbb{R}^{J}$ and where $F_{n}$ is the $n\times J$ matrix
whose $(i,j)$ entry is $f_{j}(x_{i})$, assumed to be positive definite; see
below. Given a prior distribution $\widetilde{\Pi}$ for $\beta$—which
determines a prior $\Pi$ for $\theta$ through the aforementioned mapping—we
can first construct the Gibbs posterior for $\beta$ as in (5) with the pseudo-
likelihood $\beta\mapsto\exp\\{-\omega nr_{n}(\beta)\\}$. If we write
$\widetilde{\Pi}_{n}$ for this Gibbs posterior for $\beta$, then the
corresponding Gibbs posterior for $\theta$ is given by

$\Pi_{n}(A)=\widetilde{\Pi}_{n}(\\{\beta:\theta_{\beta}\in A\\}),\quad
A\subseteq\Theta.$ (40)

Therefore, the concentration properties of $\Pi_{n}$ are determined by those
of $\widetilde{\Pi}_{n}$.

We can now proceed very much like we did before, but the details are slightly
more complicated in the present inid case. Taking expectation with respect to
the joint distribution of $(Y_{1},\ldots,Y_{n})$ is, as usual, the average of
marginal expectations; however, since the data are not iid, these marginal
expectations are not all the same. Therefore, the expected empirical risk
function is

$\displaystyle\bar{r}_{n}(\beta)=P^{n}r_{n}(\beta)$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}P_{i}\\{Y_{i}-\theta_{\beta}(x_{i})\\}^{2}$
(41)

where $P_{i}=P_{x_{i}}$ is the marginal distribution of $Y_{i}$ and where
$F_{n,i}$ is the $i^{\text{th}}$ row of $F_{n}$. Since the expected empirical
risk function depends on $n$, through $(x_{1},\ldots,x_{n})$, so too does the
risk minimizer

$\beta^{\dagger}_{n}=\arg\min_{\beta}\bar{r}_{n}(\beta).$ (42)

If $P_{i}$ has finite variance, then $\bar{r}_{n}(\beta)$ differs from
$\\{\theta^{\star}(x_{1:n})-F_{n}\beta\\}^{2}$ by only an additive constant
not depending on $\beta$, and this becomes a least-squares problem, with
solution

$\beta_{n}^{\dagger}=(F_{n}^{\top}F_{n})^{-1}F_{n}^{\top}\,\theta^{\star}(x_{1:n}),$
(43)

where $\theta^{\star}(x_{1:n})$ is the $n$-vector
$(\theta^{\star}(x_{1}),\ldots,\theta^{\star}(x_{n}))^{\top}$. Our expectation
is that the Gibbs posterior $\widetilde{\Pi}_{n}$ for $\beta$ will suitably
concentrate around $\beta_{n}^{\dagger}$, which implies that the Gibbs
posterior $\Pi_{n}$ for $\theta$ will suitably concentrate around
$\theta_{\beta_{n}^{\dagger}}$. Finally, if the above holds and the basis
representation is suitably flexible, then $\theta_{\beta_{n}^{\dagger}}$ will
be close to $\theta^{\star}$ in some sense and, hence, we achieve the desired
concentration.

The flexibility of the basis representation depends on the dimension $J$.
Since $\theta^{\star}$ need not be of the form $\theta_{\beta}$, a good
approximation will require that $J=J_{n}$ be increasing with $n$. How fast
$J=J_{n}$ must increase depends on the smoothness of $\theta^{\star}$. Indeed,
if $\theta^{\star}$ has smoothness index $\alpha>0$ (made precise below), then
many systems of basis functions—including Fourier series and B-splines—have
the following approximation property: there exists an $H>0$ such that for
every $J$

$\text{there exists $\beta\in\mathbb{R}^{J}$ such that $\|\beta\|_{\infty}<H$
and $\|\theta_{\beta}-\theta^{\star}\|_{\infty}\lesssim J^{-\alpha}$}.$ (44)

Then the idea is to set the approximation error in (44) equal to the target
rate of convergence, which depends on $n$ and on $\alpha$, and then solve for
$J=J_{n}$.

For Gibbs posterior concentration at or near the optimal rate, we need the
prior distribution for $\beta$ to be sufficiently concentrated in a bounded
region of the $J$-dimensional space in the sense that

$\widetilde{\Pi}(\\{\beta:\|\beta-\beta^{\prime}\|_{2}\leq\varepsilon\\})\gtrsim(C\varepsilon)^{J},\quad\text{for
all $\beta^{\prime}\in\mathbb{R}^{J}$ with $\|\beta^{\prime}\|_{\infty}\leq
H$},$ (45)

for the same $H$ as in (44), for a small constant $C>0$, and for all small
$\varepsilon>0$.

######  Assumption 4.

  1. 1.

The function $\theta^{\star}:[0,1]\to\mathbb{R}$ belongs to a class
$\Theta=\Theta(\alpha,L)$ of Hölder smooth functions parametrized by
$\alpha\geq 1/2$ and $L>0$. That is, $\theta^{\star}$ satisfies

$|\theta^{\star([\alpha])}(x)-\theta^{\star([\alpha])}(x^{\prime})|\leq
L|x-x^{\prime}|^{\alpha-[\alpha]},\quad\text{for all $x,x^{\prime}\in[0,1]$},$

where the superscript “$(k)$” means $k^{\text{th}}$ derivative and $[\alpha]$
is the integer part of $\alpha$;

  2. 2.

for a given $x$, the response $Y$ is sub-Gaussian and with variance and
variance proxy—both of which can depend on $x$—uniformly upper bounded by
$\sigma^{2}$;

  3. 3.

the eigenvalues of $F_{n}^{\top}F_{n}$ are bounded away from zero and
$\infty$;

  4. 4.

the approximation property (44) holds; and

  5. 5.

the prior distribution $\widetilde{\Pi}$ for $\beta$ satisfies (45) and has a
bounded density on the $J_{n}-$dimensional parameter space.

The bounded variance assumption is implied, for example, if the variance of
$Y$ is a smooth function of $x$ in $[0,1]$, which is rather mild. And assuming
the eigenvalues of $F_{n}^{\top}F_{n}$ are bounded is not especially strong
since, in many cases, the basis functions would be orthonormal. In that case,
the diagonal and off-diagonal entries of $F_{n}^{\top}F_{n}$ would be
approximately 1 and 0, respectively, and the bounds are almost trivial. The
conditions on the prior distribution are weak; as we argue in the proof, it
can be satisfied by taking the joint prior density to be the product of $J$
independent prior densities on the components of $\beta$, and where each
component density is strictly positive.

######  Proposition 4.

If Assumption 4 holds, and the learning rate $\omega$ is a sufficiently small
constant, then the Gibbs posterior $\Pi_{n}$ for $\theta$ concentrates at
$\theta^{\star}$ with rate $\varepsilon_{n}=n^{-\alpha/(1+2\alpha)}$ with
respect to the empirical $L_{2}$ norm $\|\theta-\theta^{\star}\|_{n,2}$, where

$\|f\|_{n,2}^{2}=\frac{1}{n}\sum_{i=1}^{n}f^{2}(x_{i}).$

We should emphasize that the quantity of interest, $\theta$, is high-
dimensional, and the rate $\varepsilon_{n}$ given in Proposition 4 is optimal
for the given smoothness $\alpha$. The proof of Proposition 4 uses a slight
refinement of Theorem 1—essentially the same argument used to prove Theorem 2.
This highlights both the strength and flexibility of that result.

We have chosen to present a high-dimensional regression problem in this
section in order to push Theorem 1 to its limits. But it is also possible to
apply our results to the fixed-dimensional case with constant $J$. In that
case, we can simultanesously weaken the requirement on the response $Y$ in
Assumption 4.2 from sub-Gaussian to sub-exponential, and strengthen the
conclusion to a root-$n$ concentration rate.

###  5.5 Binary classification

Let $Y\in\\{0,1\\}$ be a binary response variable and
$X=(X_{0},X_{1},\ldots,X_{q})^{\top}$ a $(q+1)$-dimensional predictor
variable. We consider classification rules of the form

$\phi_{\theta}(X)=1\\{X^{\top}\theta>0\\}=1\\{\alpha
X_{0}+(X_{1},\ldots,X_{q})^{\top}\beta>0\\},\quad\theta=(\alpha,\beta)\in\mathbb{R}^{q+1},$
(46)

and the goal is to learn the optimal $\theta$ vector, i.e.,
$\theta^{\star}=\arg\min_{\theta}R(\theta)$, where
$R(\theta)=P\\{Y\neq\phi_{\theta}(X)\\}$ is the misclassification error
probability, and $P$ is the joint distribution of $(X,Y)$. This optimal
$\theta^{\star}$ is such that $\eta(x)>\frac{1}{2}$ if
$x^{\top}\theta^{\star}>0$ and $\eta(x)<\frac{1}{2}$ if
$x^{\top}\theta^{\star}\leq 0$, where $\eta(x)=P(Y=1\mid X=x)$ is the
conditional probability function. Below we construct a Gibbs posterior
distribution that concentrates around this optimal $\theta^{\star}$ at rate
that depends on certain local features of that $\eta$ function.

Suppose our data consists of iid copies $(X_{i},Y_{i})$, $i=1,\ldots,n$, of
$(X,Y)$ from $P^{\star}$, and define the empirical risk function

$R_{n}(\theta)=\frac{1}{n}\sum_{i=1}^{n}1\\{Y_{i}\neq\phi_{\theta}(X_{i})\\}.$
(47)

In addition to the empirical risk function we need specify a prior $\Pi$ and
here the prior plays a significant role in the Gibbs posterior concentration
results.

A unique feature of this problem, which makes the prior specification a little
different than in a linear regression problem, is that the scale of $\theta$
does not affect classification performance, e.g., replacing $\theta$ with
$1000\theta$ gives exactly the same classification performance. To fix a
scale, we follow Jiang and Tanner, (2008), and

  * •

assume that the $x_{0}$ component of $x$ is of known importance and always
included in the classifier,

  * •

and constrain the corresponding coefficient, $\alpha$, to take values $\pm 1$.

This implies that the $\alpha$ and $\beta$ components of the $\theta$ vector
should be handled very differently in terms of prior specification. In
particular, $\alpha$ is a scalar with a discrete prior—which we take here to
be uniform on $\pm 1$—and $\beta$, being potentially high-dimensional, will
require setting-specific prior considerations.

The characteristic that determines the difficulty of a classification problem
is the distribution of $\eta(X)$ or, more specifically, how concentrated
$\eta(X)$ is near the value $\frac{1}{2}$, where one could do virtually no
worse by classifying according to a coin flip. The set
$\\{x:\eta(x)=\frac{1}{2}\\}$ is called the margin, and conditions that
control the concentration of the distribution of $\eta(X)$ around
$\frac{1}{2}$ are generally called margin conditions. Roughly, if $\eta$ has a
“jump” or discontinuity at the margin, then classification is easier and
$\eta(X)$ does not have to be so tightly concentrated around $\frac{1}{2}$. On
the other hand, if $\eta$ is smooth at the margin, then the classification
problem is more challenging in the sense that more data near the margin is
needed to learn the optimal classifier, hence, tighter concentration of
$\eta(X)$ near $\frac{1}{2}$ is required.

In the two subsubsections that follow, we consider two such margin conditions,
namely, the so-called Massart and Tsybakov conditions. The first is relatively
strong, corresponding to a jump in $\eta$ at the margin, and the result we we
establish in Proposition 5a is accordingly strong. In particular, we show that
the Gibbs posterior achieves the optimal and adaptive concentration rate in a
class of high-dimensional problems ($q\gg n$) under a certain sparsity
assumption on $\theta^{\star}$. The Tsybakov margin condition we consider
below is weaker than the first, in the sense that $\eta$ can be smooth near
the “$\eta=\frac{1}{2}$” boundary and, as expected, the Gibbs posterior
concentration rate result is not as strong as the first. See Section 5.5.2 for
those details.

####  5.5.1 Massart’s noise condition

Here we allow the dimension $q+1$ of the coefficient vector
$\theta=(\alpha,\beta)$ to exceed the sample size $n$, i.e., we consider the
so-called high-dimensional problem with $q\gg n$. Accurate estimation and
inference is not possible in high-dimensional settings without imposing some
low-dimensional structure on the inferential target, $\theta^{\star}$. Here,
as is typical in the literature on high-dimensional inference, we assume that
$\theta^{\star}$ is sparse in the sense that most of its entries are exactly
zero, which corresponds to most of the predictor variables being irrelevant to
classification. Below we construct a Gibbs posterior distribution for $\theta$
that concentrates around the unknown sparse $\theta^{\star}$ at a (near)
optimal rate.

Since the sparsity in $\theta^{\star}$ is crucial to the success of any
statistical method, the prior needs to be chosen carefully so that sparsity is
encouraged in the posterior. The prior $\Pi$ for $\theta$ will treat $\alpha$
and $\beta$ independent, and the prior for $\beta$ will be defined
hierarchically. Start with the reparametrization $\beta\to(S,\beta_{S})$,
where $S\subseteq\\{1,2,\ldots,q\\}$ denotes the configuration of zeros and
non-zeros in the $\beta$ vector, and $\beta_{S}$ denotes the $|S|$-vector of
non-zero values. Following Castillo et al., (2015), for the marginal prior
$\pi(S)$ for $S$, we take

$\pi(S)=\textstyle\binom{q}{|S|}^{-1}\,f(|S|),$

where the $f$ is a prior for the size $|S|$ and the first factor on the right-
hand side is the uniform prior for $S$ of the given size $|S|$. Various
choices of $f$ are possible, but here we take the complexity prior
$f(s)\propto(cq^{a})^{-s},\,s=0,1,\ldots,q$, a truncated geometric density,
where $a$ and $c$ are fixed (and here arbitrary) hyperparameters; a similar
choice is also made in Martin et al., (2017). Second, for the conditional
prior of $\beta_{S}$, given $S$, again following Castillo et al., (2015), we
take its density to be

$g_{S}(\beta_{S})=\prod_{k\in S}\tfrac{\lambda}{2}e^{-\lambda|\beta_{k}|},$

a product of $|S|$ many Laplace densities with rate $\lambda$ to be specified.

######  Assumption 5a.

  1. 1.

The marginal distribution of $X$ is compactly supported, say, on
$[-1,1]^{q+1}$.

  2. 2.

The conditional distribution of $X_{0}$, given
$\tilde{X}=(X_{1},\ldots,X_{q})$, has a density with respect to Lebesgue
measure that is uniformly bounded.

  3. 3.

The rate parameter $\lambda$ in the Laplace prior satisfies
$\lambda\lesssim(\log q)^{1/2}$.

  4. 4.

The optimal $\theta^{\star}=(\alpha^{\star},\beta^{\star})$ is sparse in the
sense that $|S^{\star}|\log q=o(n)$, where $S^{\star}$ is the configuration of
non-zero entries in $\beta^{\star}$, and $\|\beta^{\star}\|_{\infty}=O(1)$.

  5. 5.

There exists $h\in(0,\frac{1}{2})$ such that $P(|\eta(X)-\frac{1}{2}|\leq
h)=0$.

The first two parts of Assumption 5a correspond to Conditions $0^{\prime}$ and
$0^{\prime\prime}$ in
Jiang and Tanner, (2008). and Assumption 5a (5) above is precisely the margin
condition imposed in Equation (5) of Massart and Nedelec, (2006); see, also,
Mammen and Tsybakov, (1999) and Koltchinskii, (2006). This concisely states
that there is either a jump in the $\eta$ function at the margin or that the
marginal distribution of $X$ is not supported near the margin; in either case,
there is separation between the $Y=0$ and $Y=1$ cases, which makes the
classification problem relatively easy. With these assumptions, we get the
following Gibbs posterior asymptotic concentration rate result. Note that, in
order to preserve the high-dimensionality in the asymptotic limit, we let the
dimension $q=q_{n}$ increase with the sample size $n$. So, the data sequence
actually forms a triangular array but, as is common in the literature, we
suppress this formulation in our notation.

######  Proposition 5a.

Consider a classification problem as described above, with $q\gg n$. Under
Assumption 5a, the Gibbs posterior, with sufficiently small constant learning
rate, concentrates at rate $\varepsilon_{n}=(n^{-1}|S^{\star}|\log q)^{1/2}$
with respect to
$d(\theta,\theta^{\star})=\\{R(\theta)-R(\theta^{\star})\\}^{1/2}$.

This result shows that, even in very high dimensional settings, the Gibbs
posterior concentrates on the optimal rule $\theta^{\star}$ at a fast rate.
For example, suppose that the dimension $q$ is polynomial in $n$, i.e., $q\sim
n^{b}$ for any $b>0$, while the “effective dimension,” or complexity, is sub-
linear, i.e., $|S^{\star}|\sim n^{a}$ for $a<1$. Then we get that
$\\{\theta:R(\theta)-R(\theta^{\star})\lesssim n^{-(1-a)}\log n\\}$ has Gibbs
posterior probability converging to 1 as $n\to\infty$. That is, rates better
than $n^{-1/2}$ can easily be achieved, and even arbitrarily close to $n^{-1}$
is possible. Compare this to the rates in Propositions 2–3 in Jiang and
Tanner, (2008), also in terms of risk difference, that cannot be faster than
$n^{-1/2}$.

Further, the concentration rate in Proposition 5a is nearly the optimal rate
corresponding to an oracle who has knowledge of $S^{\star}$. That is, the
Gibbs posterior concentrates at nearly the optimal rate adaptively with
respect to the unknown complexity.

####  5.5.2 Tsybakov’s margin condition

Next, we consider classification under the more general Tsybakov margin
condition (e.g., Tsybakov,, 2004). The problem set up remains the same as in
the previous example, except that we here consider the simpler low-dimensional
case, with the number of predictors $(q+1)$ small relative to $n$. Below we
briefly mention why we do not consider the high-dimensional case under the
Tsybakov condition. Since the dimension is no longer large, prior
specification is much simpler. We will continue to assume, as before, that the
$x_{0}$ component of $x$ has a constrained coefficient $\alpha\in\\{\pm 1\\}$,
to which we assign a discrete uniform prior. Otherwise, we simply require the
prior $\Pi$ is such that it has a (well-behaved) marginal density, $\pi$, for
$\beta$, with respect to Lebesgue measure on $\mathbb{R}^{q}$.

######  Assumption 5b.

  1. 1.

The marginal prior density $\pi$ for $\beta$ is continuous and bounded away
from 0 in a neighborhood of $\beta^{\star}$.

  2. 2.

There exists $c>0$ and $\gamma>0$ such that $P(|2\eta(X)-1|\leq h)\leq
ch^{\gamma}$ for all sufficiently small $h>0$.

The concentration of the marginal distribution of $\eta(X)$ around
$\frac{1}{2}$ controls the difficulty of the classification problem, and
Condition 5b(2) concerns exactly this. Note that smaller $\gamma$ implies
$\eta(X)$ is less concentrated around $\frac{1}{2}$, so we expect our Gibbs
posterior concentration rate, say, $\varepsilon_{n}=\varepsilon_{n}(\gamma)$,
to be a decreasing function of $\gamma$. The following proposition confirms
this.

######  Proposition 5b.

Suppose Assumption 5b holds and, for the specified $\gamma>0$, let
$\varepsilon_{n}=(\log n)^{\gamma/(2+2\gamma)}n^{-\gamma/(3+2\gamma)}$. Then
the Gibbs posterior distribution, with learning rate
$\omega_{n}=\varepsilon_{n}^{1/\gamma}$, concentrates at rate
$\varepsilon_{n}$ with respect to
$d(\theta,\theta^{\star})=\\{R(\theta)-R(\theta^{\star})\\}^{1/2}$.

Note that Massart’s condition from Section 5.5.1 corresponds to Tsybakov’s
condition above with $\gamma=\infty$. In that case, the Gibbs posterior
concentration rate we recover from Proposition 5b is $\varepsilon_{n}=(\log
n)^{1/2}n^{-1/2}$, which is achieved with a suitable constant learning rate.
This is within a logarithmic factor of the optimal rate for finite-dimensional
problems. Moreover, for both the $\gamma<\infty$ and $\gamma=\infty$ cases, we
expect that the logarithmic factor could be removed following an approach like
that described in Theorem 2, but we do not explore this possible extension
here.

We should emphasize that this case is unusual because the learning rate
$\omega_{n}$ depends on the (likely unknown) smoothness exponent $\gamma$.
This means the rate in Proposition 5b is not adaptive to the margin. However,
this dependence is not surprising, as it also appears in Grünwald and Mehta,
(2020), Section 6. The reason the learning rate depends on $\gamma$ is that
the Tsybakov margin condition in Assumption 5b(2) implies the Bernstein
condition in (20) takes the form

$v(\theta,\theta^{\star})\lesssim
m(\theta,\theta^{\star})^{\gamma/(1+\gamma)}.$

Therefore, in order to verify Condition 1 using the strategy in Section 3.4,
we need $\omega_{n}m(\theta,\theta^{\star})$ and
$\omega_{n}^{3}v(\theta,\theta^{\star})$ to have the same order when
$d(\theta,\theta^{\star})\geq M_{n}\varepsilon_{n}$. This requires that the
learning rate depends on $\gamma$, in particular,
$\omega_{n}=\varepsilon_{n}^{1/\gamma}$.

Admittedly, it is undesirable for the concentration rate to depend on the
unknown smoothness $\gamma$, but at least in this case $\varepsilon_{n}$
depends on no other unknowns. For the sparse classification example in Section
5.5.1, concentration adapts to the unknown dimension $|S^{\star}|$, a
consequence of the special choice of prior distribution. A similar argument
could be made in that example to accommodate the Tsybakov margin condition,
but it would require the learning rate to depend on both $\gamma$ and, worse,
on $|S^{\star}|$, which is impractical.

###  5.6 Quantile regression curve

In this section we revisit inference on a conditional quantile, covered in
Section 5.1. The $\tau^{\text{th}}$ conditional quantile of a response $Y$
given a covariate $X=x$ is modeled by a linear combination of basis functions
$f(x)=(f_{1}(x),...,f_{J}(x))^{\top}$:

$Q_{Y|X=x}(\tau)=\beta^{\top}f(x),\quad\beta\in\mathbb{R}^{J}.$

In Section 5.1, we made the rather restrictive assumption that the true
conditional quantile function $\theta^{\star}(x)$ belonged to the span of a
fixed set of $J$ basis functions. In practice, it may not be possible to
identify such a set of functions, which is why we considered using a sample-
size dependent sequence of sets of basis functions in Section 5.4 to model a
smooth function, $\theta^{\star}$. When the degree of smoothness, $\alpha$, of
$\theta^{\star}$ is known we can choose the number of basis functions to use
in order to achieve the optimal concentration rate. But, in practice, $\alpha$
may not be known, which creates a challenge because, as mentioned, the number
of terms needed in the basis function expansion modeling $\theta^{\star}$
depends on this unknown degree of smoothness.

To achieve optimal concentration rates adaptive to unknown smoothness, the
choice of prior is crucial. In particular, the prior must support a very large
model space in order to guarantee it places sufficient mass near
$\theta^{\star}$. Our approximation of $\theta^{\star}$ by a linear
combination of basis functions suggests a hierarchical prior for
$\theta\equiv(J,\beta_{J})$, similar to Section 5.5.1, with a marginal prior
$\pi$ for the number of basis functions $J$ and a conditional prior
$\widetilde{\Pi}_{J}$ for $\beta_{J}$, given $J$. The resulting prior for
$\theta$ is given by a mixture,

$\Pi(A)=\sum_{j=1}^{\infty}\pi(j)\,\widetilde{\Pi}_{j}(\\{\beta_{j}\in\mathbb{R}^{j}:\beta_{j}^{\top}f\in
A\\}),\quad A\subseteq\Theta.$ (48)

Then, in order for $\Pi$ to place sufficient mass near $\theta^{\star}$, it is
sufficient the marginal and conditional priors satisfy the following
conditions: the marginal prior $\pi$ for $J$ satisfies for some $c_{1}>0$ for
every $J=j$

$\pi(j)\geq e^{-c_{1}j\log j};$ (49)

and, the conditional prior $\widetilde{\Pi}$ for $\beta_{J}$ given $J=j$
satisfies for every $j$

$\widetilde{\Pi}(\\{\beta:\|\beta-\beta^{\prime}\|_{2}\leq\varepsilon\\})\gtrsim
e^{-Cj\log(1/\varepsilon)},\quad\text{for all
$\beta^{\prime}\in\mathbb{R}^{j}$ with $\|\beta^{\prime}\|_{\infty}\leq H$},$
(50)

for the same $H$ as in (44) and for some constant $C>0$ for all sufficiently
small $\varepsilon>0$. Fortunately, many simple choices of
$(\pi,\widetilde{\Pi}_{J})$ are satisfactory for obtaining adaptive
concentration rates, e.g., a Poisson prior on $J$ and a $J$-dimensional normal
conditional prior for $\beta$, given $J$; see Conditions (A1) and (A2) and
Remark 1 in Shen and Ghosal, (2015). Besides the conditions in (49) and (50)
we need to make a minor modification of $\Pi$ to make it suitable for our
proof of Gibbs posterior concentration; see below.

Similar to our choice in Section 5.1 we link the data and parameter through
the check loss function

$\ell_{\theta}(u)=\tfrac{1}{2}(|\theta(x)-y|-|y|)+(1-\tau)\theta(x),$

where $\theta(x)=\beta^{\top}f(x)$. See Koltchinskii, (1997) for a proof that
$P\ell_{\theta}$ is minimized at $\theta^{\star}$. It is straightforward to
show the check loss $\theta\mapsto\ell_{\theta}(u)$ is $L$-Lipschitz with
$L<1$. From there, if $Y$ were bounded we could use Condition 1 to compute the
concentration rate. However, to handle an unbounded response we need the
flexibility of Condition 2 and Theorem 5. To verify (26), we found it
necessary to limit the complexity of the parameter space by imposing a
constraint on the prior distribution, namely that the sequence of prior
distributions places all its mass on the set
$\Theta_{n}:=\\{\theta:\|\theta\|_{\infty}\leq\Delta_{n}\\}$ for some
diverging sequence $\Delta_{n}$; see Assumption 6.4. This constraint implies
the prior depends on $n$, and we refer to this sequence of prior distributions
by $\Pi^{(n)}$. Given the hierarchical prior $\Pi$ in (48) one straightforward
way to define a sequence of prior distributions satisfying the constraint is
to restrict and renormalize $\Pi$ to $\Theta_{n}$, i.e., define $\Pi^{(n)}$ as

$\displaystyle\Pi^{(n)}(A)=\Pi(A)/\Pi(\Theta_{n}),\quad
A\subseteq\Theta\cap\Theta_{n}$ (51)

This particular construction of $\Pi^{(n)}$ in (51) is not the only way to
define a sequence of priors satisfying the restriction to $\Theta_{n}$, but it
is convenient. That is, if $\Pi$ places mass $\eta$ on a sup-norm neighborhood
of $\theta^{\star}$ (see the proof of Proposition 6), then, by construction,
$\Pi^{(n)}$ in (51) places at least mass $\eta$ on the same neighborhood.

We should emphasize this restriction of the prior to $\Theta_{n}$ is only a
technical requirement needed for the proof, but it is not unreasonable. Since
the true function $\theta^{\star}$ is bounded, it is eventually in the growing
support of the prior $\Pi$. Similar assumptions have been used in the
literature on quantile curve regression; for example, Theorem 6 in Takeuchi et
al., (2006) requires that the parameter space consists only of bounded
functions, which is a stricter assumption than ours here.

######  Assumption 6.

  1. 1.

The function $\theta^{\star}:\mathbb{X}\mapsto\mathbb{R}$ is Hölder smooth
with parameters $(\alpha,L)$ (see Assumption 4.1);

  2. 2.

the basis functions satisfy the approximation property in (44);

  3. 3.

the covariate space $\mathbb{X}$ is compact and there exists a $\delta>0$ such
that the conditional density of $Y$, given $X=x$, is continuous and bounded
away from $0$ by a constant $\beta>0$ in the interval
$(\theta^{\star}(x)-\delta,\,\theta^{\star}(x)+\delta)$ for every $x$; and,

  4. 4.

the sequence $\Pi^{(n)}$ of prior distributions satisfies (51) for a sequence
of subsets of the parameter space
$\Theta_{n}:=\\{\theta:\|\theta\|_{\infty}\leq\Delta_{n}\\}$ for some sequence
$\Delta_{n}>0$, for $\Pi$ as defined in (48), and for marginal and conditional
priors $(\pi,\widetilde{\Pi})$ for $J$ and $\beta_{J}$ given $J=j$ that
satisfy (49) and (50).

######  Proposition 6.

Define $\varepsilon_{n}=(\log n)^{1/2}\Delta_{n}^{2}n^{-\alpha/(1+2\alpha)}$.
If the learning rate satisfies $\omega_{n}=c\Delta_{n}^{-2}$ for some
$0<c<1/2$ and Assumption 6 holds, then the Gibbs posterior distribution
concentrates at rate $\varepsilon_{n}$ with respect to
$d(\theta,\theta^{\star})=\|\theta-\theta^{\star}\|_{L_{2}(P)}$.

Since the mathematical statement does not give sufficient emphasis to the
adaptation feature, we should follow up on this point. That is,
$n^{-\alpha/(2\alpha+1)}$ is the optimal rate (Shen and Ghosal,, 2015) for
estimating an $\alpha$-Hölder smooth function, and it is not difficult to
construct an estimator that achieves this, at least approximately, if $\alpha$
is known. However, $\alpha$ is unknown in virtually all practical situations,
so it is desirable for an estimator, Gibbs posterior, etc. to adapt to the
unknown $\alpha$. The concentration rate result in Proposition 6 says that the
Gibbs posterior achieves nearly the optimal rate adaptively in the sense that
it concentrates at nearly the optimal rate as if $\alpha$ were known.

The concentration rate in Proposition 6 depends on the complexity of the
parameter space as determined by $\Delta_{n}$ in Assumption 6(3). For example,
if the sup-norm bound on $\theta^{\star}$ were known, then $\Delta_{n}$ and
the learning rate $\omega_{n}$ could be taken as constants and the rate would
be optimal up to a $(\log n)^{1/2}$ factor. On the other hand, if greater
complexity is allowed, e.g., $\Delta_{n}=(\log n)^{p}$ for some power $p>0$,
then the concentration rate takes on an additional $(\log n)^{2p}$ factor,
which is not a serious concern.

##  6 Application: personalized MCID

###  6.1 Problem formulation

In the medical sciences, physicians who investigate the efficacy of new
treatments are challenged to determine both statistically and practically
significant effects. In many applications some quantitative effectiveness
score can be used for assessing the statistical significance of the treatment,
but physicians are increasingly interested also in patients’ qualitative
assessments of whether they believed the treatment was effective. The aim of
the approach described below is to find the cutoff on the effectiveness score
scale that best separates patients by their reported outcomes. That cutoff
value is called the minimum clinically important difference, or MCID.

In this section we follow up on the MCID problem discussed in Syring and
Martin, (2017) with a covariate-adjusted, or _personalized_ , versions. The
point is that, in medicine, there is a trend away from the classical “one size
fits all” treatment procedures, to treatments that are tailored more-or-less
to each individual. Along these lines, naturally, doctors would be interested
to understand how that threshold for practical significance depends on the
individual, hence there is interest in a so-called personalized MCID. Some key
references include Hedayat et al., (2015) and Zhou et al., (2020).

###  6.2 Smooth personalized MCID

Let the data $U^{n}=(U_{1},\ldots,U_{n})$ be iid $P$, where each observation
is a triple $U_{i}=(X_{i},Y_{i},Z_{i})$ denoting the patient’s diagnostic
measurement, their self-reported effectiveness outcome $Y_{i}\in\\{-1,1\\}$,
and covariate value $Z_{i}\in\mathbb{Z}\subseteq\mathbb{R}^{q}$, for
$i=1,\ldots,n$ and $q\geq 1$. In practice the diagnostic measurement has
something to do with the effectiveness of the treatment so one can imagine
examples including blood pressure, blood glucose level, and viral load.
Examples of covariates include a patient’s age, weight, and gender. The idea
is that the $x$-scale cutoff for practical significance would depend on the
covariate $z$, hence the MCID is a function, say, $\theta(z)$, and the goal is
to learn this function.

The true MCID $\theta^{\star}$ is defined as the solution to an optimization
problem. That is, if

$\ell_{\theta}(x,y,z)=\tfrac{1}{2}[1-y\,\mathrm{sign}\\{x-\theta(z)\\}],\quad(x,y,z)\in\mathbb{R}\times\\{-1,1\\}\times\mathbb{Z},$
(52)

then the expected loss is $R(\theta)=P[Y\neq\mathrm{sign}\\{X-\theta(Z)\\}]$,
and the true MCID function is defined as the minimizer
$\theta^{\star}=\arg\min_{\theta\in\Theta}R(\theta)$, where the minimum is
taken over a class $\Theta$ of functions on $\mathbb{Z}$. Alternatively, as in
Section 5.5, the true $\theta^{\star}$ satisfies $\eta_{z}(x)>\tfrac{1}{2}$ if
$x>\theta^{\star}(z)$ and $\eta_{z}(x)\leq\tfrac{1}{2}$ if
$x\leq\theta^{\star}(z)$, where $\eta_{z}(x)=P(Y=1\mid X=x,\,Z=z)$ is the
conditional probability function.

###  6.3 Gibbs posterior construction

As described in Section 2, the Gibbs posterior distribution is based on an
empirical risk function which, in the present case, is given by

$R_{n}(\theta)=\frac{1}{2n}\sum_{i=1}^{n}[1-Y_{i}\,\mathrm{sign}\\{X_{i}-\theta(Z_{i})\\}],\quad\theta\in\Theta.$
(53)

In order to put this theory into practice, it is necessary to give the
function space $\Theta$ a lower-dimensional parametrization. In particular, we
consider a true MCID function $\theta^{\star}$ belonging to a Hölder class as
in Assumption 4 but with unknown smoothness, as in Section 5.6. And, we model
$\theta^{\star}$ by a linear combination of basis functions
$\theta(z)=\theta_{J,\beta}(z):=\textstyle\sum_{j=1}^{J}\beta_{j}f_{j}(z)$,
for basis functions $f_{j}$, $j=1,\ldots,J$. Then, each $\theta$ is identified
by a pair $(J,\beta)$ consisting of a positive integer $J$ and a
$J$-dimensional coefficient vector $\beta$. We use cubic B-splines in the
numerical examples in Section 6.5, but any basis capable of approximating
$\theta^{\star}$ will work, and see (44).

The prior setup is similar to the setup in Section 5.6 above in (48). That is,
the prior is specified hierarchically with a marginal prior $\pi$ for $J$ and
a suitable conditional prior $\widetilde{\Pi}_{J}$ for $\beta_{J}$, given $J$.
And, as mentioned before, very simple choices of the marginal and conditional
priors achieve the desired adaptive rates.

###  6.4 Concentration rate result

Assumption 7 below concerns the smoothness of $\theta^{\star}$ and requires
the chosen basis satisfy the approximation property used previously; it also
refers to the same mild assumptions on random series priors used in Section
5.6 sufficient to ensure adequate prior mass is assigned to a neighborhood of
$\theta^{\star}$; finally, it assumes a margin condition on the classifier
like that used in Section 5.5.1 and Assumption 5a(5). These conditions are
sufficient to establish a Gibbs posterior concentration rate.

######  Assumption 7.

  1. 1.

The true MCID function $\theta^{\star}:\mathbb{Z}\to\mathbb{R}$ for a compact
subset $\mathbb{Z}$ of $\mathbb{R}$ and $\theta^{\star}$ is Hölder smooth with
parameters $(\alpha,L)$ (see Assumption 4.1);

  2. 2.

the basis functions satisfy the approximation property in (44);

  3. 3.

the prior distribution $\Pi$ for $\theta$ is defined hierarchically as in (48)
with marginal and conditional priors $(\pi,\widetilde{\Pi})$ for $J$ and
$\beta_{J}$ given $J=j$ that satisfy (49) and (50); and,

  4. 4.

there exists $h\in(0,1)$ such that $P\\{|2\eta_{Z}(X)-1|\leq h\\}=0$; and,

  5. 5.

the conditional distribution, $P_{z}$, of $X$, given $Z=z$, has a density with
respect to Lebesgue measure that is uniformly bounded away from infinity.

######  Proposition 7.

Suppose Assumption 7 holds, with $\alpha$ as defined there, and set
$\varepsilon_{n}=(\log n)n^{-\alpha/(1+\alpha)}$. For any fixed $\omega>0$ the
Gibbs posterior distribution concentrates at rate $\varepsilon_{n}$ with
respect to the divergence

$\displaystyle d(\theta,\theta^{\star})$
$\displaystyle=P\\{\theta(Z)\wedge\theta^{\star}(Z)\leq
X\leq\theta(Z)\vee\theta^{\star}(Z)\\}$
$\displaystyle=\int_{\mathbb{Z}}\int_{\theta(z)\wedge\theta^{\star}(z)}^{\theta(z)\vee\theta^{\star}(z)}P_{z}(dx)\,P(dz).$
(54)

###### Proof.

See Appendix D. ∎

The Gibbs posterior distribution we have defined for the personalized MCID
function achieves the concentration rate in Proposition 7 adaptively to the
unknown smoothness $\alpha$ of $\theta^{\star}$. Mammen and Tsybakov, (1995)
consider estimation of the boundary curve of a set, and they show that the
minimax optimal rate is $n^{-\alpha/(\alpha+1)}$ when the boundary curve is
$\alpha$-Hölder smooth and distance is measured by the Lebesgue measure of the
set symmetric difference. In our case, if $(X,Z)$ has a joint density, bounded
away from 0, then our divergence measure $d(\theta,\theta^{\star})$ is
equivalent to

$\text{Leb}(\\{(x,z):x\leq\theta(z)\\}\,\triangle\,\\{(x,z):x\leq\theta^{\star}(z)\\}),$

in which case our rate is within a logarithmic factor of the minimax optimal
rate.

Hedayat et al., (2015) also study the personalized MCID and derive a
convergence rate for an M-estimator of $\theta^{\star}$ based on a smoothed
and penalized version of the misclassification error loss function we defined
in Section 6.2. It is difficult to compare our result with theirs, for
instance, because their rate depends on two user-controlled sequences related
to the smoothing and penalization of their loss. But, as mentioned above, our
rate is near optimal in certain cases, so the asymptotic results in Hedayat et
al., (2015) cannot be appreciably better than our rate in Proposition 7.

###  6.5 Numerical illustrations

We performed two simulation examples to investigate the performance of the
Gibbs posterior for the personalized MCID. In both examples we use a constant
learning rate $\omega=1$, but we generally recommend data-driven learning
rates; and see Syring and Martin, (2019).

For the first example we sample $n=100$ independent observations of $(X,Y,Z)$.
The covariate $Z$ is sampled from a uniform distribution on the interval
$[0,3]$. Given $Z=z$, the diagnostic measure $X$ is sampled from a normal
distribution with mean $z^{3}-3z^{2}+5$ and variance $1$, and the patient-
reported outcome $Y$ is sampled from a Rademacher distribution with
probability

$\eta_{z}(x)=\begin{cases}\Phi(x;z^{3}-3z^{2}+5-0.05,1/2),&x>z^{3}-3z^{2}+5\\\
\Phi(x;z^{3}-3z^{2}+5+0.05,1/2),&x\leq z^{3}-3z^{2}+5,\end{cases}$ (55)

where $\Phi(x;\mu,\sigma)$ denotes the ${\sf N}(\mu,\sigma)$ distribution
function. The addition of $\pm 0.05$ in the formula of $\eta_{z}(x)$ is to
meet the margin condition in Assumption 7.4. As mentioned above, we
parametrize the MCID function by piecewise polynomials, specifically, cubic
B-splines. For highly varying MCID functions, a reversible-jump MCMC algorithm
that allows for changing numbers of and break points in the piecewise
polynomials may be helpful; see Syring and Martin, (2020) for such an
approach. However, for this example we fix the parameter dimension to just six
B-spline functions, which allows us to use a simple Metropolis–Hastings
algorithm to sample from the Gibbs posterior distribution. Since the dimension
is fixed, the prior is only needed for the B-spline coefficients, and for
these we use diffuse independent normal priors with mean zero and standard
deviation of $6$. Over $250$ replications, the average empirical
misclassification rate is 16% using the Gibbs posterior mean MCID function
compared to 13% using the true MCID function when applying these two
classifiers to a hold-out sample of $100$ data points.

The left pane of Figure 1 shows the results for one simulated data set under
the above formulation. Even with only $n=100$ samples, the Gibbs posterior
does a good job of centering on the true MCID function. The right pane
displays the pointwise Gibbs posterior mean MCID function for each of $250$
repetitions of simulation 1, along with the overall pointwise mean of these
functions, and the true MCID function. The Gibbs posterior mean function is
stable across repetitions of the simulation.

Figure 1: Left: the posterior mean function (dashed), true MCID function
(solid), and data ($Y=1$ black points, $Y=-1$ gray points) for one replication
of the first simulation. Right: the Gibbs posterior mean MCID functions (solid
gray) for each of $250$ repetitions of the first simulation, the overall mean
function across those repetitions (dashed black), and the true MCID function
(solid black).

The second example we consider includes a vector covariate similar to that in
Example 1 of Scenario 2 in Hedayat et al., (2015). We sample $n=1000$
independent observations of $(X,Y,Z)$, where $Z=(Z_{1},Z_{2})$ has a uniform
distribution on the square $[0,3]^{2}$. Given $Z=z$, the diagnostic measure
$X$ has a normal distribution with mean $z_{1}+2z_{2}$ and variance $1$, and
the patient-reported outcome $Y$ is a Rademacher random variable with
probability

$\eta_{z}(x)=\begin{cases}\Phi(x;z_{1}+2z_{2}-0.05,1),&x>z_{1}+2z_{2}\\\
\Phi(x;z_{1}+2z_{2}+0.05,1),&x\leq z_{1}+2z_{2}.\end{cases}$ (56)

In practice it is common to have more than one covariate, so this second
example is perhaps more realistic than the first. However, it is much more
difficult to visualize the MCID function for more than one covariate, so we do
not display any figures for this example. We use tensor product B-splines with
$8$ fixed B-spline functions ($16$ coefficients) to parametrize the MCID
function. Again, we use independent diffuse normal priors with zero mean and
standard deviation equal to $6$ for each coefficient. Over $100$ repetitions
of this simulation we observe an average empirical misclassification rate of
24% using the Gibbs posterior mean MCID function compared to 23% using the
true MCID function when applied to a hold-out sample of $1000$ data points.

##  7 Concluding remarks

In this paper we focus on developing some simple, yet general, techniques for
establishing various asymptotic concentration properties for Gibbs posteriors.
A key take-away is that the robustness to model misspecification offered by
the Gibbs framework does not come at the expense of slower concentration
rates. Indeed, in the examples presented here—and others presented elsewhere,
e.g., Syring and Martin, (2020)—the rates achieved are the same as those
achieved by traditional Bayesian posteriors and (nearly) minimax optimal.
Another main point is that Gibbs posterior distributions are not inherently
challenging to analyze; on the contrary, the proofs presented herein are
concise and transparent. An additional novelty to the analysis presented here
is that we consider cases where the learning rate can be non-constant, i.e.,
either a vanishing sequence or data-dependent, and prove corresponding
posterior concentration rate results.

While a number of examples were presented herein, the results should be
applicable to many other situations. In particular, problems concerning high-
dimensional parameters subject to sparsity constraints are very common in the
recent literature, and the Gibbs posterior framework seems well-suited to
their study. Recent work on estimating high-dimensional sparse precision
matrices has advanced the understanding of priors and/or penalties sufficient
for obtaining the best concentration rates (Gan et al.,, 2019). Many versions
of this problem minimize a convex loss, so the ideas presented in Sections 5.1
and 5.4 may be helpful for developing Gibbs posterior distributions for these
parameters. In a number of statistical models involving stochastic processes
the likelihood is not available, so estimation of unknown parameters has to
proceed using other means. For example, M-estimators have recently been
formulated to estimate the so-called Hurst parameter in processes modeled by
tempered fractional Brownian motion (Boniece et al.,, 2016), and the Gibbs
framework could be applied there in order to incorporate prior knowledge about
the Hurst parameter into the data analysis.

Here the focus has been on deriving Gibbs posteriors with the best possible
concentration rates, and selection of learning rates has proceeded with these
asymptotic properties in mind. In other works (see, e.g., Syring and Martin,,
2019), random learning rates are derived for good uncertainty quantification
in finite-samples. We conjecture, however, that the two are not mutually-
exclusive, and that learning rates arising as solutions to the calibration
algorithm in Syring and Martin, (2019) also have the desirable concentration
rate properties. Proving this conjecture seems challenging, but Section 3.3
provides a first step in this direction.

Additionally, it is reasonable to expect Condition 1 and its variations may be
valuable in the study of M-estimators. Recently, we have made some connections
between these conditions and the properties of hypothesis tests constructed
from empirical risk functions, and work tying these topics together is already
underway.

## Acknowledgments

The authors thank the anonymous reviewers for their helpful feedback on a
previous version of the manuscript. This work is partially supported by the
U.S. National Science Foundation, DMS–1811802.

## References

  * Alquier, (2008) Alquier, P. (2008).  PAC-Bayesian bounds for randomized empirical risk minimizers.  Math. Methods Statist. 17(4):279–304.
  * Alquier et al., (2016) Alquier, P., Ridgway, J., and Chopin, N. (2016).  On the properties of variational approximations of Gibbs posteriors.  J. Mach. Learn. Res. 17:1–41.
  * Barron et al., (1999) Barron, A., Schervish, M. J., and Wasserman, L. (1999).  The consistency of posterior distributions in nonparametric problems.  Ann. Statist. 27(2):536–561.
  * Bhattacharya and Martin, (2020) Bhattacharya, I., and Martin, R. (2020).  Gibbs posterior inference on multivariate quantiles.  arXiv:2002.01052.
  * Bissiri et al., (2016) Bissiri, P.G., Holmes, C.C., and Walker, S.G. (2016).  A general framework for updating belief distributions.  J. R. Stat. Soc. Ser. B. Stat. Methodol. 78:1103–1130.
  * Boniece et al., (2016) Boniece, B. C., Didier, G., and Sabzikar, F. (2020+).  Tempered fractional Brownian motion: wavelet estimation, modeling and testing.  To appear in Appl. Comput. Harmon. Anal.
  * Boucheron et al., (2012) Boucheron, S., Lugosi, G., and Massart, P. (2012).  Concentration Inequalities: A Nonasymptotic Theory of Independence. Clarendon Press, Oxford.