# Energy dissipation in high speed impact on granular media

Manish Kumar Mandal1, Saikat Roy1, Corresponding author:
<EMAIL_ADDRESS>1 Department of Chemical Engineering, Indian Institute
of Technology Ropar, Rupnagar, Punjab, India 140001

###### Abstract

In this work, we thoroughly investigate the impact process on the granular
media in the limit when the ratio of the impact velocity to the acoustic speed
becomes of the order of $0.01$-$1$, which is far greater than the existing
literature ($0.0001-0.001$). We show that the energy dissipation is largely
due to the energy cost associated with the exploration between different
metastable states via large scale reorganization of the force chain network.
In this regime, the conventional drag force models break down, and the drag
force can not be decomposed into a depth dependent static pressure and a depth
independent inertial drag as proposed in the existing literature. The high
dynamical stress generates acoustic pulses, which propagate longer distances
rather than decaying exponentially, as observed in the previous works. In the
latter stage of the impact process, the boundary also plays an essential role
in the reorganization of the force chains as the reflected acoustic pulses
interact with the original impact pulses. Furthermore, we study the scaling of
the early stage peak forces with the impact velocity and find that spatial
dimensionality strongly influences the scaling.

##  I Introduction

Granular material is a special class of complex systems composed of many
interacting constituents behaving collectively. On top of the inherent complex
behaviour of granular material, the response of granular media under high
speed impact is notoriously difficult to model analytically since the impact
process never reaches a steady state (steady-state velocity different from
zero) except when a heavy intruder hits a superlight granular media Pacheco-
Vázquez  _et al._ (2011). The existing constitutive laws for the granular
materials are applicable only for the steady, fully developed flow conditions,
whereas the granular impact process leads to unsteady and complex flow.
Scarcity of governing equations and insufficient force and trajectory data at
the level of grain led to the development of numerous phenomenological
modelsCiamarra  _et al._ (2004); de Bruyn and Walsh (2004); Lohse  _et al._
(2004); Hou  _et al._ (2005); Ambroso  _et al._ (2005); Goldman and Umbanhowar
(2008) to describe the scaling of crater morphology, collision time,
penetration depth with that of impact velocity and grain and intruder
properties.

Following these works, various researchers Allen  _et al._ (1957); Katsuragi
and Durian (2007); Umbanhowar and Goldman (2010) proposed that the granular
drag force term, $F_{d}$ can be decomposed into a static depth dependent
friction term and a velocity dependent inertial drag term. Although the drag
model was shown to be valid in the past works, but almost all of the
experimental and simulational studies in the literature focused on the impact
velocities, $V_{0}$ ($1$ to $5$ $m/s$) that are far below than the velocity
scale set by the acoustic speed, $V_{a}$ ($2000$ to $5000$ $m/s$) in the same
media. Recently, Clark  _et al._ (2012, 2015) cleverly reduced the stiffness
of the grains to bring down the force propagation speed and consequently made
the low impact velocity approach the force propagation speed. The nature of
force propagation was shown to depend on a dimensionless parameter, $B$, which
is the ratio of the collision time scale ($t_{col}$) and the time scale set by
the intruder impact velocity ($D_{p}/V_{0}$), where $D_{p}$ is the grain
diameter and the collision time, $t_{col}$ can be calculated based on the
interaction law Clark  _et al._ (2015); SM . Impact pulses propagate through
the sparse force chains when $B\approx 0.1$, whereas $B\to 1$ leads to a dense
space filling network with a homogeneous front. Although this study made some
interesting observations, but it does not address the real scenario where the
high speed impact requires greater amount of energy to be dissipated, and the
energy dissipation mechanism can be quite different compared to the low speed
impact. Very recent investigation Krizou and Clark (2020) suggests some
universal scaling of the early stage peak forces with the impact velocity, and
the scaling turns out to be insensitive to the spatial dimension and many
other system parameters. This is a puzzling observation and begs for a
detailed study on the nature of the initial forces during an impact in
granular media. The process of the granular impact and crater formation has
rich physics with wider application in many disciplines like
ballisticsForrestal and Luk (1992); Glößner  _et al._ (2017),
astrophysicsRobbins and Hynek (2012), wind-blown transport of sands via
granular splashing Kamath  _et al._ (2022) and earth sciencesMelosh (1989). At
present, there is very scant literature available on the high speed
(comparable to the force propagation speed) impact due to the technological
limitations. The nature of the drag forces and the energy dissipation
mechanism are completely unknown in this regime. In this work, we employ
extensive numerical simulations (both in $2D$ and $3D$) to comprehend the
physics of the impact process in the high speed limit. The applicability of
the existing drag force models is also tested for a wide range of
$V_{0}/V_{a}$, from $0.008$ to $0.25$. The mechanism of energy transfer and
its eventual dissipation in granular media during high speed impact is
unveiled via spatio-temporal monitoring of the displacement field, velocity
field and complex force chain networks. The scaling of the early stage peak
forces with the impact velocity is thoroughly investigated. Also, the effect
of the boundary and its importance in transmitting or holding the impact
stress is explored in detail.

##  II Simulation methodology

###  II.1 Contact interaction

Frictional granular material is used as a model system for studying the high-
speed impact cratering. Discrete element method (DEM ) simulation is employed
to keep track of the particles with frictional interaction taking into account
both the normal($F_{n}$) and the path dependent tangential forces($F_{t}$).
Simulation is performed both in two and three dimensions. Open source codes,
Large-scale Atomic/Molecular Massively Parallel Simulator[LAMMPS] Plimpton
(1995); Silbert  _et al._ (2001); Brilliantov  _et al._ (1996) are used and
customized to carry out the numerical simulation. In 2D as well as 3D the
particle-particle interactions are modeled as linear spring-dashpot model with
a velocity dependent damping and static friction. We have also used non-linear
Hertzian interaction between particles in $3D$. Implementation of static
friction is done through the tracking of the elastic part of the shear
displacement from the time contact was first made. Particles $i$ and $j$, with
position vectors given by ${{\underline{r}}_{i},{\underline{r}}_{j}}$, have
linear velocities ${{\underline{v}}_{i},{\underline{v}}_{j}}$ and angular
velocities ${{\underline{\omega}}_{i},{\underline{\omega}}_{j}}$ respectively.
Grains will experience a normal force, ${\underline{F}}^{(n)}_{ij}$ whenever
there is a relative normal compression on contact given by
$\Delta_{ij}=|{\underline{r}}_{ij}-D_{ij}|$, where ${\underline{r}}_{ij}$
denotes the vector joining the centers of mass and $D_{ij}=R_{i}+R_{j}$ with
$R_{i}$ and $R_{j}$ being radii of particles. The normal force is modeled as a
Hookean spring like interaction, whereas the tangential force is given by
similar linear elastic relation upto the sliding pointBandi  _et al._ (2018) .
The force magnitudes are given as,

${\underline{F}}^{(n)}_{ij}=k_{n}\Delta_{ij}{\underline{n}}_{ij}-\frac{\gamma_{n}}{2}{\underline{v}}_{n_{ij}}$
(1)
${\underline{F}}^{(t)}_{ij}=-k_{t}{\underline{t}}_{ij}-\frac{\gamma_{t}}{2}{\underline{v}}_{t_{ij}}$
(2)

where $\Delta_{ij}$ and $t_{ij}$ denote normal and tangential displacements
respectively; ${\underline{n}}_{ij}$ denotes the normal unit vector given by
${\underline{r}}_{ij}/|{\underline{r}}_{ij}|$. $k_{n}$ and $k_{t}$ are
respectively stiffness of the springs for the normal and tangential mode of
elastic displacement. For the Hertzian case, contact normal force is given as,
$F_{Hertzian}=F_{Hookean}\sqrt{\Delta_{ij}}\sqrt{\frac{R_{i}R_{j}}{R_{i}+R_{j}}}$.
Viscoelastic damping constant for normal and tangential deformation are
denoted by $\gamma_{n}$ and $\gamma_{t}$ respectively and
${\underline{v_{n}}}_{ij}$ as well as ${\underline{v_{t}}}_{ij}$ designate the
normal and tangential component of the relative velocity between two grains.
The relative normal and tangential velocity are given as:

$\displaystyle{\underline{v}}_{n_{ij}}$ $\displaystyle=$
$\displaystyle({\underline{v}}_{ij}.{\underline{n}}_{ij}){\underline{n}}_{ij}$
(3) $\displaystyle{\underline{v}}_{t_{ij}}$ $\displaystyle=$
$\displaystyle{\underline{v}}_{ij}-{\underline{v}}_{n_{ij}}-\frac{1}{2}({\underline{\omega}}_{i}+{\underline{\omega}}_{j})\times{\underline{r}}_{ij}.$
(4)

where ${\underline{v}}_{ij}={\underline{v}}_{i}-{\underline{v}}_{j}$. Elastic
tangential displacement ${\underline{t}}_{ij}$ is set to zero when the contact
develops for the first time between two particles and is computed using
$\frac{d{\underline{t}}_{ij}}{dt}={\underline{v}}_{t_{ij}}$ . The simulation
also account for the rigid body rotation around the contact point to make sure
that ${\underline{t}}_{ij}$ always remains in the local tangential plane of
the contact. The gravitational forces are also accounted for in the
simulation. The translational and rotational degrees of freedom of the
particles are computed using Newton’s second law; total forces and torques on
particle $i$ are given as:

$\displaystyle{\underline{F}}^{(tot)}_{i}$ $\displaystyle=$ $\displaystyle
m_{i}{\underline{g}}+\sum_{j}\left({\underline{F}}^{(n)}_{ij}+{\underline{F}}^{(t)}_{ij}\right)$
(5) $\displaystyle{\underline{\tau}}^{(tot)}_{i}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{j}{\underline{r}}^{ij}\times{\underline{F}}^{(t)}_{ij}.$
(6)

Note that, the tangential force follows a linear relationship with the
relative tangential displacement at the contact point as long as the
tangential force is below the limit set by the Coulomb friction,

$F^{(t)}_{ij}\leq\mu F^{(n)}_{ij}\ ,$ (7)

where $\mu$ stands for friction coefficient. Upon exceeding this, the contact
slips in a dissipative fashion and the tangential displacement is truncated
accordingly to satisfy the Coulomb criterion. Simulation also incorporates the
effect of inelastic collision for both normal and tangential mode of relative
movement via the viscoelastic damping coefficients ($\gamma_{n,t}$) which are
related to the coefficient of restitution($\epsilon_{n,t}$) and the collision
time as below:

$\epsilon_{n,t}=exp(-\gamma_{n,t}t_{col}/2),$ (8)

where the collision time $t_{col}$ is given as:

$t_{col}=\pi(2k_{n}/m-\gamma_{n}^{2}/4)^{-0.5}$ (9)

In our simulation, $\epsilon_{n,t}$ is taken asSilbert  _et al._ (2001) $0.9$
since the coefficient of restitution for dry sand falls in the similar range.
In order to capture the dynamics at the time scale of collision, the time step
is set as $t_{col}/50$. $t_{col}$ was calculated for the simulation parameters
shown in Table 1.

Table 1: Simulation parameters used for both 2D and 3D simulations $k_{n}$ | $\gamma_{n}$ | $\frac{k_{t}}{k_{n}}$ | $\frac{\gamma_{t}}{\gamma_{n}}$ | $\mu$ | $\epsilon$ | $g$
---|---|---|---|---|---|---
$2\times 10^{5}$ | $600$ | $2/7$ | $1$ | $0.5$ | $0.90$ | $9.8$

In Table 1 parameters are shown only for a reference simulation, we also
varied the system parameters for different simulations and describe the same
in the text whenever the parameter values change with respect to the reference
values. Mass per unit area for the smallest particle (in $2D$ bi-disperse
particles are used) $m_{g}$ is $0.133$. Accordingly, the acoustic speed based
on the properties of the smallest diameter particle is approximately
$1200$($V_{a}=\sqrt{\frac{k_{n}}{m_{g}}}$). Also, parameter,
$B$($=\frac{t_{col}V_{0}}{D_{p}}$) can be calculated by estimating the
collision time based on the initial impact velocity, $V_{0}$ and the form of
the potential ($F=k_{n}\Delta^{\beta}$) assuming no viscous dissipation. Here,
$D_{p}$ is the diameter of the smallest particle and $\beta=1$ and $1.5$ for
Hookean and Hertzian, respectively. Consequently, $B$ is given as
$P(\beta)\left(\frac{V_{0}}{V_{a}}\right)^{\left(\frac{2}{\beta+1}\right)}$,
where $P(\beta)$ is as followsClark  _et al._ (2015),

$P(\beta)=(\pi(\beta+1)/16)^{\frac{1}{\beta+1}}\frac{4\sqrt{\pi}\Gamma(1+\frac{1}{\beta+1})}{\Gamma(\frac{1}{2}+\frac{1}{\beta+1})}$
(10)

For the Hookean case, $P(\beta)=3.9374$. Note that all the parameter values
are reported in SI units.

###  II.2 Bed preparation

Figure 1: Left: Stable granular bed in two dimension, Right: Corresponding
force chain network before the impact

For $3D$ simulation of high speed impact, we first make a three dimensional
box with periodic boundary conditions in the $x$ and $y$ direction, and a
fixed wall at the bottom plane whose outward unit normal is $\widehat{e_{z}}$.
The box has a dimension of $30D_{p}\times 30D_{p}\times 500D_{p}$. The length
in the $z$ direction is kept long enough to prevent atom loss due to the high
speed impact. Next, $N=50000$ (larger system size $N=250000$ is also
investigated) mono-disperse spherical particles having diameter, $D_{p}=0.2m$
are dropped under gravity and following that, the bed is allowed sufficient
time to attain the mechanically stable state with an average force and torque
on each particle of the order of $10^{-9}$ with negligible kinetic energy.
This stabilization process is extremely important to ensure that we start with
a stable bed rather than some fragile configuration which may lead to spurious
results. After a stable bed is created, the impacting ball of spherical shape
of diameter $10D_{p}$ is placed very close to the free surface and is given an
initial impact velocity in the $z$ direction which is also the direction of
the gravity. The average stable bed height is $46D_{p}$ and the corresponding
volume fraction, $\phi$ of the prepared bed is $0.62$, which is very close to
the random close packingScott and Kilgour (1969) state.
We also performed simulation in two dimension because visualization of the
grain scale phenomena is much easier in $2D$ compared to higher dimensions.
Similar to $3D$, we first define a simulation box of dimension
$100D_{p1}\times 1000D_{p1}$ in the $xy$ plane, followed by pouring of $10000$
(we also simulate a large system,$N=40000$) bi-disperse granular particles,
half of which is having diameter, $D_{p1}=0.2m$ and the other half has
diameter, $D_{p2}=0.28m$, into the box under gravity. Selection of bi-disperse
particles was done to prevent crystallization which is spontaneous for mono-
disperse particles in two dimension. The particles were also poured layer by
layer instead of pouring them at one go to achieve a stable configuration with
sufficient mechanical equilibrium. Pouring all the particle at one go would
have caused large collisional stress leading to very long computational time
before the system can be relaxed. Following the similar analysis as done on
$3D$ bed, the bed height in $2D$ is calculated to be $139D_{p1}$ having
packing fraction of $0.83$. In Fig. 1 , representative schematic of the
initial $2D$ granular bed is shown with and without the force chains. The
force chain figure is the visual representation of contact forces between
particles giving us an idea about the gradient of the stress created due to
the gravity. The contact force data was used to draw the force chain of a
certain thickness that scales with the magnitude of the force. As expected,
the force chains are more dense and thick at the bottom due to the
gravitational stress whereas the force chain network is very sparse close to
the surface. Force chains represent the stress transmission paths in the
granular media and it plays a crucial role in the propagation of any
disturbance through the granular materials.

##  III Results and discussions

Figure 2: (a) Temporal variation of drag force in $2D$ for
$V_{0}/V_{a}$=$0.03$, $0.04$, $0.06$ and $0.08$: inset shows the corresponding
variation of velocity as a function of the intruder depth. (b) Drag force in
$3D$ for the similar range of velocities, inset: intruder velocity vs intruder
depth. (c) Scaling of the dimensionless peak force, $F_{max}/k_{n}D$ with the
non-dimensional impact velocity, $V_{0}/V_{a}$ in $2D$ for different stiffness
and system size. The solid line shows a power law fit with an exponent $\sim
1.33$. (d) The same scaling in $3D$ for different system size, interaction law
and stiffness. The solid fit corresponds to a power law scaling with an
exponent $\sim 1.5$. (e) and (f) show the force chain evolution (for Hookean
interaction) in the initial stages of the impact for $B=0.01$ and $B=0.1$,
respectively. Here, the line thickness is scaled according to the contact
force magnitude normalized by the mean force.

Figure 3: (a) Displacement field near the bottom wall showing the arrival of
the compression front. Here, $V_{0}=4.4$ m/s and displacement magnitude is
magnified 40 times the original. (b) Elastic unloading of the bed leading to
the reversal of the displacement field.(c)-(d) Temporal variation of the force
on the wall and the intruder in $3D$ are plotted simultaneously for
$V_{0}=4.4$ and $40$ m/s, respectively. We also observe a similar behavior in
$2D$ (not shown).

###  III.1 Phenomenology of impact

We begin our study by analyzing the force on the intruder exerted by the
granular media in Fig. 2 (a)-(b) for both $2D$ and $3D$. In line with the low
speed impactKrizou and Clark (2020), the drag force attains a maximum very
quickly, followed by slow relaxations dominated by fluctuations. The
fluctuations decay very fast in the three dimension compared to the two
dimensional case since the extra dimension gives the system additional
direction to relax the effect of the impact. Also, the response in $3D$ in
terms of the peak force becomes stronger than $2D$ with the increase in impact
velocity. Fig. 2 (c)-(d) depicts the power law scaling of the peak force with
$V_{0}$. Interestingly, the exponent in $3D$ is higher ($\sim 1.5$) than in
2D, contradicting the recent results where the scaling exponent ($\sim 1.33$)
is independent of the spatial dimension. Since the volume of the
configurational space of stress paths (force chains) is greater in $3D$ than
$2D$, the number of acoustic events or pulses that carry the intruder energy
into the medium also increases proportionately in three dimension compared to
the two-dimensional case. Hence, greater resistance to impact is observed in
$3D$ than in $2D$ for comparable initial impact speeds. To explicitly show
that the observed scaling is not an artifact of the finite size, we vary the
system size in both two and three dimension and find that the scaling exponent
remains invariant of the system size. We have also checked in all our
simulations that the peak force occurs well before the acoustic pulses get
reflected from the boundary (Please see supplemental(SM, ) video 1 showing the
pulse propagation in a large system, $N=40000$). The exponent of the power law
scaling is also insensitive to the interaction potential, and the stiffness,
in line with the recent observations Krizou and Clark (2020).

Note that during the attainment of the peak force, the velocity does not
reduce significantly and the intruder hardly penetrates. After reaching the
peak, the retarding force starts to relax and the ball starts penetrating
significantly into the media with “stop and go” kind of motion where at some
moments the ball is falling freely under gravity even inside the bed and this
phenomenon is observed for all ranges of impact speed. The stronger
fluctuations in $2D$ is also reflected in the velocity-depth trajectories (see
the inset of Fig. 2(a) and (b)). The shape of the velocity-depth curve is
concave upward, which presents striking dissimilarity with the existing
literature in the low speed limit Huang  _et al._ (2020); Ambroso  _et al._
(2005); Umbanhowar and Goldman (2010), where the shape is concave downward and
can be reproduced by solving the conventional drag models Huang  _et al._
(2020); Katsuragi and Durian (2007); Pacheco-Vázquez  _et al._ (2011). This
observation is evocative of the possible breakdown of the known macroscopic
drag force models in the high speed limit and presents the possibility of
unexplored rich grain scale physics.

###  III.2 Grain scale picture

Before we test the drag force models explicitly, we turn our attention to the
grain scale picture of the impact process in terms of the spatio-temporal
variations of the complex force networks, displacement field and the velocity
field. Despite, the recent experiments Clark  _et al._ (2012, 2015) with
photoelastic disks presented some interesting grain scale picture of the
impact process, the exact and the complete understanding of the force network
evolution and its effect on the intruder motion is not well understood.
Photoelastic measurements have a resolution of $256\times 584$ pixels Clark
_et al._ (2012) at high speed and thus give only a measure of the total
photoelastic intensity in an image. In contrast, simulations can provide
better insights into the nature of the vectorial contact forces. In Fig. 2
(e)-(f), we show the evolution of the force chain networks in the early stages
of the impact for the two values of $B=0.01$ and $B=0.1$. For both cases,
before the impact, gravity sets the gradient of the pressure for which the
force chains look denser at the bottom. As soon as the intruder strikes the
bed, the large dynamical stress dictates the gradient of pressure, and the
force chains look denser close to the impact point. The impact energy gets
propagated in the form of acoustic pulses, which reach the end of the system
boundary even though the system size is large enough to avoid boundary effects
(see the Supplemental videos 2-3SM ). This very fast large length scale
propagation of disturbances is evocative of the collective motion of the
granular particles, which are correlated upto long range even before the
impact.

We observe the reflection of the acoustic pulses from the boundary and the
sideways scattering and branching of the pulses. Reflected pulses also
interact with the original pulses emitted from the intruder and give rise to
continuous large scale reorganization of the force chains, which results in
the temporally fluctuating force on the intruder. None of the existing
experimental studies capture such long range propagation of disturbances due
to the low resolution of the photoelastic measurements that are typically used
to characterize the forces in the experiments. The existing literature,
without any physical explanation, suggests an exponential decay of the pulses
meaning the pulses decay almost immediately after traveling only a few
particle diameter, which is at variance with our simulation observations. We
observe that the force propagation happens via a well-defined compression
front for $B=0.01$, which is far below that reported ($B=0.6$) in the recent
$2D$ experiments Clark  _et al._ (2015). Even for $V_{0}/V_{a}=0.03$
($B=0.1$), we see a dense compression front propagating through the media
whereas the previous observations showed sparse chain-like force propagation
for the same value of $B$. We speculate that the setup used in the experiment
had strong side-wall friction and boundary effects which led to the quick
damping of the energy pulses. It is also possible that the limited resolution
of the photoelastic response at low stress levels makes the determination of
the signal propagation far beneath the intruder difficult.

We also simultaneously analyze the particle displacement field near the bottom
(see Fig.3 (a)-(b) and also supplemental video 4 SM ) and observe that a
compression front indeed reaches the bottom, and a strong elastic resistance
is provided by the bottom wall leading to the flip in the particle
displacement field. During the whole process, the particles are moving
cooperatively, and the flipping of the displacement field takes longer than
the collision time, meaning a large length scale reorganization is inevitable.
Intriguingly, the phenomenon of compression and decompression keeps repeating
until the intruder comes to rest, which also gives rise to large fluctuations
in the force time series. We also show the velocity field at different stages
of the impact in the supplemental videos 5-6SM  for low and high impact
velocities, and the long distance propagation of an acoustic pulse is vividly
observed. The disturbance propagation speed can also be estimated by
monitoring the force on the wall (See Fig. 3 (c)-(d)) and measuring the time
taken for the disturbance to reach the wall (time of flight measurement). The
force propagation speed is almost of the order of the acoustic speed
($1200m/s$), and it is independent of the impact velocity for the linear
interaction. Surprisingly, the temporal variation of the force on the wall
looks very similar to the force-time series of the intruder, albeit with
lesser fluctuations since the wall is in contact with large number of force
chains. As the wall and the intruder are far apart, a similar temporal
response at distant points suggests that the large scale reorganization of the
force chain networks dictates the response.

###  III.3 Force network reorganization, dissipation and breakdown of inertial
drag models

Figure 4: (a) Temporal variation of the force anisotropy, $a_{n}$ (Right
$y$-axis) is plotted simultaneously with the temporal variation of the force,
$F_{i}$ on the intruder (Left $y$-axis) for $B=0.2$: The inset shows the
change in the orientations of the force network as a function of time (b)
Similar to (a) except $B=0.5$.

For the quantitative description of the force chain reorganization, we now
monitor the anisotropy of the force network and its preferred orientations. As
earlier studies showed that friction does not significantly influence the
dynamic impact process, we focus on the force skeletons formed by the contact
normal forces only. The normal force anisotropy and its preferred direction
are defined by $a_{n}$ and $\theta_{f}$, respectively. The calculation of
these parameters from the discrete simulation data is performed by introducing
a second order tensor,
$\xi_{ij}\approx\frac{1}{N_{g}}\sum_{\theta_{g}}\bar{f_{n}}n_{i}n_{j}$, where
$N_{g}$ denotes the number of orientation intervals spanning from $0$ to
$2\pi$, $\theta_{g}$ is the average orientation of a group and the
corresponding average normal force is denoted by $\bar{f_{n}}$, $n_{i}$
denotes the Cartesian components of the contact unit normal vector. Anisotropy
parameters, $a_{n}$ and $\theta_{f}$ are related to the invariants of
$\xi_{ij}$ and its principal directions:
$a_{n}=\frac{2\sqrt{(\xi_{11}-\xi_{22})^{2}+4\xi_{12}^{2}}}{\xi_{11}+\xi_{22}};\;tan\,2\theta_{f}=\frac{2\xi_{12}}{\xi_{11}-\xi_{22}}$.

Figure 5: (a) Net acceleration, $a+g$ versus square of the velocity, $V^{2}$
is plotted at seven fixed depths($d$) for different initial impact velocities
($V_{0}/V_{a}=0.008$ to $0.25$) in $2D$ , (b) the same is plotted for the
three dimensional case for impact velocities ranging from $V_{0}/V_{a}=0.02$
to $0.13$.

Fig. 4 depicts the temporal variation of both the anisotropy and its principal
direction, along with the force time series of the intruder for different
impact velocities. Before the impact, the force chains are organized mostly in
the direction of gravity and $a_{n}=0.2$. Upon impact, the force chains start
to reorganize (see the inset of Fig. 4) as reflected by the change in the
principal direction ($\theta_{f}$) of the force network. Also, the force
network becomes progressively anisotropic to support the sudden impact load
and hits a peak, which in turn gives rise to a maximum force on the intruder.
After the peak, the force anisotropy decreases in a manner similar to the
decrease in the force on the intruder. In the later stages of the relaxation,
temporal variation of the force anisotropy decorrelates from the force-time
series of the intruder. Furthermore, we find a strong correlation between the
orientation of the normal force network and the temporal evolution of the
force on the intruder. We also observe a time lag between the force network
reorganization and its effect to be felt on the intruder force. The time lag
decreases with the impact velocity, suggestive of a decreasing length scale
upto which the reorganization occurs. In summary, under high speed impact, the
granular media constantly traverses between different fragile states via large
scale reorientation of the force networks, and the force on the intruder is
the consequence of this large scale reorganization. These transient
rearrangements of the force network lead to plastic dissipation and are the
principal energy loss mechanism during the high speed impact. Finally, we test
the validity of the existing drag models Katsuragi and Durian (2007); Allen
_et al._ (1957) in both $2D$ and $3D$ by monitoring the net acceleration and
the speed of the intruder at different fixed depths for different impact
velocities. If the depth-independent inertial drag were to apply to our high
speed regime, net acceleration would be quadratic in speed resulting in
parallel straight lines when $a+g$ is plotted against $V^{2}$ for different
depths. Fig. 5 instead presents an entirely contrasting picture; depth-
dependent quadratic profiles are obtained when net acceleration is plotted
against the square of the speed at a fixed depth for different trajectories.
Hence, the conventional depth independent inertial drag models are unable to
capture the force on the intruder in the high speed limit; rather we observe
that the net acceleration varies linearly with the velocity with a depth
dependent slope, though this scaling needs to be checked extensively with
large data set.

##  IV Conclusion

In summary, we employed large scale numerical simulations to understand the
response of the granular media under a high speed impact. Although a large
volume of the work on granular impact exists in the literature, but most of
the approaches to tackle such a complex problem are heuristic with
insufficient grain scale understanding of the highly dynamic impact
phenomenon. This work presents a detailed microscopic length scale picture of
the impact process in terms of the evolution of the inhomogeneous force chain
networks, displacement field, and the velocity field as the impact progresses.
These particle scale information proved to be quite useful in understanding
the dissipation mechanism in the granular materials which are neither solid
nor fluid. Contrary to the previous works showing the exponential spatial
decay of the acoustic pulses, we vividly demonstrate a large-length scale
propagation of disturbances that get reflected from the boundary, interfering
with the original pulses. These acoustic pulses, in turn, induce large scale
reorganization of the force chain network, and the granular media constantly
explores different fragile states to support the impact load. Reorientation of
the force chains leads to plastic dissipation and the eventual absorption of
the impact energy. The large scale temporal evolution of the force chain
networks dictates the force on the intruder. Consequently, this novel energy
dissipation picture does not corroborate with the conventional drag models and
hence, the breakdown of the depth independent inertial drag forces.
Furthermore, the power law scaling of the early stage peak forces with the
impact velocity shows a dependence on the spatial dimensionality, which is at
variance with the past works. The result of this work begs for the development
of a novel theoretical framework to explain the drag force on the intruder in
the high speed limit. It would also be interesting to study the effect of
cohesive interactions on the different aspects of the impact
processRalaiarisoa  _et al._ (2022) and the scaling of the peak forces since,
in a natural setup, attractive forces are expected to be present due to van
der Waals forces, humidity, moisture, etc.

Acknowledgement-S.R. acknowledges the support of SERB under Grant No.
SRG/2020/001943 and the IIT Ropar under ISIRD grant.

## References

  * Pacheco-Vázquez  _et al._ (2011) F. Pacheco-Vázquez, G. Caballero-Robledo, J. Solano-Altamirano, E. Altshuler, A. Batista-Leyva, and J. Ruiz-Suárez, Physical review letters 106, 218001 (2011).
  * Ciamarra  _et al._ (2004) M. P. Ciamarra, A. H. Lara, A. T. Lee, D. I. Goldman, I. Vishik, and H. L. Swinney, Physical review letters 92, 194301 (2004).
  * de Bruyn and Walsh (2004) J. R. de Bruyn and A. M. Walsh, Canadian Journal of Physics 82, 439 (2004).
  * Lohse  _et al._ (2004) D. Lohse, R. Bergmann, R. Mikkelsen, C. Zeilstra, D. Van Der Meer, M. Versluis, K. Van Der Weele, M. van der Hoef, and H. Kuipers, Physical review letters 93, 198003 (2004).
  * Hou  _et al._ (2005) M. Hou, Z. Peng, R. Liu, K. Lu, and C. Chan, Physical review E 72, 062301 (2005).
  * Ambroso  _et al._ (2005) M. Ambroso, R. D. Kamien, and D. J. Durian, Physical Review E 72, 041305 (2005).
  * Goldman and Umbanhowar (2008) D. I. Goldman and P. Umbanhowar, Physical Review E 77, 021308 (2008).
  * Allen  _et al._ (1957) W. A. Allen, E. B. Mayfield, and H. L. Morrison, Journal of Applied Physics 28, 370 (1957).
  * Katsuragi and Durian (2007) H. Katsuragi and D. J. Durian, Nature physics 3, 420 (2007).
  * Umbanhowar and Goldman (2010) P. Umbanhowar and D. I. Goldman, Physical review E 82, 010301 (2010).
  * Clark  _et al._ (2012) A. H. Clark, L. Kondic, and R. P. Behringer, Physical review letters 109, 238302 (2012).
  * Clark  _et al._ (2015) A. H. Clark, A. J. Petersen, L. Kondic, and R. P. Behringer, Physical review letters 114, 144502 (2015).
  * (13) See Supplemental Material at for additional information .
  * Krizou and Clark (2020) N. Krizou and A. H. Clark, Physical Review Letters 124, 178002 (2020).
  * Forrestal and Luk (1992) M. Forrestal and V. Luk, International Journal of Impact Engineering 12, 427 (1992).
  * Glößner  _et al._ (2017) C. Glößner, S. Moser, R. Külls, S. Heß, S. Nau, M. Salk, D. Penumadu, and N. Petrinic, Experimental Mechanics 57, 261 (2017).
  * Robbins and Hynek (2012) S. J. Robbins and B. M. Hynek, Journal of Geophysical Research: Planets 117 (2012).
  * Kamath  _et al._ (2022) S. Kamath, Y. Shao, and E. J. Parteli, Geophysical Research Letters 49, e2022GL097767 (2022).
  * Melosh (1989) H. J. Melosh, New York: Oxford University Press; Oxford: Clarendon Press (1989).
  * Plimpton (1995) S. Plimpton, Journal of computational physics 117, 1 (1995).
  * Silbert  _et al._ (2001) L. E. Silbert, D. Ertaş, G. S. Grest, T. C. Halsey, D. Levine, and S. J. Plimpton, Physical Review E 64, 051302 (2001).
  * Brilliantov  _et al._ (1996) N. V. Brilliantov, F. Spahn, J.-M. Hertzsch, and T. Pöschel, Physical review E 53, 5382 (1996).
  * Bandi  _et al._ (2018) M. Bandi, H. G. E. Hentschel, I. Procaccia, S. Roy, and J. Zylberg, EPL (Europhysics Letters) 122, 38003 (2018).
  * Scott and Kilgour (1969) G. D. Scott and D. M. Kilgour, Journal of Physics D: Applied Physics 2, 863 (1969).
  * Huang  _et al._ (2020) K. Huang, D. Hernández-Delfin, F. Rech, V. Dichtl, and R. C. Hidalgo, Scientific reports 10, 1 (2020).
  * Ralaiarisoa  _et al._ (2022) V. Ralaiarisoa, P. Dupont, A. O. El Moctar, F. Naaim-Bouvet, L. Oger, and A. Valance, Physical Review E 105, 054902 (2022).


# A comparison of learning rate selection methods in generalized Bayesian
inference

Pei-Shien Wu111Department of Statistics, North Carolina State University;
<EMAIL_ADDRESS><EMAIL_ADDRESS>and Ryan Martin∗

###### Abstract

Generalized Bayes posterior distributions are formed by putting a fractional
power on the likelihood before combining with the prior via Bayes’s formula.
This fractional power, which is often viewed as a remedy for potential model
misspecification bias, is called the learning rate, and a number of data-
driven learning rate selection methods have been proposed in the recent
literature. Each of these proposals has a different focus, a different target
they aim to achieve, which makes them difficult to compare. In this paper, we
provide a direct head-to-head comparison of these learning rate selection
methods in various misspecified model scenarios, in terms of several relevant
metrics, in particular, coverage probability of the generalized Bayes credible
regions. In some examples all the methods perform well, while in others the
misspecification is too severe to be overcome, but we find that the so-called
generalized posterior calibration algorithm tends to outperform the others in
terms of credible region coverage probability.

_Keywords and phrases:_ coverage probability; generalized posterior
calibration algorithm; model misspecification; SafeBayes algorithm.

##  1 Introduction

Specification of a sound model is a critical part of an effective statistical
analysis. This is especially true for a Bayesian approach, since the
statistical model or likelihood is explicitly used to construct the posterior
distribution from which inferences will be drawn. However, it is common in
applications to know relatively little about the phenomenon under
investigation, which impacts our ability to specify a sound statistical model.
For this reason, the effects of model misspecification have received
considerable attention; in the Bayesian literature, this includes Berk,
(1966), Bunke and Milhaud, (1998), Diaconis and Freedman, 1986b ; Diaconis and
Freedman, 1986a , Kleijn and van der Vaart, (2006, 2012), Walker, (2013), De
Blasi and Walker, (2013), Ramamoorthi et al., (2015), and Grünwald and van
Ommen, (2017). In the most general case, misspecification implies that there
is no “true” parameter value that the posterior could concentrate around.
Instead, under suitable conditions, the posterior will concentrate around a
“best” parameter value, one that minimizes the Kullback–Leibler divergence of
the posited model from the true data-generating distribution. But even in
those relatively nice cases, where the best parameter value around which the
posterior concentrates could be meaningful, or even equal to the real quantity
of interest, there is reason for concern. Kleijn and van der Vaart, (2012)
showed that misspecification can also affect the posterior spread, which means
that the actual frequentist coverage probability of the Bayesian posterior
credible region can be arbitrarily far below the advertised/nominal level.

Real coverage probabilities differing significantly from advertised levels is
a serious concern (Fraser, 2011; Martin, 2019). A gap between real and
advertised coverage probabilities can can have various causes, but here we
focus on model misspecification. To avoid this misspecification bias, there
are a few options: first, to take an approach that does not depend explicitly
on a statistical model; second, to work with a model that is sufficiently
broad that misspecification is virtually impossible; and third, to make some
adjustments to correct for potential model misspecification. From a Bayesian
perspective, the first fix is not available, since Bayes’s formula requires a
likelihood function. The second fix amounts to the use of Bayesian
nonparametrics but, when the quantity of interest is a low-dimensional feature
of the full distribution, introducing an infinite-dimensional parameter, with
the computational and statistical challenges that entails, would be overkill.
This leaves only the third option, but what kind of adjustments might the
Bayesian consider? Recently, Grünwald and van Ommen, (2017) argued that a
certain adjustment to the usual Bayesian posterior distribution could repair
inconsistencies resulting from model misspecification. The goal of the present
paper is to investigate the extent to which Grünwald and van Ommen’s proposed
adjustment—and other related adjustments in the literature—can close the gap
between the real and advertised coverage probabilities of Bayesian posterior
credible regions affected by model misspecification.

More specifically, here we will be working in the so-called generalized Bayes
framework, which differs from the traditional Bayes framework only in that a
learning rate parameter, a power $\eta>0$ on the likelihood function, is
introduced. That is, if we have data $D^{n}$ and a posited statistical model
$P_{\theta}^{n}$, indexed by a parameter $\theta$ in $\Theta$, then the
generalized Bayes posterior distribution for $\theta$ is

$\Pi_{n}^{(\eta)}(d\theta)\propto
L_{n}^{\eta}(\theta)\,\Pi(d\theta),\quad\theta\in\Theta,$ (1)

where $\theta\mapsto L_{n}(\theta)=L(\theta;D^{n})$ is the likelihood function
and $\Pi$ is a prior distribution on $\Theta$. Among the first papers to adopt
such an approach is Walker and Hjort, (2001), followed up on by Zhang, (2006).
Bissiri et al., (2016) showed that (1) is the principled way to update prior
beliefs when the model is potentially misspecified, and that the appearance of
a non-trivial learning rate is a necessary by-product. A different connection
between robustness and learning rate $\eta<1$ was made recently in Miller and
Dunson, (2019). More details about model misspecification and generalized
Bayes are given in Section 2.

Grünwald and van Ommen’s claim is that, for a sufficiently small learning rate
$\eta$, certain model misspecification biases can be repaired. Of course, the
threshold defining “sufficiently small” cannot be known in practice, so some
data-driven choices are required. Grünwald, (2012, 2018), Grünwald and van
Ommen, (2017), and Heide et al., (2020) developed a so-called SafeBayes
algorithm to choose the learning rate $\eta$, based on minimizing a sequential
risk measure. A number of other learning rate selection methods have been
proposed recently, including the two distinct information matching strategies
in Holmes and Walker, (2017) and Lyddon et al., (2019), and bootstrap-
motivated calibration methods of Syring and Martin, (2019). Since the role
played by the learning rate is relatively unfamiliar and since the various
methods differ significantly in terms of their motivations and
implementations, it would be beneficial to see a head-to-head comparison in
terms of some standard metrics, for example, the validity and efficiency of
the corresponding generalized Bayes credible regions. This paper aims to fill
this gap.

The remainder of this paper is organized as follows. In Section 2, we discuss
the behavior of Bayesian posterior distribution under a misspecified model and
define and review the literature on generalized Bayes posteriors. For the
latter, the choice of learning rate is essential, and we provide details for
four recently proposed learning rate selection methods in Section 3. Then, in
Section 4, we show a simple illustrative example to give some intuition about
how the different methods perform and, in particular, this suggests that the
methods which are not designed specifically to calibrate the credible region’s
coverage probability may not be able to achieve the nominal level in general.
Simulation results are presented in Sections 5–6, for linear and binary
regression models, and the take-away message is that the method of Syring and
Martin, (2019) is more stable in achieving the coverage probability across
different sample sizes and misspecification degrees compared to the others.
Some concluding remarks are given in Section 7.

##  2 Background

###  2.1 Model misspecification

Suppose we have data $D^{n}$ which, for simplicity, we assume consists of
independent and identically distributed observations: either response
variables $Y_{i}$ only or predictor and response variables pairs
$(X_{i},Y_{i})$, $i=1,\ldots,n$. To analyze these data, we posit a statistical
model $\mathscr{P}=\\{P_{\theta}:\theta\in\Theta\\}$, a collection of
probability measures on the sample space, indexed by a parameter $\theta$
taking values in the parameter space $\Theta$. From this model and the
observed $D^{n}$, we obtain a likelihood function $L_{n}$. The likelihood
summarizes the information in the data relative to the posited model, which
can be combined with prior information encoded in a distribution $\Pi$ for
$\theta$ on $\Theta$ via Bayes’s formula:

$\Pi_{n}(d\theta)\propto L_{n}(\theta)\,\Pi(d\theta),\quad\theta\in\Theta.$

In the Bayesian paradigm, inferences about $\theta$ are drawn based on the
posterior distribution $\Pi_{n}$, e.g., degrees of belief about the
truthfulness of an assertion “$\theta\in A$,” for $A\subset\Theta$, are
summarized by the posterior probability $\Pi_{n}(A)$.

Let $P^{\star}$ denote the true distribution of $Y_{1}$ or of $(X_{1},Y_{1})$.
If the model is correctly specified, then there exists a
$\theta^{\star}\in\Theta$ such that $P^{\star}=P_{\theta^{\star}}$. In that
case, under suitable regularity conditions, inference based on the posterior
distribution will be valid, at least asymptotically. That is, $\Pi_{n}$ will
concentrate its mass around $\theta$ as $n\to\infty$ and, moreover, the
Bernstein–von Mises theorem (e.g., van der Vaart, 2000, Theorem ??), states
that $\Pi_{n}$ is approximately a normal distribution, centered at the maximum
likelihood estimator $\hat{\theta}_{n}$, with covariance matrix proportional
to the inverse of the Fisher information matrix at $\theta^{\star}$. This
implies, among other things, that credible regions derived from $\Pi_{n}$
closely resemble those asymptotic confidence regions based on likelihood
theory. Therefore, asymptotically, the Bayesian posterior credible regions
will have frequentist coverage probability close to the advertised level.

If the model is incorrectly specified, in the sense that
$P^{\star}\not\in\mathscr{P}$, then there are several challenges. First, there
is no “true” $\theta^{\star}$, which creates some challenges in
interpretation. Indeed, the maximum likelihood estimator, Bayes posterior, or
any other model-based procedure will identify the Kullback–Leibler projection
of $P^{\star}$ onto the model, i.e.,

$\theta^{\dagger}=\arg\min_{\theta}K(P^{\star},P_{\theta}),$

where $K(P^{\star},P_{\theta})=\int\log(dP^{\star}/dP_{\theta})\,dP^{\star}$
is the Kullback–Leibler divergence of $P_{\theta}$ from $P^{\star}$. In
general, $\theta^{\dagger}$ does not have a real-life interpretation but, in
some cases, certain features of $P^{\star}$ can be identified based on a
misspecified model. For example, if $\mathscr{P}$ is an exponential family,
then the mean function of the exponential family model, evaluated at
$\theta^{\dagger}$, equals the mean of $P^{\star}$ (Bunke and Milhaud, 1998,
Example 2). Another similar case is considered in Section 5. The second
challenge is that, even in the case where $\theta^{\dagger}$ has a (limited)
real-life interpretation, misspecification can still negatively impact
posterior inferences. Kleijn and van der Vaart, (2012) established a version
of the Bernstein–von Mises theorem under model misspecification which states
that, under certain regularity conditions, the posterior $\Pi_{n}$ will be
approximately normal, with mean equal to the maximum likelihood (or M-)
estimator $\hat{\theta}$ and covariance matrix $V_{\theta^{\dagger}}^{-1}$,
where

$V_{\theta^{\dagger}}=\int\Bigl{(}\frac{\partial^{2}\log
p_{\theta}}{\partial\theta\partial\theta^{\top}}\Bigr{|}_{\theta=\theta^{\dagger}}\Bigr{)}\,dP^{\star},$

and $p_{\theta}$ is the density function corresponding to $P_{\theta}$. The
problem, of course, is that $V_{\theta^{\dagger}}^{-1}$ is not the asymptotic
covariance matrix of $\hat{\theta}_{n}$; the latter, as shown by Huber, (1967)
and van der Vaart, (2000), has the famous sandwich formula
$V_{\theta^{\dagger}}^{-1}\Lambda_{\theta^{\dagger}}V_{\theta^{\dagger}}^{-1}$,
where

$\Lambda_{\theta^{\dagger}}=\int\Bigl{(}\frac{\partial\log
p_{\theta}}{\partial\theta}\Bigr{|}_{\theta=\theta^{\dagger}}\Bigr{)}\Bigl{(}\frac{\partial\log
p_{\theta}}{\partial\theta}\Bigr{|}_{\theta=\theta^{\dagger}}\Bigr{)}^{\top}\,dP^{\star}.$

The implication of this covariance mismatch is that, even if the quantity of
interest can be identified under the misspecified model, the frequentist
coverage probability of the Bayes posterior credible sets could be arbitrarily
far from the advertised level. The question is: can something be done to
correct this problematic behavior?

###  2.2 Generalized Bayes

Modifying the usual Bayesian update with a learning rate $\eta$ as in (1) is a
simple change, but it has some unexpected consequences. In particular, Walker
and Hjort, (2001) showed that, for a correctly specified model, consistency of
the generalized Bayes posterior $\Pi_{n}^{(\eta)}$ in (1) could be established
for any $\eta<1$, with only local conditions on the prior—as opposed to the
local and global conditions required for consistency with $\eta=1$ (e.g.,
Ghosal et al., 1999; Barron et al., 1999). The intuition given by Walker et
al., (2005) is that inconsistencies result from the posterior over-fitting or
tracking the data too closely, and the fractional power discounts the data
slight to prevent this over-fitting. The Walker–Hjort result has been extended
to cover posterior concentration rates, where the removal of the global prior
conditions—usually formulated in terms of metric entropy (cf. Ghosal et al.,
2000; Ghosal and van der Vaart, 2017)—leads to simpler proofs and generally
(at least slightly) faster rates. See Zhang, (2006) for one of the first
papers exploring these ideas, and Bhattacharya et al., (2019) and Grünwald and
Mehta, (2020) for more recent contributions. The fractional power has also
been employed recently in work on high-dimensional problems using an empirical
or data-driven prior (e.g., Martin and Walker, 2019; Martin et al., 2017;
Martin and Tang, 2020) where, again, the fractional power is motivated by the
desire to prevent over-fitting; see, also, Martin, (2017) and Martin and Ning,
(2020) for some potential benefits of $\eta<1$ to uncertainty quantification.

When the model is misspecified, however, the learning rate is less about
convenience and more about necessity. Bissiri et al., (2016) showed that the
generalized Bayes update (1) is fundamental from a decision-theoretic point of
view. Moreover, they argue that the learning rate $\eta$ naturally emerges
since, roughly, the parameter $\theta^{\dagger}$ being estimated is defined by
minimizing the expectation of a loss function $\theta\mapsto\log p_{\theta}$,
and since that minimization problem is invariant to scalar multiples of the
loss, the learning rate should appear in the posterior (1). In fact, the loss
function interpretation makes their result much more general. In many cases,
it is more natural to formulate the inference problem with a loss function
rather than a statistical model. These are often referred to as Gibbs
posterior distributions; see Syring and Martin, (2017, 2019); Syring and
Martin, 2020b ; Syring and Martin, 2020a , Bhattacharya and Martin, (2020),
and Wang and Martin, (2020).

Beyond recognizing the importance of the learning rate parameter, an actual
value for $\eta$ needs to be set in practical applications. Several recent
papers—including Grünwald and van Ommen, (2017), Holmes and Walker, (2017),
Lyddon et al., (2019), and Syring and Martin, (2019)—have proposed data-driven
choices for the learning rate, with different motivations. Section 3 describes
these methods. The remainder of the paper is focused on a comparison of these
different learning rate methods.

##  3 Learning rate selection methods

###  3.1 Grünwald’s SafeBayes

Grünwald and van Ommen, (2017) observe that, when the model is non-convex and
misspecified, there is a chance for hyper-compression. This is the term they
use to describe the seemingly paradoxical result that the Bayesian predictive
distribution can be closer, in a Kullback–Leibler sense, to the true
$P^{\star}$ then the within-model Kullback–Leibler minimizer
$P_{\theta^{\dagger}}$. What makes this possible is non-convexity: the
predictive distribution is an average of in-model distributions $P_{\theta}$
which, without convexity, could end up outside the model and potentially
closer to $P^{\star}$ than is $P_{\theta^{\dagger}}$. Besides being counter-
intuitive, hyper-compression also reveals a practical problem, namely,
inconsistency—that the posterior distribution is not concentrating its mass
near $\theta^{\dagger}$ as expected. To overcome this, Grünwald and van Ommen,
(2017) suggest to work with a new (hypothetical) model, with densities

$p_{\theta}^{(\eta)}(x,y)=p^{\star}(x,y)\bigl{\\{}p_{\theta}(y\mid
x)/p_{\theta^{\dagger}}(y\mid x)\bigr{\\}}^{\eta},$

indexed by a parameter $\eta>0$. We say this model is “hypothetical” because
it depends on $p^{\star}$ and $\theta^{\dagger}$, two ingredients that are not
available to the data analyst. However, if $\eta$ is sufficiently small, in
the sense that $\int p_{\theta}^{(\eta)}(x,y)\,dx\,dy$ is strictly less than
1, then this indeed defines a genuine statistical model, with two interesting
properties:

  * •

it is not misspecified, i.e., the Kullback–Leibler minimizer is
$\theta^{\dagger}$ and $p_{\theta^{\dagger}}^{(\eta)}=p^{\star}$;

  * •

and the Bayesian posterior based on this new model is precisely the
generalized Bayes posterior $\Pi_{n}^{(\eta)}$, with learning rate $\eta$, as
in (1).

Since this new model is not misspecified, hyper-compression and inconsistency
of the $\eta$-generalized Bayes posterior can be avoided. So: how to choose
$\eta$ sufficiently small?

Grünwald and van Ommen, (2017), building on work in, e.g., Grünwald, (2012),
argue that the so-called SafeBayes algorithm will select a learning rate
$\eta$ that is sufficiently small in the sense above. Define the cumulative
expected log-loss under the $\eta$-generalized Bayes posterior distribution,
as a function of $\eta$:

$\eta\mapsto\sum_{i=1}^{n}\int-\log p_{\theta}(Y_{i}\mid
X_{i})\,\Pi_{i-1}^{(\eta)}(d\theta).$

The SafeBayes algorithm returns the minimizer, $\hat{\eta}$, of this function
over the range $\eta\in[0,1]$. Grünwald, (2012) presents an argument for why
the SafeBayes choice of $\hat{\eta}$ works in the sense of being sufficiently
small as in the discussion above.

What we have described here is one of two versions of the SafeBayes algorithm
presented in Grünwald and van Ommen, (2017), namely, the “R-SafeBayes”
version. In our examples below, we found that the “R” version outperformed the
other—namely, the “I-SafeBayes” version—so here we only discuss the former.

###  3.2 Holmes and Walker, (2017)

Following Bissiri et al., (2016), the Bayesian and generalized Bayesian
frameworks can be considered simply as rules for using data to update prior
beliefs to posterior beliefs. As such, it makes sense to consider how much
information has been gained from the update, by comparing the prior to the
posterior. Of course, this information gain depends on both the updating rule
and on the data, and Holmes and Walker, (2017) proposed a procedure for
selecting the learning rate $\eta$ based on matching the expected information
gain between Bayes and generalized Bayes updates.

More formally, if $I_{\eta}(x,y)$ denotes the information gain in the
generalized Bayes update from prior to posterior based on learning rate $\eta$
and data values $(x,y)$, then Holmes and Walker, (2017) propose to set $\eta$
such that

$\int I_{\eta}(x,y)\,P^{\star}(dx,dy)=\int
I_{1}(x,y)\,P_{\theta^{\dagger}}(dx,dy),$ (2)

where $I_{1}(\cdot)$ denotes the information gain in the standard Bayesian
update. The specific choice of information measure they recommend is the
Fisher divergence

$I_{\eta}(x,y)=\int\\{\nabla\log\pi_{x,y}^{(\eta)}(\theta)-\nabla\log\pi(\theta)\\}^{2}\,\pi(\theta)\,d\theta,$

where $\pi_{x,y}^{(\eta)}$ denotes the generalized Bayes posterior based on
data $(x,y)$ and learning rate $\eta$, and $\nabla$ is the gradient operator.
Then it is straightforward to check that $I_{\eta}(x,y)=\eta^{2}I_{1}(x,y)$
and, therefore, by (2), an “oracle” learning rate is given by

$\eta^{\star}=\Bigl{\\{}\frac{\int
I_{1}(x,y)\,P_{\theta^{\dagger}}(dx,dy)}{\int
I_{1}(x,y)\,P^{\star}(dx,dy)}\Bigr{\\}}^{1/2}.$

Of course, both $P^{\star}$ and $P_{\theta^{\dagger}}$ are unknown, so
$\eta^{\star}$ cannot be evaluated, but the expectations can be estimated with
the actual data $\\{(X_{i},Y_{i}):i=1,\ldots,n\\}$. That is,

$\hat{\eta}=\Bigl{\\{}\frac{\int I_{1}(x,y)\,P_{\hat{\theta}_{n}}(dx,dy)}{\int
I_{1}(x,y)\,\mathbb{P}_{n}(dx,dy)}\Bigr{\\}}^{1/2},$

where $\hat{\theta}_{n}$ is the maximum likelihood estimator of $\theta$ under
the model—which is an estimate of $\theta^{\dagger}$—and $\mathbb{P}_{n}$ is
the empirical distribution of the data.

###  3.3 Lyddon et al., (2019)

The learning rate selection strategy presented in Lyddon et al., (2019) is
motivated by the weighted likelihood bootstrap approach of Newton and Raftery,
(1994), which was shown to generate bootstrap samples that have the same
asymptotic distribution as Bayesian posterior distribution under a correctly
specified model. For the case of a misspecified model, Lyddon et al., (2019)
proposed a modified the weighted likelihood bootstrap approach which replaces
the ordinary bootstrap with the Bayesian bootstrap, and establish its
asymptotic limiting distribution. Then following a strategy similar to that in
Holmes and Walker, (2017) described above, they propose to choose $\eta$ in
order to match the limiting $\eta$-generalized Bayes posterior to that of this
modified likelihood bootstrap. They then show that, using the notation defined
at the end of Section 2.1, an “oracle” learning rate is

$\eta^{\star}=\frac{\mathrm{tr}(V_{\theta^{\dagger}}\Lambda_{\theta^{\dagger}}^{-1}V_{\theta^{\dagger}})}{\mathrm{tr}(V_{\theta^{\dagger}})}.$

Again, since $\theta^{\dagger}$ and $P^{\star}$ are unknown, this oracle value
cannot be evaluated. However, a data-driven choice $\hat{\eta}$ can be
obtained by replacing $\theta^{\dagger}$ the maximum likelihood estimator and
the expectations with respect to $P^{\star}$ in $V_{\theta}$ and
$\Lambda_{\theta}$, respectively, with expectations with respect to the
empirical distribution $\mathbb{P}_{n}$.

###  3.4 Syring and Martin, (2019)

The three previous subsections describe principled learning rate selection
strategies, but none of those are tailored so that the generalized posterior
distribution achieves any specific and desirable properties. Since the
learning rate’s effect on the posterior is to control the spread, Syring and
Martin, (2019) proposed to tune the learning rate such that posterior credible
sets approximately achieve the nominal frequentist coverage probability.

The coverage probability function is given by

$c_{\alpha}(\eta\mid
P^{\star})=P^{\star}\\{C_{\alpha}^{(\eta)}(D^{n})\ni\theta^{\dagger}\\},$

where $C_{\alpha}^{(\eta)}$ is the $\eta$-generalized Bayes $100(1-\alpha)$%
credible region for $\theta$, e.g., a highest posterior density region,
$\theta^{\dagger}$ is the Kullback–Leibler minimizer in the model, treated as
a functional of $P^{\star}$, and $D^{n}=\\{(X_{i},Y_{i}):i=1,\ldots,n\\}$ is
the iid data set from $P^{\star}$. Then the goal is to find $\eta$ such that
$c_{\alpha}(\eta\mid P^{\star})=1-\alpha$. Of course, lots of the quantities
involved in this equation are unknown, but they can be estimated. In
particular, if $P^{\star}$ is replaced by the empirical distribution
$\mathbb{P}_{n}$, then the new equation is

$c_{\alpha}(\eta\mid\mathbb{P}_{n}):=\mathbb{P}_{n}\\{C_{\alpha}^{(\eta)}(D^{n})\ni\hat{\theta}_{n}\\}=1-\alpha,$

where $\hat{\theta}_{n}$ is the maximum likelihood estimator based on the
observed data, i.e., the “$\theta^{\dagger}$-functional” applied to
$\mathbb{P}_{n}$. Even this alternative coverage probability function cannot
be evaluated, since it requires enumerating all $n^{n}$ possible with-
replacement samples from the observed data, but a bootstrap approximation is
possible. That is, for $B$ bootstrap samples
$\tilde{D}_{1}^{n},\ldots,\tilde{D}_{B}^{n}$, calculate

$\hat{c}_{\alpha}(\eta\mid\mathbb{P}_{n})=\frac{1}{B}\sum_{b=1}^{B}1\\{C_{\alpha}^{(\eta)}(\tilde{D}_{b}^{n})\ni\hat{\theta}_{n}\\}.$

To solve the equation, $\hat{c}_{\alpha}(\eta\mid\mathbb{P}_{n})=1-\alpha$,
Syring and Martin, (2019) recommend a stochastic approximation scheme that
iteratively defines a learning rate sequence $(\eta_{t})$ as

$\eta_{t+1}=\eta_{t}+k_{t}\\{\hat{c}_{\alpha}(\eta\mid\mathbb{P}_{n})-(1-\alpha)\\},\quad
t\geq 1,$

where $k_{t}$ is a sequence such that $\sum_{t}k_{t}=\infty$ and
$\sum_{t}k_{t}^{2}<\infty$. When the $\eta_{t}$ sequence effectively
converges, the limit is the suggested learning rate $\hat{\eta}$. This is what
Syring and Martin, (2019) refer to as the generalized posterior calibration
(GPC) algorithm. Like SafeBayes, it is relatively expensive
computationally—these algorithms require posterior computations for multiple
learning rates and data sets—but the benefit is having a posterior
distribution with meaningful spread, even in finite samples.

##  4 Learning rates in a toy example

Before we consider the effect of different learning rate selection methods in
some non-trivial real-world problems, it helps to consider a simple example,
one where some of the calculations can be done by hand, to develop some
intuition about what to expect.

Suppose that the posited model for iid data $Y^{n}=(Y_{1},\ldots,Y_{n})$ is
$P_{\theta}={\sf N}(\theta,\sigma^{2})$, with $\sigma>0$ fixed, but that the
true distribution is $P^{\star}={\sf N}(\theta^{\star},\sigma^{\star 2})$,
where $\sigma^{\star}>0$ is potentially different from $\sigma$. It is easy to
confirm that the Kullback–Leibler minimizer satisfies
$\theta^{\dagger}=\theta^{\star}$, but the misspecified variance can still
cause problems, as we now demonstrate.

It is intuitively clear that the generalized Bayes framework could completely
resolve the model misspecification if the learning rate was chosen as
$\eta^{\star}=(\sigma/\sigma^{\star})^{2}$. More formally, Heide et al.,
(2020) show that, if the learning rate is no larger than this ratio, the
generalized Bayes posterior will enjoy fast root-$n$ rate convergence
properties, while Syring and Martin, (2019) argue that learning rate equal to
the ratio in order to achieve exact coverage of credible sets. In any case, of
course, one cannot make this learning rate choice in practice because it
depends on the unknown value of true variance. But this intuition tells us
what the different learning rate selection methods’ target should be.

To evaluate the performance of the different learning rate methods, we
simulate $1000$ data sets, for each of several different sample sizes $n$ and
values of $\eta^{\star}=(\sigma/\sigma^{\star})^{2}$, and compare the average
estimated learning rate against $\eta^{\star}$; see Figure 1. If the estimated
$\eta$ is close to the diagonal line $\eta=\eta^{\star}$, then the generalized
Bayesian credible sets have coverage probability near the nominal level. To
confirm this, see Figure 2. When the degree of misspecification is relatively
mild, all the methods perform well. As the misspecification degree increases,
however, we find that SafeBayes and the Holmes and Walker method have
decreasing coverage probability, quickly falling below any reasonable
tolerance. On the other hand, both the Lyddon et al. and Syring and Martin
methods are able to achieve the target 95% coverage probability over the
entire range of settings.

(a) $n=100$

(b) $n=200$

(c) $n=400$

(d) $n=800$

Figure 1: Average learning rate $\eta$, across 1000 replications, versus the
optimal $\eta^{\star}=(\sigma/\sigma^{\star})^{2}$, the closer to the diagonal
line the better. True (black), GPC (red), R-SafeBayes (blue), Holmes and
Walker (green), Lyddon et al. (orange).

(a) $n=100$

(b) $n=200$

(c) $n=400$

(d) $n=800$

Figure 2: Average coverage probability of the nominal $95\%$ generalized Bayes
credible intervals, across 1000 replications. GPC (red), SafeBayes (blue),
Holmes and Walker (green), Lyddon et al. (orange).

##  5 Learning rates in linear regression

###  5.1 Model setup

Consider a linear regression model of the form

$y_{i}=x_{i}^{\top}\beta+\sigma\varepsilon_{i},\quad i=1,\ldots,n,$ (3)

where the pairs $(x_{1},y_{1}),\ldots,(x_{n},y_{n})$, taking values in
$\mathbb{R}^{p}\times\mathbb{R}$, are independent, $\beta\in\mathbb{R}^{p}$ is
an unknown vector of coefficients, $\sigma$ is an unknown scale parameter, and
$\varepsilon_{1},\ldots,\varepsilon_{n}$ are random error terms. As is most
common, here we will consider a model that assumes the errors
$\varepsilon_{1},\ldots,\varepsilon_{n}$ are iid ${\sf N}(0,1)$, independent
of $x_{1},\ldots,x_{n}$. If it happens that the true distribution is different
from this posited model then, in general, we can expect an ordinary Bayes
posterior to suffer from misspecification bias. The goal here is to
investigate how the different learning rate methods can help the generalized
Bayes posterior to correct for this misspecification bias.

The two key assumptions behind the textbook linear regression model are that
the errors are (a) independent of covariates and (b) normally distributed.
Here we present results for two types of misspecification, namely, Dependent
Errors and Non-normal Errors. The specific form and degree of these
misspecifications will be described in the following subsections. For the
comparison, the metrics we consider are

  * •

mean value of the learning rate estimates, $\hat{\eta}$;

  * •

coverage probability of the $\hat{\eta}$-generalized Bayes 95% highest
posterior density credible sets for the full $\beta$ vector;

  * •

mean square error of the $\hat{\eta}$-generalized Bayes posterior mean of
$\beta$;

  * •

the average of the marginal $\hat{\eta}$-generalized Bayes posterior variances
for each coordinate of $\beta$.

We are specifically interested in the learning rate and its effect on the
coverage probability of the generalized posterior credible sets, so the first
two metrics are clear. The mean square error of the generalized posterior mean
acts like an overall measure of bias, i.e., how far does the center of the
posterior tend to be from the true parameter values. In the examples that
follow, we find that the mean square error does not vary much relative to the
learning rate selection method, which confirms our intuition that the learning
rate really only impacts the posterior spread. The fourth metric is an overall
measure of the spread of the generalized Bayes posterior, and we expect that
those learning rate selection methods whose credible regions tend to under-
cover will have smaller total variance.

###  5.2 Dependent errors

In the linear regression model described above, we sample the predictors
$x_{i}$ independently from a multivariate normal distribution with mean zero,
unit variance, and a first-order autocorrelation structure, i.e.,
$\mathsf{E}(x_{ij}x_{ik})=\rho^{|j-k|}$, with correlation $\rho=0.2$. We also
set the true coefficient vector to be $\beta=(1,1,2,-1)^{\top}$. The errors
$\varepsilon_{1},\ldots,\varepsilon_{n}$ are independent standard normal, but
the error standard deviation is misspecified in the sense that the presumed
constant $\sigma$ is not a constant, in particular, its value depends on the
individual $x_{i}$ in the following way. Let $\hat{\xi}_{0.05}$ and
$\hat{\xi}_{0.95}$ denote the sample $5^{\text{th}}$ and $95^{\text{th}}$
percentiles of $x_{11},\ldots,x_{n1}$. Then define the case-specific standard
deviation as

$\sigma_{i}=\begin{cases}s_{\text{small}}&\text{if
$x_{i1}<\hat{\xi}_{0.05}$}\\\ s_{\text{mod}}&\text{if $\hat{\xi}_{0.05}\leq
x_{i1}\leq\hat{\xi}_{0.95}$}\\\ 1&\text{if
$x_{i1}>\hat{\xi}_{0.95}$},\end{cases}$

where the small and moderate values, $s_{\text{small}}$ and $s_{\text{mod}}$,
control the degree of the departures from constant variance. We consider three
different degrees of misspecification.

Degree 1.

$s_{\text{small}}=0.25$ and $s_{\text{mod}}=0.50$;

Degree 2.

$s_{\text{small}}=0.05$ and $s_{\text{mod}}=0.25$;

Degree 3.

$s_{\text{small}}=0.01$ and $s_{\text{mod}}=0.10$.

A summary of the different learning rate selection procedures, across the
different misspecification degrees and sample sizes $n\in\\{100,200,400\\}$,
is presented in Table 1, based on 1000 data sets for each combination. In
Degree 1, where the misspecification is relatively mild, we see that all four
learning rate selections perform well and similarly in terms of both the
learning rates chosen—all near 1—and in the coverage probabilities. As
expected, however, as the misspecification gets more severe, in Degrees 2 and
3, the more disparity we see between the selected learning rates and, in turn,
in the coverage probabilities. Only GPC is able to achieve the nominal
coverage probability in the more severe misspecification settings, while the
performance of other methods can be quite poor, especially under Degree 3 with
small sample sizes. The mean square errors are more or less the same for the
methods within each sample size–degree combination; and the fact that these
values are small indicates the the posterior is generally centered around the
target $\beta$ values. As for the posterior spread, there is not much
difference between the results in the Degree 1 case with only mild
misspecification. However, in Degrees 2 and 3, where the misspecification is
more severe, we see greater difference in the posterior variance. As expected,
those methods whose posterior variance tends to be small are those who tend to
have credible sets that under-cover, in many cases severely.

Degree | $n$ | Method | $\hat{\eta}$ | Coverage | MSE | Variance
---|---|---|---|---|---|---
1 | $100$ | GPC | $0.95$ | $0.95$ | $0.05$ | $0.012$
|  | SafeBayes | $0.92$ | $0.94$ | $0.05$ | $0.014$
|  | Holmes and Walker | $1.00$ | $0.93$ | $0.05$ | $0.011$
|  | Lyddon et al. | $1.18$ | $0.89$ | $0.05$ | $0.010$
| $200$ | GPC | $0.95$ | $0.93$ | $0.02$ | $0.006$
|  | SafeBayes | $0.92$ | $0.93$ | $0.02$ | $0.007$
|  | Holmes and Walker | $0.99$ | $0.92$ | $0.02$ | $0.006$
|  | Lyddon et al. | $1.06$ | $0.90$ | $0.02$ | $0.005$
| $400$ | GPC | $0.94$ | $0.95$ | $0.01$ | $0.003$
|  | SafeBayes | $0.93$ | $0.94$ | $0.01$ | $0.003$
|  | Holmes and Walker | $0.99$ | $0.94$ | $0.01$ | $0.003$
|  | Lyddon et al. | $0.99$ | $0.94$ | $0.01$ | $0.003$
2 | $100$ | GPC | $0.79$ | $0.95$ | $0.06$ | $0.015$
|  | SafeBayes | $0.90$ | $0.90$ | $0.06$ | $0.014$
|  | Holmes and Walker | $0.98$ | $0.89$ | $0.06$ | $0.012$
|  | Lyddon et al. | $1.33$ | $0.76$ | $0.06$ | $0.009$
| $200$ | GPC | $0.75$ | $0.95$ | $0.03$ | $0.008$
|  | SafeBayes | $0.92$ | $0.90$ | $0.03$ | $0.006$
|  | Holmes and Walker | $0.97$ | $0.89$ | $0.03$ | $0.006$
|  | Lyddon et al. | $1.11$ | $0.84$ | $0.03$ | $0.005$
| $400$ | GPC | $0.74$ | $0.94$ | $0.01$ | $0.004$
|  | SafeBayes | $0.93$ | $0.89$ | $0.01$ | $0.003$
|  | Holmes and Walker | $0.96$ | $0.88$ | $0.01$ | $0.003$
|  | Lyddon et al. | $0.97$ | $0.88$ | $0.01$ | $0.003$
3 | $100$ | GPC | $0.54$ | $0.98$ | $0.07$ | $0.023$
|  | SafeBayes | $0.75$ | $0.87$ | $0.07$ | $0.018$
|  | Holmes and Walker | $0.94$ | $0.80$ | $0.07$ | $0.012$
|  | Lyddon et al. | $2.45$ | $0.38$ | $0.07$ | $0.005$
| $200$ | GPC | $0.53$ | $0.95$ | $0.04$ | $0.011$
|  | SafeBayes | $0.76$ | $0.86$ | $0.04$ | $0.008$
|  | Holmes and Walker | $0.91$ | $0.79$ | $0.04$ | $0.006$
|  | Lyddon et al. | $1.74$ | $0.53$ | $0.04$ | $0.003$
| $400$ | GPC | $0.53$ | $0.95$ | $0.02$ | $0.005$
|  | SafeBayes | $0.78$ | $0.84$ | $0.02$ | $0.004$
|  | Holmes and Walker | $0.89$ | $0.81$ | $0.02$ | $0.003$
|  | Lyddon et al. | $1.25$ | $0.69$ | $0.02$ | $0.002$
Table 1: Comparison of average learning rate estimates ($\hat{\eta}$),
estimated coverage probabilities (Coverage), mean square error (MSE), and
total posterior variance (Variance) across different sample sizes and
misspecification degrees in the Dependent Errors example.

###  5.3 Non-normal errors

Next, we consider departures from the specified model in terms of the
distribution of the error terms. It turns out that the performance of the
learning rate selection methods was less sensitive to departures from
normality compared to departures from the constant-error-variance assumption.
Here we present the results for only one kind of departure from normality,
namely, with heavy-tailed errors. In particular, consider errors
$\varepsilon_{1},\ldots,\varepsilon_{n}$ iid from a Student-t distribution
with degrees of freedom $\nu$. As before, we consider three degrees of
misspecification, each sufficiently light-tailed that the variance exists.

Degree 1.

$\nu=5$;

Degree 2.

$\nu=4$;

Degree 3.

$\nu=3$.

Table 2 summarizes the results just like in the previous subsection. Here,
however, the differences in performance across different learning rate
selection methods, sample sizes, and misspecification degrees is much smaller.
Overall the methods return similar learning rate estimates and hit the target
coverage probability on the mark. The method of Lyddon et al. tends to select
a learning rate that is too large, leading to under-coverage, but its
performance tends to improve as the sample size increases.

Degree | $n$ | Method | $\hat{\eta}$ | Coverage | MSE | Variance
---|---|---|---|---|---|---
1 | $100$ | GPC | $0.98$ | $0.95$ | $0.07$ | $0.019$
|  | SafeBayes | $0.90$ | $0.95$ | $0.07$ | $0.023$
|  | Holmes and Walker | $1.00$ | $0.94$ | $0.07$ | $0.019$
|  | Lyddon et al. | $1.28$ | $0.87$ | $0.07$ | $0.014$
| $200$ | GPC | $0.99$ | $0.96$ | $0.04$ | $0.009$
|  | SafeBayes | $0.90$ | $0.96$ | $0.04$ | $0.011$
|  | Holmes and Walker | $0.98$ | $0.96$ | $0.04$ | $0.009$
|  | Lyddon et al. | $1.15$ | $0.92$ | $0.04$ | $0.008$
| $400$ | GPC | $1.00$ | $0.95$ | $0.02$ | $0.005$
|  | SafeBayes | $0.92$ | $0.95$ | $0.02$ | $0.005$
|  | Holmes and Walker | $0.98$ | $0.95$ | $0.02$ | $0.005$
|  | Lyddon et al. | $1.07$ | $0.94$ | $0.02$ | $0.004$
2 | $100$ | GPC | $0.97$ | $0.96$ | $0.08$ | $0.024$
|  | SafeBayes | $0.89$ | $0.96$ | $0.08$ | $0.028$
|  | Holmes and Walker | $0.99$ | $0.96$ | $0.08$ | $0.023$
|  | Lyddon et al. | $1.38$ | $0.87$ | $0.08$ | $0.016$
| $200$ | GPC | $0.98$ | $0.96$ | $0.04$ | $0.011$
|  | SafeBayes | $0.91$ | $0.96$ | $0.04$ | $0.013$
|  | Holmes and Walker | $0.97$ | $0.96$ | $0.04$ | $0.011$
|  | Lyddon et al. | $1.20$ | $0.91$ | $0.04$ | $0.009$
| $400$ | GPC | $0.99$ | $0.96$ | $0.02$ | $0.006$
|  | SafeBayes | $0.92$ | $0.96$ | $0.02$ | $0.007$
|  | Holmes and Walker | $0.96$ | $0.96$ | $0.02$ | $0.006$
|  | Lyddon et al. | $1.12$ | $0.92$ | $0.02$ | $0.005$
3 | $100$ | GPC | $0.92$ | $0.97$ | $0.13$ | $0.044$
|  | SafeBayes | $0.87$ | $0.97$ | $0.13$ | $0.047$
|  | Holmes and Walker | $0.93$ | $0.96$ | $0.13$ | $0.041$
|  | Lyddon et al. | $1.69$ | $0.81$ | $0.13$ | $0.020$
| $200$ | GPC | $0.96$ | $0.97$ | $0.06$ | $0.018$
|  | SafeBayes | $0.90$ | $0.96$ | $0.06$ | $0.021$
|  | Holmes and Walker | $0.92$ | $0.96$ | $0.06$ | $0.020$
|  | Lyddon et al. | $1.37$ | $0.86$ | $0.06$ | $0.011$
| $400$ | GPC | $0.97$ | $0.96$ | $0.03$ | $0.009$
|  | SafeBayes | $0.91$ | $0.96$ | $0.03$ | $0.010$
|  | Holmes and Walker | $0.86$ | $0.96$ | $0.03$ | $0.012$
|  | Lyddon et al. | $1.26$ | $0.89$ | $0.03$ | $0.006$
Table 2: Comparison of average learning rate estimates ($\hat{\eta}$),
estimated coverage probabilities (Coverage), mean square error (MSE), and
total posterior variance (Variance) across different sample sizes and
misspecification degrees in the Non-normal Errors example.

###  5.4 Other experiments

Finally, we considered other types of misspecification in addition to those
presented above. These results are not presented here because all four
learning rate selection methods performed similarly and displaying a table of
similar numbers is not a good use of space. But it is worth mentioning in what
cases these methods perform comparably, and below is a brief summary of our
findings.

  * •

In cases where the heteroscedasticity is less extreme than in Section 5.2
above, in particular, with errors having non-constant variance but independent
of $x$, we found that all four learning rate selection methods performed well.
That is, the learning rate estimates were all similar and the credible regions
all had coverage probability near the nominal 95% level.

  * •

The example in Section 5.3 considered heavy-tailed error distributions. We
also considered cases where the error distribution was asymmetric, e.g., skew-
normal (Pérez-Rodríguez et al., 2018). Apparently, misspecification in the
shape of the error distribution has little effect because, as above, the
learning rate selection methods all performed well in these cases.

##  6 Learning rates in logistic regression

###  6.1 Model setup

An important problem in medical statistics is estimation of the so-called
minimum clinically important difference (MCID) that assesses the practical as
opposed to statistical significance of a treatment. In words, the MCID is the
threshold on the diagnostic measure scale such that improvements beyond that
level are associated with patients feeling better after the treatment; see,
e.g., Hedayat et al., (2015) and the references therein. To set the scene, let
$X\in\mathbb{R}$ denote the patient’s diagnostic measure, e.g., the pre-
treatment minus post-treatment difference in blood pressure, and let
$Y\in\\{-1,+1\\}$ denote the patient-reported indicator of whether they felt
the treatment was effective, with “$y=+1$” indicating effective. The quantity
of interest, $\theta$, the MCID, is the cutoff on the $X$ scale such that the
indicator $1\\{X>\theta\\}$ is most highly associated with $Y$. More
precisely, the MCID is defined as

$\theta=\arg\min_{\vartheta}P\\{Y\neq\text{sign}(X-\vartheta)\\},$

where $\text{sign}(0)=1$. Clearly, $\theta$ depends on the unknown joint
distribution $P$ of $(X,Y)$.

Towards inference on the MCID, it is natural to introduce a statistical model
for $P$. It would be difficult to develop a model for which $\theta$ is
directly a model parameter, but one idea would be to use a logistic regression
model with $Y$ as the binary response and $X$ as a continuous predictor. That
is, the logistic regression model states that

$(Y\mid X=x)\sim{\sf Rad}\bigl{(}F(\beta_{0}+\beta_{1}x)\bigr{)},$

where ${\sf Rad}(p)$ denotes a Rademacher distribution, i.e., a binary
distribution on $\\{-1,+1\\}$, with probability mass $p$ assigned to the value
$+1$, and $F$ is a logistic distribution function with
$F(u)=\\{1+e^{-u}\\}^{-1}$, for $u\in\mathbb{R}$. The logistic regression
model is determined by the unknown parameters $(\beta_{0},\beta_{1})$. Since
the MCID $\theta$ also satisfies $P(Y=+1\mid X=\theta)=\tfrac{1}{2}$, if the
above model is assumed, then

$\theta=-\beta_{0}/\beta_{1}.$

Given independent observations $(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})$ from this
model, a posterior distribution for $(\beta_{0},\beta_{1})$, generalized Bayes
or otherwise, can be obtained. From this, one can readily obtain the
corresponding posterior distribution of $\theta$ via the identity above.

Of course, this model could easily be misspecified. So it is of interest to
investigate what happens with the the generalized Bayes posterior with
suitably chosen learning rates when the logistic link function $F$ is
incorrectly specified.

###  6.2 Results

We fit a misspecified logistic regression model, i.e., where the diagnostic
measure $X$ comes from a normal mixture model with distribution function

$F^{\star}(x)=0.7\,\Phi(x\mid 5,1)+0.3\,\Phi(x\mid\mu,1),$

and the patient reported effectiveness indicator is $(Y\mid X=x)\sim{\sf
Rad}(F^{\star}(x))$. Here the quantity $\mu$ controls the degree of
misspecification, with $\mu$ closer to 5 corresponding to “less
misspecification”—relative to the logistic link function $F$ above—compared to
$\mu$ further from 5. The three specific degrees considered are:

Degree 1.

$\mu=7$;

Degree 2.

$\mu=8$;

Degree 3.

$\mu=9$.

For the posited logistic regression model, we follow Robert and Casella,
(2004, Example 7.11) and take a default prior distribution for
$(\beta_{0},\beta_{1})$ to be

$\pi(\beta_{0},\beta_{1})=\hat{b}^{-1}\exp\\{\beta_{0}-\hat{b}^{-1}e^{\beta_{0}}\\},$

which is simply a flat prior for $\beta_{1}$ and an exponential prior for
$e^{\beta_{0}}$ with scale $\hat{b}=\exp(\hat{\beta}_{0}+\gamma)$, where
$\hat{\beta}_{0}$ is the maximum likelihood estimator and $\gamma\approx
0.5772$ is Euler’s constant. This default prior is used here, rather than,
say, the Pólya–gamma prior of Polson et al., (2013), because embedding its
more sophisticated posterior sampling scheme inside the learning rate
selection procedures was simply too expensive.

The goal is, as in the previous section, to investigate the extent to which
the learning rate selection methods can help the generalized Bayes posterior
distribution to overcome the model misspecification, and Table 3 summarizes
the results. There we present the average learning rate value, the coverage
probability of 95% credible intervals for $\theta$, the average length of
those credible intervals, and mean square error, all based on 500
replications, for each pair of $\mu$ and sample size $n$. Here we see that, in
the Degree 1 case where misspecification is relatively mild, the methods
perform reasonably well in terms of coverage probability, but things get worse
as sample size increases, a symptom of the model misspecification bias. For
the Degree 2–3 cases with even more model misspecification, all the methods
perform quite poorly. Apparently none of the learning rate selection methods
can help the posterior overcome the relatively severe model misspecification
bias in this example.

Degree | $n$ | Method | $\hat{\eta}$ | Coverage | Length | MSE
---|---|---|---|---|---|---
1 | $100$ | GPC | $0.904$ | $0.938$ | $0.883$ | $0.051$
|  | SafeBayes | $0.790$ | $0.953$ | $0.987$ | $0.053$
|  | Holmes and Walker | $0.999$ | $0.927$ | $0.834$ | $0.051$
|  | Lyddon et al. | $1.003$ | $0.923$ | $0.830$ | $0.051$
| $200$ | GPC | $0.977$ | $0.914$ | $0.575$ | $0.029$
|  | SafeBayes | $0.913$ | $0.916$ | $0.599$ | $0.029$
|  | Holmes and Walker | $0.999$ | $0.902$ | $0.568$ | $0.030$
|  | Lyddon et al. | $1.003$ | $0.923$ | $0.830$ | $0.030$
| $400$ | GPC | $0.910$ | $0.912$ | $0.418$ | $0.015$
|  | SafeBayes | $0.822$ | $0.926$ | $0.450$ | $0.015$
|  | Holmes and Walker | $0.999$ | $0.890$ | $0.397$ | $0.015$
|  | Lyddon et al. | $0.988$ | $0.892$ | $0.401$ | $0.015$
2 | $100$ | GPC | $0.786$ | $0.866$ | $1.148$ | $0.139$
|  | SafeBayes | $0.890$ | $0.893$ | $1.283$ | $0.138$
|  | Holmes and Walker | $1.000$ | $0.836$ | $1.071$ | $0.137$
|  | Lyddon et al. | $0.986$ | $0.838$ | $1.085$ | $0.137$
| $200$ | GPC | $0.970$ | $0.788$ | $0.752$ | $0.083$
|  | SafeBayes | $0.906$ | $0.816$ | $0.786$ | $0.082$
|  | Holmes and Walker | $1.001$ | $0.776$ | $0.741$ | $0.082$
|  | Lyddon et al. | $0.974$ | $0.727$ | $0.742$ | $0.087$
| $400$ | GPC | $0.901$ | $0.622$ | $0.539$ | $0.067$
|  | SafeBayes | $0.833$ | $0.646$ | $0.574$ | $0.067$
|  | Holmes and Walker | $0.999$ | $0.564$ | $0.511$ | $0.067$
|  | Lyddon et al. | $0.969$ | $0.588$ | $0.518$ | $0.067$
3 | $100$ | GPC | $0.955$ | $0.742$ | $1.412$ | $0.345$
|  | SafeBayes | $0.891$ | $0.750$ | $1.482$ | $0.346$
|  | Holmes and Walker | $1.001$ | $0.726$ | $1.377$ | $0.344$
|  | Lyddon et al. | $0.970$ | $0.731$ | $1.385$ | $0.343$
| $200$ | GPC | $0.951$ | $0.532$ | $0.976$ | $0.266$
|  | SafeBayes | $0.887$ | $0.562$ | $1.024$ | $0.264$
|  | Holmes and Walker | $1.004$ | $0.502$ | $0.946$ | $0.265$
|  | Lyddon et al. | $0.955$ | $0.510$ | $0.963$ | $0.261$
| $400$ | GPC | $0.820$ | $0.244$ | $0.693$ | $0.247$
|  | SafeBayes | $0.893$ | $0.290$ | $0.743$ | $0.248$
|  | Holmes and Walker | $1.000$ | $0.228$ | $0.657$ | $0.248$
|  | Lyddon et al. | $0.951$ | $0.236$ | $0.671$ | $0.247$
Table 3: Summary of learning rate selection method performance in the
misspecified logistic regression example based on 500 replications.

###  6.3 A Gibbs posterior

The generalized Bayes posterior is not able to overcome this apparently rather
severe form of misspecification bias. As an alternative, we can consider a
different type of posterior construction, the so-called Gibbs posterior. What
distinguishes a Gibbs from a generalized Bayes posterior is that the former is
defined by a loss function while the latter is defined by a likelihood. In the
present context, the Gibbs posterior is a more appropriate approach since the
MCID is not naturally defined as a model parameter in a likelihood. Therefore,
it is possible to directly define the Gibbs posterior for the MCID, $\theta$,
as opposed to indirectly through a likelihood for $\beta$ and marginalizing to
$\theta$.

Define the loss function
$\ell_{\theta}(x,y)=\frac{1}{2}\\{1-y\,\text{sign}(x-\theta)\\}$ and the
corresponding risk (expected loss) $R(\vartheta)=P\ell_{\vartheta}$. As
Hedayat et al., (2015), showed, the MCID is the minimizer of $R$, i.e.,
$\theta^{\star}=\arg\min_{\theta}R(\theta)$. So the goal is to construct an
empirical version of the risk function, and then a sort of posterior
distribution that will concentrate around values that make the empirical risk
small. For the empirical risk, let

$R_{n}(\theta)=\frac{1}{n}\sum_{i=1}^{n}\ell_{\theta}(X_{i},Y_{i}).$

Then the Gibbs posterior distribution for $\theta$ has a density function
defined as

$\pi_{n}^{(\eta)}(\theta)\propto e^{-\eta nR_{n}(\theta)}\,\pi(\theta),$

where $\eta>0$ is, as before, the learning rate. In principle, all the
different learning rate selection methods considered above can be applied to
the Gibbs posterior framework to choose an appropriate value of $\eta$. Here,
however, the loss function is not differentiable, which creates a challenge
for the methods of Holmes and Walker, (2017) and Lyddon et al., (2019).
Therefore, in what follows, we only compare GPC and SafeBayes.

The results in Table 4 are to illustrate the performance of finding the
learning rate using Gibbs posterior with scaling algorithm in Syring and
Martin, (2019) and Grünwald, (2018). Since the true MCID is almost certain to
be in the range of observed $X$ values, the results here for both GPC and
SafeBayes are based on uniform prior on $[X_{(1)},X_{(n)}]$, the sample range.
Here we observe that the GPC is able to choose the learning rate such that the
desired 95% coverage target for each sample size. SafeBayes, on the other
hand, tends to choose too large of a learning rate, leading to (sometimes
severe) under-coverage.

Degree | $n$ | Method | $\hat{\eta}$ | Coverage | Length | MSE
---|---|---|---|---|---|---
1 | $100$ | GPC | $0.496$ | $0.948$ | $1.579$ | $0.103$
|  | SafeBayes | $0.982$ | $0.810$ | $0.910$ | $0.111$
| $200$ | GPC | $0.396$ | $0.932$ | $1.177$ | $0.060$
|  | SafeBayes | $0.986$ | $0.700$ | $0.587$ | $0.073$
| $400$ | GPC | $0.292$ | $0.952$ | $0.971$ | $0.033$
|  | SafeBayes | $0.975$ | $0.588$ | $0.378$ | $0.048$
2 | $100$ | GPC | $0.444$ | $0.972$ | $2.251$ | $0.228$
|  | SafeBayes | $0.966$ | $0.830$ | $1.216$ | $0.201$
| $200$ | GPC | $0.339$ | $0.950$ | $1.797$ | $0.143$
|  | SafeBayes | $0.967$ | $0.700$ | $0.750$ | $0.123$
| $400$ | GPC | $0.246$ | $0.970$ | $1.396$ | $0.073$
|  | SafeBayes | $0.964$ | $0.592$ | $0.490$ | $0.066$
3 | $100$ | GPC | $0.408$ | $0.966$ | $3.076$ | $0.548$
|  | SafeBayes | $0.953$ | $0.804$ | $1.536$ | $0.372$
| $200$ | GPC | $0.313$ | $0.958$ | $2.452$ | $0.351$
|  | SafeBayes | $0.964$ | $0.692$ | $0.892$ | $0.205$
| $400$ | GPC | $0.231$ | $0.964$ | $1.953$ | $0.187$
|  | SafeBayes | $0.965$ | $0.618$ | $0.544$ | $0.080$
Table 4: Summary of GPC and SafeBayes learning rate selection method
performance using a Gibbs posterior, based on 500 replications.

##  7 Conclusion

In this paper, we investigated the performance of several recently proposed
procedures for choosing the learning rate parameter in generalized Bayes
models. The goal was to see which, if any of these methods, are able to
overcome the model misspecification bias and give valid posterior uncertainty
quantification. While there are some models that are too severely misspecified
for a simple learning rate adjustment to accommodate, but we did find that
this can, in fact, be successful when misspecification is only relatively
mild.

A take-away message is that, among the learning rate selection methods
considered here, the GPC algorithm of Syring and Martin, (2019) seems to be
best suited overall for calibrating the generalized Bayes credible regions.
This is not surprising, given that is precisely what the GPC algorithm is
designed to do. GPC is computationally more expensive than, say, the method of
Lyddon et al., (2019), but our results here suggest that the extra time/effort
is well spent. Although GPC has been shown to have very good empirical
performance here and in a number of other references, there is still no formal
proof that it does indeed provide valid posterior uncertainty quantification.

Finally, our focus here was exclusively on inference, but it would be of
interest to see if/how different learning rate selection methods might assist
in generalized Bayes prediction. After all, the prediction problem is one
where it is possible to perform well even without a model, so developing a
learning rate selection method that would correct for certain kinds of model
misspecification, e.g., misspecified tails, should be within reach. That is,
can a suitable choice of learning rate ensure that quantiles of the posterior
predictive distribution achieve the nominal prediction coverage probability?

## Acknowledgments

This work is partially supported by the U.S. National Science Foundation,
DMS–1811802.

## References

  * 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.
  * Berk, (1966) Berk, R. H. (1966).  Limiting behavior of posterior distributions when the model is incorrect.  Ann. Math. Statist., 37(1):51–58.
  * Bhattacharya et al., (2019) Bhattacharya, A., Pati, D., and Yang, Y. (2019).  Bayesian fractional posteriors.  Ann. Statist., 47(1):39–66.
  * 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. Series B Stat. Methodol., 78(5):1103–1130.
  * Bunke and Milhaud, (1998) Bunke, O. and Milhaud, X. (1998).  Asymptotic behavior of Bayes estimates under possibly incorrect models.  Ann. Statist., 26(2):617–644.
  * De Blasi and Walker, (2013) De Blasi, P. and Walker, S. G. (2013).  Bayesian asymptotics with misspecified models.  Statist. Sinica, 23:169–187.
  * (8) Diaconis, P. and Freedman, D. (1986a).  On inconsistent Bayes estimates of location.  Ann. Statist., 14(1):68–87.
  * (9) Diaconis, P. and Freedman, D. (1986b).  On the consistency of Bayes estimates.  Ann. Statist., 14(1):1–26.
  * Fraser, (2011) Fraser, D. A. S. (2011).  Is Bayes posterior just quick and dirty confidence?  Statist. Sci., 26(3):299–316.
  * Ghosal et al., (1999) Ghosal, S., Ghosh, J. K., and Ramamoorthi, R. V. (1999).  Posterior consistency of Dirichlet mixtures in density estimation.  Ann. Statist., 27(1):143–158.
  * Ghosal et al., (2000) Ghosal, S., Ghosh, J. K., and van der Vaart, A. W. (2000).  Convergence rates of posterior distributions.  Ann. Statist., 28(2):500–531.
  * Ghosal and van der Vaart, (2017) Ghosal, S. and van der Vaart, A. (2017).  Fundamentals of Nonparametric Bayesian Inference, volume 44 of Cambridge Series in Statistical and Probabilistic Mathematics.  Cambridge University Press, Cambridge.
  * Grünwald, (2012) Grünwald, P. (2012).  The safe Bayesian: learning the learning rate via the mixability gap.  In Algorithmic Learning Theory, volume 7568 of Lecture Notes in Comput. Sci., pages 169–183. Springer, Heidelberg.
  * Grünwald, (2018) Grünwald, P. (2018).  Safe probability.  J. Statist. Plann. Inference, 195:47–63.
  * Grünwald and van Ommen, (2017) Grünwald, P. and van Ommen, T. (2017).  Inconsistency of Bayesian inference for misspecified linear models, and a proposal for repairing it.  Bayesian Anal., 12(4):1069–1103.
  * Grünwald and Mehta, (2020) Grünwald, P. D. and Mehta, N. A. (2020).  Fast rates for general unbounded loss functions: from ERM to generalized Bayes.  J. Mach. Learn. Res., 21(56):1–80.
  * Hedayat et al., (2015) Hedayat, A., Wang, J., and Xu, T. (2015).  Minimum clinically important difference in medical studies.  Biometrics, 71(1):33–41.
  * Heide et al., (2020) Heide, R., Kirichenko, A., Grunwald, P., and Mehta, N. (2020).  Safe-Bayesian generalized linear regression.  In International Conference on Artificial Intelligence and Statistics, pages 2623–2633. PMLR.
  * Holmes and Walker, (2017) Holmes, C. C. and Walker, S. G. (2017).  Assigning a value to a power likelihood in a general Bayesian model.  Biometrika, 104(2):497–503.
  * Huber, (1967) Huber, P. J. (1967).  The behavior of maximum likelihood estimates under nonstandard conditions.  In Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, volume 1, pages 221–233. University of California Press.
  * Kleijn and van der Vaart, (2006) Kleijn, B. J. K. and van der Vaart, A. (2006).  Misspecification in infinite-dimensional bayesian statistics.  Ann. Statist., 34(2):837–877.
  * Kleijn and van der Vaart, (2012) Kleijn, B. J. K. and van der Vaart, A. W. (2012).  The Bernstein-von-Mises theorem under misspecification.  Electronic Journal of Statistics, 6:354–381.
  * Lyddon et al., (2019) Lyddon, S. P., Holmes, C. C., and Walker, S. G. (2019).  General Bayesian updating and the loss-likelihood bootstrap.  Biometrika, 106(2):465–478.
  * Martin, (2017) Martin, R. (2017).  Invited comment on the article by van der Pas, Szabó, and van der Vaart.  Bayesian Anal., 12(4):1254–1258.
  * Martin, (2019) Martin, R. (2019).  False confidence, non-additive beliefs, and valid statistical inference.  Internat. J. Approx. Reason., 113:39–73.
  * Martin et al., (2017) Martin, R., Mess, R., and Walker, S. G. (2017).  Empirical Bayes posterior concentration in sparse high-dimensional linear models.  Bernoulli, 23(3):1822–1847.
  * Martin and Ning, (2020) Martin, R. and Ning, B. (2020).  Empirical priors and coverage of posterior credible sets in a sparse normal mean model.  Sankhyā A., 82:477–498.  Special issue in memory of Jayanta K. Ghosh.
  * Martin and Tang, (2020) Martin, R. and Tang, Y. (2020).  Empirical priors for prediction in sparse high-dimensional linear regression.  J. Mach. Learn. Res., 21(144):1–30.
  * Martin and Walker, (2019) Martin, R. and Walker, S. G. (2019).  Data-dependent priors and their posterior concentration rates.  Electron. J. Stat., 13(2):3049–3081.
  * Miller and Dunson, (2019) Miller, J. W. and Dunson, D. B. (2019).  Robust Bayesian inference via coarsening.  J. Amer. Statist. Assoc., 114(527):1113–1125.
  * Newton and Raftery, (1994) Newton, M. A. and Raftery, A. E. (1994).  Approximate Bayesian inference with the weighted likelihood bootstrap.  J. R. Stat. Soc. Series B Stat. Methodol., 56(1):3–26.
  * Pérez-Rodríguez et al., (2018) Pérez-Rodríguez, P., Acosta-Pech, R., Pérez-Elizalde, S., Cruz, C. V., Espinosa, J. S., and Crossa, J. (2018).  A Bayesian genomic regression model with skew normal random errors.  G3: Genes, Genomes, Genetics, 8(5):1771–1785.
  * Polson et al., (2013) Polson, N. G., Scott, J. G., and Windle, J. (2013).  Bayesian inference for logistic models using Pólya-Gamma latent variables.  J. Amer. Statist. Assoc., 108(504):1339–1349.
  * Ramamoorthi et al., (2015) Ramamoorthi, R. V., Sriram, K., and Martin, R. (2015).  On posterior concentration in misspecified models.  Bayesian Anal., 10:759–789.
  * Robert and Casella, (2004) Robert, C. and Casella, G. (2004).  Monte Carlo Statistical Methods.  Springer, New York, 2nd edition.
  * Syring and Martin, (2017) Syring, N. and Martin, R. (2017).  Gibbs posterior inference on the minimum clinically important difference.  J. Statist. Plann. Inference, 187:67–77.
  * Syring and Martin, (2019) Syring, N. and Martin, R. (2019).  Calibrating general posterior credible regions.  Biometrika, 106(2):479–486.
  * (39) Syring, N. and Martin, R. (2020a).  Gibbs posterior concentration rates under sub-exponential type losses.  arXiv:2012.04505.
  * (40) Syring, N. and Martin, R. (2020b).  Robust and rate-optimal Gibbs posterior inference on the boundary of a noisy image.  Ann. Statist., 48(3):1498–1513.
  * van der Vaart, (2000) van der Vaart, A. W. (2000).  Asymptotic Statistics, volume 3.  Cambridge university press.
  * Walker and Hjort, (2001) Walker, S. and Hjort, N. (2001).  On Bayesian consistency.  J. R. Stat. Soc. Series B Stat. Methodol., 63(4):811–821.
  * Walker, (2013) Walker, S. G. (2013).  Bayesian inference with misspecified models.  J. Statist. Plann. Inference, 143(10):1621–1633.
  * Walker et al., (2005) Walker, S. G., Lijoi, A., and Prünster, I. (2005).  Data tracking and the understanding of Bayesian consistency.  Biometrika, 92(4):765–778.
  * Wang and Martin, (2020) Wang, Z. and Martin, R. (2020).  Model-free posterior inference on the area under the receiver operating characteristic curve.  J. Statist. Plann. Inference, 209:174–186.
  * Zhang, (2006) Zhang, T. (2006).  From $\epsilon$-entropy to KL-entropy: analysis of minimum information complexity density estimation.  Ann. Statistx, 34(5):2180–2210.