,

# Breakdown of Boltzmann-type Models for the Alignment of Self-propelled Rods

Patrick Murphy111These authors contributed equally to this project
222corresponding author Department of Bioengineering, Rice University,
Houston, TX 77005  Center for Theoretical Biological Physics, Rice University,
Houston, TX 77005  Misha Perepelitsa${}^{\textrm{a}}$ Department of
Mathematics, University of Houston, TX 77204<EMAIL_ADDRESS>Ilya
Timofeyev  Department of Mathematics, University of Houston, TX 77204  Matan
Lieber-Kotz  Department of Bioengineering, Rice University, Houston, TX 77005
Brandon Islas  Department of Computational and Applied Mathematics, Rice
University, Houston, TX 77005  Oleg A. Igoshin${}^{\textrm{b,}}$ Department of
Bioengineering, Rice University, Houston, TX 77005  Center for Theoretical
Biological Physics, Rice University, Houston, TX 77005  Department of
Chemistry, Rice University, Houston, TX 77005  Department of Biosciences, Rice
University, Houston, TX 77005<EMAIL_ADDRESS>
###### Abstract

Studies in the collective motility of organisms use a range of analytical
approaches to formulate continuous kinetic models of collective dynamics from
rules or equations describing agent interactions. However, the derivation of
these kinetic models often relies on Boltzmann’s hypothesis of ”molecular
chaos”, that correlations between individuals are short-lived. While this
assumption is often the simplest way to derive tractable models, it is often
not valid in practice due to the high levels of cooperation and self-
organization present in biological systems. In this work, we illustrated this
point by considering a general Boltzmann-type kinetic model for the alignment
of self-propelled rods where rod reorientation occurs upon binary collisions.
We examine the accuracy of the kinetic model by comparing numerical solutions
of the continuous equations to an agent-based model that implements the
underlying rules governing microscopic alignment. Even for the simplest case
considered, our comparison demonstrates that the kinetic model fails to
replicate the discrete dynamics due to the formation of rod clusters that
violate statistical independence. Additionally, we show that introducing noise
to limit cluster formation helps improve the agreement between the analytical
model and agent simulations but does not restore agreement completely. These
results highlight the need to both develop and disseminate improved moment-
closure methods for modeling biological and active matter systems.

self-propelled rods, kinetic models, Boltzmann Equation, cluster formation,
statistical independence, agent-based model

##  I Introduction

Self-propelled rods are a fascinating class of active matter seen across both
living and non-living systems [27, 35, 26, 36, 37, 2, 1]. Due to their shape,
such systems are intrinsically capable of collective behavior through
realignments of rods due to physical collisions [4, 3] or longer-ranged
hydrodynamic interactions when in a fluid [35, 5]. These interactions lead to
the emergence of macroscopic collective motion such as flocking, clustering,
phase changes, and vortexes. Two notable biological examples of collective
motion in self-propelled rods are the dynamics of rod-shaped gliding bacteria
such as the soil bacterium Myxococcus xanthus [22, 43, 39], and the behavior
of groups of cellular cytoskeletal rods (such as F-actin [36] or microtubules
[37]) driven by molecular motors deposited on the surface. Collisions between
M. xanthus cells often result in the head of the colliding cell reorienting
along the length of the cell that was struck (see Fig. 1) [2]. Reorientation
itself occurs on a characteristic time scale related to the ratio of the cell
length to cell velocity. Microtubule collisions on a 2D surface also exhibit
asymmetric realignments upon collision, with the additional possibility of the
colliding microtubule stalling until the struck microtubule has passed [37].

Figure 1: An example of asymmetric alignment during the collision of two M.
xanthus cells due to cell-substrate forces. Adapted from [2].

Due to the interest surrounding systems of self-propelled rods, researchers
have developed various methods of modeling them. A popular approach to
studying the emergent behavior of these systems is by deriving a probabilistic
description [6, 12, 13, 18, 3, 30, 20, 8, 17, 33, 28]. A number of heuristic
mean-field models based on the notion of the mean nematic director have been
considered in [31, 15, 13, 14, 11, 12]. Since rod reorientation takes place
only when rods collide, a more accurate approach consists in tracing rod
collisions as they naturally occur and updating the kinetic distribution
function accordingly.

A common starting point for such derivations is an application of Boltzmann
framework that originated in the kinetic theory of gases. The framework is
based on 4 key elements: i) the geometry (shape) of interacting particles; ii)
the type of interactions, that determines how particle states change in
collisions; iii) the proper asymptotic regime in the model parameters; and iv)
the assumption of “molecular chaos”, that is the absence of two-particle
correlations. Importantly, the validity of hypothesis iv) depends on the
geometry of interacting particles and the type of interactions.

This modeling approach for self-propelled particles aligning through
collisions has been taken up by several authors [6, 18, 17, 28, 40, 19, 38].
In Bertin et al (2006) [6], authors consider spherically shaped cells that
experience polar alignment with some amount of noise in the post-collisional
orientations. Their model is based on the zero-order approximation of the
interaction operator, which amounts to treating terms of order
$\frac{Nl^{2}}{L^{2}}$ as negligible, where $N$ represents the number of
particles, $l$ denotes the particle size (radius), and $L$ signifies a
macroscopic length scale. This approximation signifies an asymptotic regime
where variations in the kinetic function across distances comparable to the
particle size can be disregarded, reminiscent of the classical Boltzmann
equation. Under these conditions and hypothesis iv), the authors derived a
kinetic equation for the distribution of cells and a reduced hydrodynamic-type
model for the first few moments of the cell density distribution, assuming
that the macroscopic motion is slow.

Cells of more realistic shape (rod-shaped cells) were considered in Hittmeir
et al (2021) [18]. The authors built a kinetic model in which a pair of
colliding co-oriented cells (rods with angle difference
$|\theta_{1}-\theta_{2}|<\pi/2$) change their orientation to their average
alignment; otherwise ($|\theta_{1}-\theta_{2}|>\pi/2$) both colliding cells
reverse their direction of motion. This complex interaction consists of a
partial polar alignment and a partial reversal. Again, the kinetic model is
obtained using the zero-order approximation of the interaction operator and
the molecular chaos assumption. For special sets of cell orientations, the
authors identify the corresponding set of equilibrium distributions and derive
a hydrodynamic-type model in the limit $\frac{Nl}{L}\to\infty.$

Over the last decade, the validity of Boltzmann formalism has come into
question. Evidence has shown that in some biological systems, this modeling
approach is not enough to reproduce the observed system dynamics due to either
weak binary interactions when stronger, multiparticle interactions are needed
[38], or due to rapid cluster formation that violates the molecular chaos
hypothesis [17, 40, 19]. These latter studies suggest that biologically-
relevant phenomena reliant on cluster formation, such as the transitions from
unordered to collective motion [30, 32, 41], may not be accurately modeled
with Boltzmann-type equations. However, results from Thüroff et al (2013) [40]
suggest that a Boltzmann approach can be appropriate for systems that align
weakly or gradually over time.

In this paper, we present a refined kinetic model that offers a more precise
representation of cell interactions, taking into account both the nature of
these interactions and the specific asymptotic regime governing the
interaction operator. Furthermore, we examine the compatibility between these
assumptions and the underlying hypothesis of molecular chaos.

We assume that cells are rod-shaped and when a given cell strikes a second
cell, it instantly turns around its “head” (the tip of the rod in the
direction of motion) to match a second cell’s orientation. Thus, in our model,
the re-orientations are truly nematic and also asymmetric, in agreement with
empirical observations of myxobacteria alignment, Fig 1. For the asymptotic
parameter regime, our model goes beyond the zero-order approximation used in
the above-mentioned papers, by including the next-order correction. In fact,
if one takes the values of model parameters from typical experiments with
myxobacteria, one can see that the second-order correction term (in a non-
dimensional equation) can be as large as 0.1 and, thus, is not negligible. We
derive our main kinetic equation for the nematic alignment under these
assumptions (section II, equation (12)). In section III, we reduce this
general kinetic equation to a system of numerically tractable PDEs by assuming
only a finite number of orientations are present. The system turns out to be
conservative and of the hyperbolic type. We then use a Lax-Friedrichs
numerical algorithm, specifically designed to treat such systems, to obtain
the numerical solutions.

Additionally, in section III we perform a comparative study of the numerical
solutions of the PDEs with the agent-based simulations to determine the
validity of the molecular chaos assumption. Our main finding here is that,
unless an appreciable amount of noise is added to the model, the molecular
chaos assumption does not hold. We show that even in simple cases the agent
model exhibits cluster formation, with the rate of clustering increasing with
the number of cells for the asymptotic regime considered.

As previously discussed, these results do not come as a complete surprise;
instead, they validate the criticisms directed towards the Boltzmann framework
when applied to scenarios involving self-propelled, rod-shaped cells engaged
in asymmetric nematic alignment.

##  II Derivation of kinetic equation

Our goal is to derive a tractable kinetic description of a system of $N$ self-
propelled rods in of length $l$ moving at constant speed $v$ in 2 dimensions,
starting from a microscopic description of collisions. Let $f({\bf
x},\theta,t)$ be the probability a rod with orientation $\theta$ has its head
located at ${\bf x}\in\mathbb{R}^{2}$ at time $t$. We will next derive a
kinetic equation for $f$ by finding an expression for the change in $f$ over a
single time step of size $\Delta t$. To this end, we make several assumptions
about rod collisions. First, we assume that collisions are binary events, with
the striking rod reorienting to match the struck rod’s orientation modulo
$\pi$. This results in the nematic alignment of rods. Motivated by M. xanthus
dynamics, we further assume that reorientation occurs through the head of the
rod. Since we assume binary collisions, we will make use of the 2-particle
distribution function $f_{2}({\bf x},\theta,{\bf x}_{1},\theta_{1},t)$ to help
calculate collision events.

Denote by $\Delta f$ the change of $f$ along the trajectory:

$\Delta f({\bf x},\theta,t)=f({\bf x}+v{\bf e}(\theta)\hskip 2.0pt\Delta
t,\theta,t+\Delta t)-f({\bf x},\theta,t),$

with $N\Delta f$ representing the change in the number of rods with a given
location and orientation. In the absence of collisions, the change is zero. If
a collision occurred it could result in the gain or loss of rods of
orientation $\theta.$ The associated geometry of rods for gain and loss is
determined by fixing a spatial location ${\bf x}$ and considering two cases
using two (non-interacting) rods with orientations $\theta$ and $\theta_{1}$
and heads at ${\bf x}$. We then look at the sets of collisions in the next
$dt$ time that will result in either a gain of a rod with orientation $\theta$
(due to a collision with a $\theta$-rod) or an equivalent loss (due to the
$\theta$-rod colliding). This yields two regions $P_{1}({\bf
x},\theta,\theta_{1})$ and $P_{2}({\bf x},\theta,\theta_{1})$ where another
rod could be located to cause a collision (see Figure 2).

Figure 2: Geometry of the collision scheme considered in our model. $P_{1}$ is
the set of heads of rods with orientation $\theta_{1}$ that will hit the side
of the rod with orientation $\theta$ having its head at ${\bf x},$ at time
$t,$ when it moves for $\Delta t$ units of time. Likewise, a rod with
orientation $\theta$ at ${\bf x}$ at time $t$ will hit the side of a rod with
orientation $\theta_{1}$ having its head in $P_{2}$ at time $t.$ $P_{1}$ is
formed by vectors $l{\bf e}(\theta)$ and $v\hskip 2.0pt\Delta t({\bf
e}(\theta_{1})-{\bf e}(\theta)).$ $P_{2}$ is formed by vectors $l{\bf
e}(\theta_{1})$ and $v\hskip 2.0pt\Delta t({\bf e}(\theta_{1})-{\bf
e}(\theta)).$

We can now use this setup to derive the change $\Delta f.$

$\begin{split}\Delta f({\bf
x},\theta,t)=&-(N-1)\int_{-\pi}^{\pi}\int_{P_{2}({\bf
x},\theta,\theta_{1})}f_{2}({\bf x},\theta,{\bf x}_{1},\theta_{1},t)\,d{\bf
x}_{1}d\theta_{1}\\\ &+(N-1)\int_{\theta-\pi/2}^{\theta+\pi/2}\int_{P_{1}({\bf
x},\theta,\theta_{1})}f_{2}({\bf x},\theta_{1},{\bf x}_{1},\theta,t)\,d{\bf
x}_{1}d\theta_{1}\\\ &+(N-1)\int_{\theta-\pi/2}^{\theta+\pi/2}\int_{P_{1}({\bf
x},\theta+\pi,\theta_{1})}f_{2}({\bf x},\theta_{1},{\bf
x}_{1},\theta+\pi,t)\,d{\bf x}_{1}d\theta_{1}\\\ \end{split}$ (1)

To get a closed-form equation for $f$ we employ a commonly used assumption
that 2-particle distribution $f_{2}({\bf x},\theta,{\bf x}_{1},\theta_{1},t)$
can be written as the product of the marginal distributions $f({\bf
x},\theta,t)$ and $f({\bf x}_{1},\theta_{1},t)$, i.e.

$f_{2}({\bf x},\theta,{\bf x}_{1},\theta_{1},t){}={}f({\bf x},\theta,t)f({\bf
x}_{1},\theta_{1},t),$ (2)

for all pairs $({\bf x},\theta)$ and $({\bf x}_{1},\theta_{1}).$ This first-
order moment closure is known as statistical independence, or molecular chaos,
as it implies that the presence of one rod does not affect the probability of
finding the other at the given position and orientation. Next, we write the
integral terms using a Taylor expansion in ${\bf x}_{1}$

$f({\bf x}_{1},\theta_{1},t)=f({\bf x},\theta_{1},t)+\gradient_{\bf x}f({\bf
x},\theta_{1},t)\cdot({\bf x}_{1}-{\bf x})+O(|{\bf x}_{1}-{\bf x}|^{2}),$ (3)

so that only an integral over $\theta_{1}$ remains. Note that $|{\bf
x}_{1}-{\bf x}|=O(l)$ since two rods must be close to enable collision. A
geometric computation shows that the regions $P_{1}$ and $P_{2}$ where
collisions can take place have area $|P_{1}|=|P_{2}|=lv\hskip 2.0pt\Delta
t|\sin(\theta-\theta_{1})|.$ Substituting the Taylor expansion into the
expression for $\Delta_{1}f$ and using the areas of $P_{1}$, $P_{2}$ with the
integrals

$\int_{P_{2}({\bf x},\theta,\theta_{1})}({\bf x}_{1}-{\bf x})\,d{\bf
x}_{1}{}={}\frac{l^{2}v\hskip 2.0pt\Delta t}{2}|\sin(\theta-\theta_{1})|\,{\bf
e}(\theta_{1})+O(l^{3}v\hskip 2.0pt\Delta t),$ (4)

and

$\int_{P_{1}({\bf x},\theta,\theta_{1})}({\bf x}_{1}-{\bf x})\,d{\bf
x}_{1}{}={}\frac{l^{2}v\hskip 2.0pt\Delta t}{2}|\sin(\theta-\theta_{1})|\,{\bf
e}(\theta)+O(l^{3}v\hskip 2.0pt\Delta t),$ (5)

we obtain

$\begin{split}&\Delta f({\bf x},\theta,t)=(N-1)lv\hskip 2.0pt\Delta
t\int\limits_{\theta-\pi/2}^{\theta+\pi/2}|\sin(\theta-\theta_{1})|f({\bf
x},\theta+\pi,t)f({\bf x},\theta_{1},t)\,d\theta_{1}\\\ &{}-{}(N-1)lv\hskip
2.0pt\Delta
t\int\limits_{\theta-\pi/2}^{\theta+\pi/2}|\sin(\theta-\theta_{1})|f({\bf
x},\theta,t)f({\bf x},\theta_{1}+\pi,t)\,d\theta_{1}\\\
&{}+{}\frac{(N-1)l^{2}v\hskip 2.0pt\Delta
t}{2}\int\limits_{\theta-\pi/2}^{\theta+\pi/2}|\sin(\theta-\theta_{1})|f({\bf
x},\theta_{1},t)\left({\bf e}(\theta)\cdot\gradient_{\bf x}f({\bf
x},\theta,t)+{\bf e}(\theta+\pi)\cdot\gradient_{\bf x}f({\bf
x},\theta+\pi,t)\right)\,d\theta_{1}\\\ &{}-{}\frac{(N-1)l^{2}v\hskip
2.0pt\Delta
t}{2}\int\limits_{\theta-\pi/2}^{\theta+\pi/2}|\sin(\theta-\theta_{1})|f({\bf
x},\theta,t)\left({\bf e}(\theta_{1})\cdot\gradient_{\bf x}f({\bf
x},\theta_{1},t)+{\bf e}(\theta_{1}+\pi)\cdot\gradient_{\bf x}f({\bf
x},\theta_{1}+\pi,t)\right)\,d\theta_{1}\\\ &{}+{}O(Nl^{3}v\hskip 2.0pt\Delta
t\max|f|\max|\gradient_{\bf x}^{2}f|){}+{}O(Nl^{2}\hskip 2.0pt\Delta
t^{2}).\end{split}$ (6)

Dividing the equation above by $\Delta t$ and letting $\Delta t\to 0,$ we
obtain the kinetic equation

$\partial_{t}f+v{\bf e}(\theta)\cdot\gradient_{\bf
x}f=(N-1)lvQ_{0}+\frac{(N-1)l^{2}v}{2}Q_{1}+O\left(Nl^{3}v\max|f|\max|\gradient_{\bf
x}^{2}f|\right),$ (7)

where

$Q_{0}{}={}\int\limits_{\theta-\pi/2}^{\theta+\pi/2}|\sin(\theta-\theta_{1})|\left(f({\bf
x},\theta+\pi,t)f({\bf x},\theta_{1},t)-f({\bf x},\theta,t)f({\bf
x},\theta_{1}+\pi,t)\right)\,d\theta_{1},$ (8)

and

$\begin{split}Q_{1}{}&={}\int\limits_{\theta-\pi/2}^{\theta+\pi/2}|\sin(\theta-\theta_{1})|f({\bf
x},\theta_{1},t){\bf e}(\theta)\cdot\gradient_{\bf x}\left(f({\bf
x},\theta,t)-f({\bf x},\theta+\pi,t)\right)\,d\theta_{1}\\\
&-\int\limits_{\theta-\pi/2}^{\theta+\pi/2}|\sin(\theta-\theta_{1})|f({\bf
x},\theta,t){\bf e}(\theta_{1})\cdot\gradient_{\bf x}(f({\bf
x},\theta_{1},t)-f({\bf x},\theta_{1}+\pi))\,d\theta_{1}.\end{split}$ (9)

In deriving $Q_{1}$ term we used the identity ${\bf e}(\theta+\pi)=-{\bf
e}(\theta).$

We nondimensionalize equation (7) by setting $\tau{}={}L/v$ and rescaling the
variables using

$\hat{{\bf x}}={\bf x}/L,\quad\hat{t}=t/\tau,\quad\hat{f}(\hat{{\bf
x}},\theta,\hat{t})=L^{2}f(\hat{{\bf x}}L,\theta,\hat{t}\tau).$ (10)

Dropping hats, the resulting scaled kinetic equation is

$\partial_{t}f+{\bf e}(\theta)\cdot\gradient_{\bf
x}f=\frac{(N-1)l}{L}Q_{0}+\frac{(N-1)l^{2}}{2L^{2}}Q_{1}+O\left(\frac{Nl^{3}}{L^{3}}\right),$
(11)

where we assumed that the variations in the kinetic density $f$ are bounded,
that is $\max|f|\max|\gradient^{2}_{\bf x}f|\leq C.$ For a typical experiment
with $N=1000$ M. xanthus bacteria ($l=5\mu m$) on a domain of size
$L=10^{3}\mu m,$ parameter $\frac{Nl}{L}=5$ and $\frac{Nl^{2}}{L^{2}}=0.025,$
and $\frac{Nl^{3}}{L^{3}}=0.000125.$ Thus, when $C$ is of order 1, a
reasonable approximation is keeping only the first two terms on the right-hand
side. We set $\kappa=\frac{Nl^{2}}{2L^{2}}.$ The parameter $\kappa$ is related
to the mean free path, $d,$ and rod length $l$ by

$\kappa\leq\frac{l}{2d}.$

See (19) of the Appendix for details. The binary collision assumption fails if
$d$ is of order $l$ or smaller, meaning that $\kappa$ should be small. In the
numerical simulations in the subsequent sections, we will restrict
$\kappa<0.5.$ The resulting kinetic equation becomes

$\partial_{t}f+{\bf e}(\theta)\cdot\gradient_{\bf x}f=\frac{2\kappa
L}{l}Q_{0}(f)+\kappa Q_{1}(f).$ (12)

In the model of alignment that we consider in this paper, a rod after an
interaction assumes the orientation of the rod it collides with. The set of
orientations $\theta_{i},$ $\theta_{i}+\pi,$ $i=1,\ldots,k$, if present
initially, will be preserved by the dynamics of transport and collision. In
such situations, the kinetic density $f$ is determined by the set of $2k$
densities of the corresponding orientations, and equation (12) can be written
as a system of $2k$ partial differential equations for the orientation
densities. An example of such a system is considered in the following
sections.

##  III Numerical solutions and results from agent-based simulations for case
study

In the absence of correlations in the underlying system, kinetic equations
derived from Boltzmann’s hypothesis should accurately model the system’s
evolution in time. However, the presence of correlations can cause
discrepancies to appear. Here we will test the accuracy of our kinetic
equation by comparing the behavior of numerical solutions to agent-based
simulations of the underlying system. We will use a test case to demonstrate
that even for simple setups, the microscopic rules of collision lead to a
buildup of orientational correlations, causing the correlation-free model to
underestimate the effects of rod-rod interactions. While this test case does
not capture the complexity present in many physical systems, it is enough to
demonstrate the need to account for correlations between rods when modeling
active matter systems such as bacterial communities.

###  III.1 A model with two orientations

Consider equation (12) for distributions with $n$ fixed orientations $f({\bf
x},\theta,t)=\sum_{i=1}^{n}\rho_{i}({\bf x},t)\delta(\theta-\theta_{i})$,
${\bf x}=[x,y]^{T}$. In general, this will result in a system of $n$ non-
conservative PDEs in two spatial dimensions, making accurate numerical
solutions over long periods of time hard to obtain. To illustrate our point
about growth of correlations, we will focus on a reduced case with two
orientation angles $\theta_{1}$ and $\theta_{2}$ where
$|\theta_{2}-\theta_{1}|\leq\pi/2.$ Then there is no nematic alignment and
thus the term $Q_{0}(f)$ vanishes. Additionally, as we will show, the system
of two PDEs can be written in conservation form under certain assumptions
about the initial density distributions. Substituting $f({\bf
x},\theta,t)=\rho_{1}({\bf x},t)\delta(\theta-\theta_{1})+\rho_{2}({\bf
x},t)\delta(\theta-\theta_{2})$ into 12 yields a system for
${\bf\rho}=(\rho_{1},\rho_{2})$ of the form

$\partial_{t}{\bf\rho}+A({\bf\rho})\partial_{x}{\bf\rho}+B({\bf\rho})\partial_{y}{\bf\rho}{}={}0$
(13)

where the matrices $A$ and $B$ are given by

$A{}={}\left[\begin{array}[]{rr}\cos(\theta_{1})+\kappa\cos(\theta_{2})\rho_{2}&\quad-\kappa\cos(\theta_{1})\rho_{1}\\\
-\kappa\cos(\theta_{2})\rho_{2}&\quad\cos(\theta_{2})+\kappa\cos(\theta_{1})\rho_{1}\end{array}\right],$
$B{}={}\left[\begin{array}[]{rr}\sin(\theta_{1})+\kappa\sin(\theta_{2})\rho_{2}&\quad-\kappa\sin(\theta_{1})\rho_{1}\\\
-\kappa\sin(\theta_{2})\rho_{2}&\quad\sin(\theta_{2})+\kappa\sin(\theta_{1})\rho_{1}\end{array}\right].$

To simplify these kinetic equations to a form easily analyzed, we make two
assumptions: that all densities vary only in the spatial direction $x$ so the
problem is effectively 1D with $\rho({\bf x},t)=\rho(x,t)$, and that
$\theta_{1}=\frac{\pi}{4},$ $\theta_{2}=\frac{3\pi}{4}$ so the system can be
put into conservation form. In this setting, the system (13) is simply

$\begin{split}\partial_{t}\rho_{1}+\frac{\sqrt{2}}{2}\partial_{x}\left(\rho_{1}\left[1-\kappa\rho_{2}\right]\right){}&={}0\\\
\partial_{t}\rho_{2}-\frac{\sqrt{2}}{2}\partial_{x}\left(\rho_{2}\left[1-\kappa\rho_{1}\right]\right){}&={}0.\end{split}$
(14)

We restrict our analysis to the regime $0<\kappa<0.5$ to enforce the mean free
path of rods is greater than the rod length $l$ (see Appendix A). We will
consider a test case with two opposing waves of rods over a constant
background density. Each wave is a Gaussian moving over a uniform background
density of rods. The initial conditions for $\rho_{1}$ and $\rho_{2}$ are of
the form

$\rho_{i}(x,0){}={}a_{i}+\frac{b_{i}}{\sqrt{2\pi}\sigma_{i}}e^{-(x-c_{i})^{2}/(2\sigma_{i}^{2})}$
(15)

with the constraints
$\int_{0}^{1}\int_{0}^{1}\rho_{1}(x,0)+\rho_{2}(x,0)dxdy=1$ from conservation
and an enforced constraint $a_{i}/b_{i}=r_{i}$ for the ratio of background
cells to the cells forming the waves. The final parameters $c_{i}$ and
$\sigma_{i}$ then control the location and width of the waves. Additionally,
the maximum amplitude $a_{i}+b_{i}/(\sqrt{2\pi}\sigma_{i})$ set below the
bound established in Appendix B. As time evolves, the linear factor
$1-\kappa\rho_{i}$ inside the spatial derivative will decrease both the left-
and right-moving waves’ speeds when they interact. Finally, we set
$c_{1}=0.25$, $c_{2}=0.75$, $\sigma_{1}=\sigma_{2}=0.0625$, and $r_{i}=7$ so
the two waves’ peaks will move towards each other and meet at $x=0.5$. We
numerically solve the system (14) using the method described in Appendix C.

###  III.2 Agent-based model setup

For the agent-based model, we set the domain size to $L=400$, the velocity to
$v=\sqrt{2}$. Each rod $i=1,\dots,N$ followed the equations of motion

$\begin{split}\frac{d{\bf
x}_{i}}{dt}{}&={}v[\cos(\theta_{i}),\;\sin(\theta_{i})]^{T}\\\
\frac{d\theta_{i}}{dt}{}&={}(\theta_{k}-\theta_{i})\delta(t-T_{i,k}).\end{split}$
(16)

where $\delta(t)$ is the Delta distribution and $T_{i,k}$ denotes the times
when the $i$th rod collides with another rod $k$ and reorients. The initial
location of each rod was determined through random sampling of the initial
density profiles $\rho_{i}(x,0)$. Collisions between rods were determined by
implementing the schematic in Figure 2. If a rod would collide with multiple
rods in a single time step, it was assumed that it collided with the first rod
in its path to ensure collisions are binary. We restricted the mean-free path
of rods to be greater than the rod length $l$ based on the physical timescale
of reorientation $l/v$ (see Appendix A). Therefore, we took time steps of
$\Delta t=l/v$. This discretization also has the advantage of reducing
computational time.

The simulations were deterministic after initial rod positions are sampled, so
we combined 1000 simulations to smooth fluctuations and recover the mean
behavior. Since the solutions of system (14) vary spatially only in the
$x$-direction, we compared them to agent-based simulations by looking only at
the $x$-coordinates of simulated rods. We then constructed a 1D kernel density
estimate (KDE) using these values. The resulting density profiles were then
scaled to probability densities on $x\in[0,\,1]$ so they could be compared to
the density profiles obtained from numerically solving the kinetic equations.

###  III.3 Discrepancy between agent-based model and kinetic equations due to
cluster formation

Our agent-based simulations deviate from the numerical solutions of the
kinetic equations, with the latter exhibiting a far less pronounced change in
the density profile shape (Figure 3 A&B). The peaks of the density profiles
for both orientations slow down slightly as they approach each other, but the
slowdown and the increase in the peak density are much greater in the agent
framework. If this discrepancy is due to not being near the limit
$N\to\infty$, $l\to 0$ with $\kappa$ held constant, increasing the number of
simulated rods while keeping $\kappa\propto Nl^{2}$ fixed should theoretically
improve the agreement. However, we see that this is not the case (Figure 3C).
The simulations using higher numbers of rods showed greater discrepancies.
This suggests the disagreement is due to one or more flaws in the assumptions
used to derive the kinetic equations.

Figure 3: A) Numerical solution of the kinetic equations for $\kappa=0.1$. B)
Density profile for agent-based simulations for $\kappa=0.1$ with $N=2856$
rods. C) Density profile for agent-based simulations for $\kappa=0.1$ with
$N=45700$ rods. (Inset) The discrepancy at the final time between the kinetic
model density profile and the density profile from agent simulations for
various numbers of rods. The results indicate that the agent simulations show
much greater changes in the density profile and do not converge to the kinetic
model in the limit $N\to\infty$, $l\to 0$.

The key assumption used in deriving the Boltzmann-type equation was the
statistical independence of the joint probability distribution (2). This
allowed us to express the probability that a pair of rods would have a spatial
configuration leading to a collision in terms of the probability of each rod
individually occupying the corresponding spatial region. This assumption is
invalid if there are correlations between rods. Such correlations appear in
clusters of aligned rods since they have a high chance of having similar
orientations due to collisions. Our collision scheme results in both rods
possessing the same orientation, so cluster formation is possible and would be
a violation of our assumption of statistical independence.

To measure the clustering present in our agent simulations, we used a standard
algorithm based on Euclidean distance to group rods. The minimum distance
between the heads of rods was picked to be the length of a rod $l$. We then
measured the proportion of rods in clusters of size 4 or more (Figure 4A). As
expected, clusters grew over time in our agent simulations, with the rate of
growth increasing with $\kappa$ for fixed $N$. This increase in $\kappa$
corresponds to a greater rod length, increasing the chance of collisions
between rods.

###  III.4 Quantifying loss of statistical independence

Since cluster formation is linked to a loss of statistical independence [17,
40, 19], we next quantified the extent of this loss. Observing a rod with
orientation $\theta_{1}$ would decrease the probability of a nearby rod having
orientation $\theta_{2}$. This is typically captured using two-particle or
higher correlation functions. If the system is spatially homogeneous or can be
approximated as such, then there are several ingenious ways to approach
calculating these quantities [21, 24, 23]. This is commonly used to analyze
systems near the onset of polar or nematic order. In our case, we are far from
spatial homogeneity, so we use a different approach.

In our work, we directly calculated metrics from agent simulation data related
to the quantity
$f_{2}(x_{1},\theta_{1},x_{2},\theta_{2},t)-f(x_{1},\theta_{1},t)f(x_{2},\theta_{2},t)$
to see if a loss of independence is present. We first divided our domain
$\Omega=[0,\,L]\times[0,\,L]$ into $2^{2m}$ square subregions $\Omega_{ij}$,
$i,\,j=1,\dots,2^{m}$ with side length $L/2^{m}$. Since subregions should be
large enough to contain small clusters of rods, we set $m=5$ so that
$L/2^{m}=12.5$. We pooled all rods appearing in $\Omega_{ij}$ across 1000
simulations with different initial agent locations ${\bf x}_{i}(0)$, then we
calculated coarse-grained approximations of the joint and marginal
distributions conditioned on rods being in $\Omega_{ij}$

$\begin{split}f_{2}^{ij}(\theta_{1},\theta_{2},t)&=\frac{\mbox{\\# pairs of
rods with orientations $(\theta_{1},\theta_{2})$ in box}\;ij}{\mbox{\\# pairs
of rods in box}\;ij}\\\ f^{ij}(\theta_{1},t)&=\frac{\mbox{\\#rods with
orientation $\theta_{1}$ in box}\;ij}{\mbox{\\#rods in box}\;ij}\\\
f^{ij}(\theta_{2},t)&=\frac{\mbox{\\# rods with orientation $\theta_{2}$ in
box}\;ij}{\mbox{\\# rods in box}\;ij}\end{split}.$ (17)

To measure the loss of statistical independence based on orientations in
$\Omega_{ij}$ at time $t$, we chose Pearson’s correlation coefficient
calculated from the estimates of the conditional joint and marginal

$r_{\theta_{1}\theta_{2}}^{ij}(t)=\frac{\left(\sum_{m,n=1,2}\theta_{n}\theta_{m}f_{2}^{ij}(\theta_{n},\theta_{m},t)\right)-M^{2}}{\sum_{n=1,2}(\theta_{n}-M)^{2}f^{ij}(\theta_{n},t)}$
(18)

where $M$ the mean orientation calculated from $f^{ij}(\theta_{1},t)$ and
$f^{ij}(\theta_{1},t)$. Averaging over all subregions then gives us an average
correlation $r_{\theta_{1}\theta_{2}}$.

As shown in Figure 4B, the metric $r_{\theta_{1}\theta_{2}}(t)$ shows similar
trends over time to the cluster formation. The loss of statistical
independence increases over time, with higher values of $\kappa$ showing a
greater loss for fixed $N$.This reflects the fact that the chance of finding a
pair of rods close together with different orientations is lowered from what
is expected if statistical independence holds.

Figure 4: A) The percentage of rods in clusters with more than 4 rods over
time. This quantity grows over time as rods collide and co-orient. B) The loss
of statistical independence measured via a local orientation correlation
coefficient $r_{\theta_{1}\theta_{2}}(t)$. Here the correlations growing over
time indicate that it becomes rarer to see nearby rods having different
orientations, e.g. cluster formation. C) Snapshots of an agent simulation at
the end of a simulation. The snapshot is taken at the center of the domain
where the two waves interact.

Since cluster formation is the result of co-alignment between rods upon
collision, it is natural to see if the increase in cluster growth with higher
values of $\kappa$ is due to a greater number of collisions. Calculating the
mean number of collisions per rod by time $t$ shows that the number does
increase as $\kappa$ increases (Figure 5A). Furthermore, there is strong
evidence of data collapse when plotting both the proportion of rods in
clusters (Figure 5B) and the loss of statistical independence versus the
cumulative number of collisions per rod (Figure 5C). There is a slight
difference in the rescaled curve for $\kappa=0.1$ compared to the other two
curves. This is the result of some rods starting in clusters at $t=0$ at
higher densities. The reasons for the deviation in the rescaled
$r_{\theta_{1}\theta_{2}}$ cure for $\kappa=0.1$ is less clear but could be
the result of cluster-cluster interactions once most rods reside in such
clusters. The rescaling is based on rod-rod collisions, so cluster-cluster
interactions might not fully be accounted for.

Figure 5: A) The total number of collisions experienced on average per rod
over time for $\kappa=0.001$ (blue), $\kappa=0.01$ (orange), and $\kappa=0.1$
(purple). B) The fraction of rods in clusters with more than 4 rods versus the
mean number of collisions experienced per rod by time $t$. This is plotted for
$\kappa=0.001$ (blue), $\kappa=0.01$ (orange), and $\kappa=0.1$ (purple) C)
The orientational correlation coefficient $r_{\theta_{1}\theta_{2}}(t)$
between nearby pairs of rods versus the mean number of collisions experienced
per rod by time $t$. Both figures show some data collapse when compared to
Figure 3, especially for low mean rod collisions per rod.

###  III.5 Improving agreement between agent simulations and kinetic equations
by inhibiting cluster formation

While the results of our simulations have indicated close links between
cluster formation, loss of independence, and the discrepancy between the
equations and agent simulations, they have not shown a strict cause and
effect. To illustrate that clustering is the main cause of the discrepancy, we
introduce diffusion in $y$-direction into our agent simulations. The addition
of this noise will cause rods forming a cluster to slowly drift apart at a
rate dependent on the noise strength. This addition will not impact the test
case we considered for the kinetic model. Since initial conditions are
constant in the $y$-direction, the system given in 14 will remain unchanged by
the vertical diffusion.

We implemented the same sets of agent simulations as before ($N=2856$,
$\kappa=0.1$) with the addition of different levels of noise $\sigma$. Agent
rods step in the $y$-direction a random distance drawn from a normal
distribution $N(0,\sigma)$ every $\Delta t=l/v$. Since agent simulations are
run with unscaled variables, we nondimensionalize the noise by using the
scaling $\sigma\rightarrow\sigma/(l/\sqrt{2})$ in order to compare the
strength of the noise to the rod length $l$. The results of these simulations
show that the addition of noise in the $y$-direction improves agreement
between the agent simulations and the kinetic equations, with the former now
resembling the latter for a sufficient level of noise (Figure 6A). Increasing
the noise strength from zero reduces the measured discrepancy between the
kinetic and agent density profiles (Figure 6B) and decreases the proportion of
rods in clusters (Figure 6C). However, the discrepancy is reduced only up to a
point. Once the strength of the scaled noise exceeds roughly $1$
(corresponding to the standard deviation of the normal distribution equaling
the projection of a rod in the y-direction), the discrepancy increases
slightly before plateauing. Why this is the case is unclear; however, the
reduction in clustering slows down at around the same level of noise. It is
possible that there are some aspects of rod correlations that the added
diffusion does not affect. For example, our collision scheme results in both
rods having similar x-coordinates. The addition of vertical diffusion does not
change this either. Therefore, successive collisions result in more rods
sharing similar x-coordinates, and potentially forming vertical bands of rods,
regardless of vertical diffusion.

Figure 6: A) Density profile for agent-based simulations with $y$-directional
scaled noise of $1.06$ for $\kappa=0.1$ with $N=2856$ rods. A scaled noise
level of 1 corresponds to an agent rod’s vertical movement over time $\Delta
t=l/v$ being distributed normally as $y(t+dt)-y(t)\sim
N(l\sin(\theta_{i}),l\sin(\theta_{i}))=N(l/\sqrt{2},l/\sqrt{2})$ in between
collisions, or equivalently to a vertical diffusion coefficient of $l^{2}/4$.
B) Percentage of rods in clusters with more than 4 rods versus noise in
$y$-direction. Increasing the noise in the vertical direction decreases the
proportion of rods in larger clusters C) Discrepancy between the kinetic model
and the agent-based simulations for different levels of white noise in the
$y$-direction. The discrepancy initially drops sharply but then increases
slightly once the strength of the noise increases past $1$, the length of a
rod projected in the y-direction. Snapshots of agent simulations near noise
levels of $0$, $1$, and $2$ are shown on the right. The snapshots are taken at
the final time in the center of the domain ($[150,\,250]\times[150,\,250]$)
where the two waves interact.

##  IV Discussion

In this paper, we developed a kinetic model for the alignment of self-
propelled hard rods where collisions result in asymmetric alignment. Using
this model, we showed Boltzmann formalism severely underestimates the change
in rod density profiles when two opposing waves of rods interact. We
explicitly measured the loss of statistical independence that invalidates the
classic assumption of molecular chaos Boltzmann-type equations rely on. Such a
loss corresponds to the formation of rod clusters due to alignment from binary
rod collisions. Our results mirror those in other studies [17, 40], however,
we have built upon these works by showing that mechanisms that destroy or
inhibit cluster formation help restore agreement between the kinetic model and
agent-based implementations of the microscopic alignment rules. As this
discrepancy occurs even in the simple setup we consider, our work highlights
the need to extend current methodologies beyond Boltzmann-type kinetic
equations in order to accurately capture the properties of active matter in
biological systems.

Boltzmann’s hypothesis can be justified when the mean free path is large
compared to the range of local interactions, however, this is rarely satisfied
at realistic densities when ordered motion is established [34]. Additionally,
several studies have indicated that rapid cluster formation can lead to a
strong violation of the molecular chaos assumption needed in the Boltzmann
approach [17, 40]. While gradual alignment results in a better match, in
principle any collision rule between particles resulting in alignment can
cause correlations to appear. Such collisions are prevalent in collective
dynamics at the cellular level due to a low Reynolds number, where cells must
actively exert energy to maintain their motion. This is in contrast to
Boltzmann gas dynamics, where collisions are assumed to be non-elastic and
conserve momentum, resulting in particles simply bouncing off each other. In
addition to violations of Boltzmann’s hypothesis, the assumption of binary
interactions used in such models is not enough to reproduce the observed
dynamics in some biological systems [38]. Binary interactions can simply be
too weak to produce the alignment seen experimentally, even when corrections
are made to account phenomenologically for correlations. These breakdowns of
Boltzmann-type models suggest they are an overly-simplistic approach to
modeling the emergence or stability of collective alignment. The emergence of
local order almost by definition involves the breakdown of statistical
independence and the growth of correlations as agents align and start moving
together, violating Boltzmann’s hypothesis.

There are several ways to extend kinetic models to incorporate orientational
correlations. Perhaps the simplest is to explicitly extend the kinetic model
to explicitly include equations for clusters of various sizes. Gain and loss
terms for these additional equations are then the natural result of collisions
among clusters and rods [29, 42]. Peruani et al (2010) developed a framework
where cluster size was explicitly tracked in a hierarchy of equations, while
Weber et al (2013) [42] considered a simplified framework with two reaction
equations for clusters and single cells. Another set of approaches for
incorporating correlations involves direct modifications to the moment closure
method. A simple example is replacing the joint distribution
$f_{2}(x_{1},\theta_{1},x_{2},\theta_{2},t)$ with
$\chi(\theta_{2}-\theta_{1})f(x_{1},\theta_{1},t)f(x_{2},\theta_{2},t)$, where
$\chi(\theta_{2}-\theta_{1})$ is a phenomenological term accounting for
correlations between different angles [38]. More sophisticated methods involve
higher-order moment closures of the BBGKY hierarchy [7, 2, 16]) using the so-
called cluster expansion that explicitly incorporates the evolution of the
joint distribution $f_{2}$ [9]. Such an approach was used by Chou and Ihle for
Vicsek-style models to extend beyond mean-field theory [10]. With the rapid
advancement of biological studies in the last decades, developing and applying
new analytical models to understand active matter in biology is of crucial
importance. Creating a tractable class of models that can capture correlations
or non-binary interactions would provide a cornerstone for this young field.

##  V Acknowledgments

Research is supported by National Science Foundation Division of Mathematical
Sciences awards 1903270 (to M.P. and I.T.) and 1903275 (to OAI and PM). We
also would like to thank prof. D. Kuzmin (TU Dortmund) for helpful discussions
on the numerical methods for the kinetic model presented in this paper.

##  Appendix A Mean free path

The mean free path (MFP) is defined as the average distance that a rod moves
between collisions. This distance can be computed as the product of the speed
$v$ and the time $\tau$ a rod moves between re-orientations. To estimate
$\tau,$ we first compute the number of rods with orientation ${\bf
e}(\theta)=(\cos\theta,\sin\theta)$ that will collide with a given rod (rod 1)
in the time $[t^{*},t^{*}+\Delta t]$, were $t^{*}$ marks the last change in
rod 1’s orientation. Without loss of generality, we consider rod 1 as moving
horizontally to the right, as in Figure 7. Suppose that a coordinate system is
chosen so that the rod is at rest. Then the velocity of rods with orientation
$\theta$ is described by the vector ${\bf w}{}={}v(\cos\theta-1,\sin\theta).$
The heads of the rods that rod 1 can collide with are located in a
parallelogram $P$ formed by vectors ${\bf w}\hskip 2.0pt\Delta t$ and $l{\bf
e}(\theta).$ The area of this parallelogram equals $|\sin\theta|vl\hskip
2.0pt\Delta t.$ Using that $Nf(x,\theta,t^{*})$ is the local number of rods
with a given orientation, we estimate the number of collisions with
$\theta$-oriented rods in time $dt$ as $vlNf(x,\theta,t^{*})|\sin\theta|\hskip
2.0pt\Delta t,$ which is maximized at $\theta=\pi/2.$ Thus, the number of
collisions over all orientations in $\Delta t$ time can be estimated as

$vlN\hskip 2.0pt\Delta t\int f(x,\theta,t^{*})\,d\theta{}\leq{}n(x)vl\hskip
2.0pt\Delta t,$

where $n(x,t^{*})$ is the local rod number density. From this, we obtain an
upper bound for the frequency of collisions as $n(x,t^{*})vl$. The time
between collisions and the MFP are then bounded by

$\tau\geq\frac{1}{n(x,t^{*})vl},\quad d\geq\frac{1}{n(x,t^{*})l}.$

Note that both $\tau$ and the MFP $d$ are local quantities depending on $x.$

Figure 7: Vector geometry used in the estimation of the MFP. A) rod at $x$
where ${\bf e}(\theta)=(\cos\theta,\sin\theta)$ moves with velocity $v{\bf
e}(\theta).$ In a reference coordinates of a horizontally moving rod with
velocity $(v,0),$ the rod velocity is ${\bf w}.$ B) a rod moving horizontally
to the right, during time interval $[t^{*},t^{*}+\Delta t]$ can collide with a
rod which has orientation ${\bf e}(\theta)$ whose center is located in the
parallelogram P, formed by vectors ${\bf w}\hskip 2.0pt\Delta t$ and $l{\bf
e}(\theta).$ The area of the parallelogram equals $|\sin\theta|vl\hskip
2.0pt\Delta t.$

For the model we consider in this paper, we restrict the MFP to $d\geq l$. The
rationale behind this is physics-based. When two rods collide, there is a
characteristic time scale for reorientation to occur given by
$\tau_{\theta}=l/v$ [1]. Therefore it makes sense to restrict the mean time
between collisions $\tau$ to the regime $\tau\geq\tau_{\theta}=l/v$. This
naturally yields $d\geq l$.

Our assumption that $d\geq l$ puts a restriction on the local density of the
form $n(x,t^{*})l^{2}\leq 1,$ and subsequently on the parameter $\kappa$ of
the form

$\kappa=\frac{Nl^{2}}{2L^{2}}\leq\frac{1}{2}l^{2}\max_{x}\,n(x,t^{*})\leq\frac{1}{2}.$
(19)

Here we use that $N/L^{2}$ is a lower bound on the maximum of $n(x,t^{*}).$
Note that we cannot have $\kappa=0$ without the local density being
identically 0.

##  Appendix B Domain of hyperbolicity of system of equations(14)

Here we will determine conditions under which system (14) is hyperbolic. The
system is given by

$\begin{split}\partial_{t}\rho_{1}+\frac{\sqrt{2}}{2}\partial_{x}\left[\rho_{1}(1-\kappa\rho_{2})\right]{}&={}0\\\
\partial_{t}\rho_{2}-\frac{\sqrt{2}}{2}\partial_{x}\left[\rho_{2}(1-\kappa\rho_{1})\right]{}&={}0,\end{split}$
(20)

where $\kappa=\dfrac{Nl^{2}}{2L^{2}}$. Denote the column vector
$U=(\rho_{1},\rho_{2})^{T}$ and column of the fluxes
$F(U)=\frac{\sqrt{2}}{2}(\rho_{1}-\kappa\rho_{1}\rho_{2},-\rho_{2}+\kappa\rho_{1}\rho_{2})^{T}.$
Then the system of equations is expressed as

$\partial_{t}U+\gradient_{x}F(U){}={}0.$

The system is hyperbolic if the eigenvalues of the gradient matrix

$\gradient_{U}F(U){}={}\left[\begin{array}[]{cc}1-\kappa\rho_{2}&-\kappa\rho_{1}\\\
\kappa\rho_{2}&-1+\kappa\rho_{1}\end{array}\right]$

are real. The eigenvalues equal

$\lambda_{\pm}=\frac{1}{\sqrt{2}}\left(\kappa(\rho_{1}-\rho_{2})\pm\sqrt{\kappa^{2}(\rho_{1}-\rho_{2})^{2}+4-4\kappa(\rho_{1}+\rho_{2})}\right).$

Thus, the system is hyperbolic whenever

$\kappa^{2}(\rho_{1}-\rho_{2})^{2}+4-4\kappa(\rho_{1}+\rho_{2})\geq 0$

with a sufficient condition

$\kappa(\rho_{1}+\rho_{2})\leq 1.$

Since (19) implies $\kappa<0.5$, our kinetic model (14) must be restricted to
cases where the sum of the maximum density profiles for $\rho_{1}$ and
$\rho_{2}$ total less than 2.

##  Appendix C Numerical Methods

We express the system in (14) in flux form

$\partial_{t}U=-\partial_{x}F,$ (21)

with $U=(\rho_{1},\rho_{2})^{T}$ and $F=A(\rho_{1},\rho_{2})^{T}$. In
particular, for two orientations $\pi/4$ and $3\pi/4$ the fluxes become
$F_{1}=\rho_{1}(1+\rho_{2})$ and $F_{2}=\rho_{2}(1+\rho_{1})$. Next, we
discretize the equation in space and time using the Lax-Friedrichs method [25]

$U_{i}^{n+1}=U_{i}^{n}-\frac{\Delta t}{\Delta
x}[\hat{F}_{i+1/2}-\hat{F}_{i-1/2}]$

and define

$\hat{F}_{i+1/2}=\frac{F(U_{i+1}^{n})+F(U_{i}^{n})}{2}-\frac{K}{2}(U_{i+1}^{n}-U_{i}^{n})$

where $K=\cos(\theta_{1})(1+|\sin(\theta_{2}-\theta_{1})|)$ is the upper bound
on the speed of propagation given $\rho_{1},\,\rho_{2}<1$. Therefore, the CFL
condition for the scheme above is $K\Delta t/\Delta x<1$, and the additional
diffusion introduced by the scheme is proportional to $D=K\Delta x/2$. The
scheme is first order in space and time, but we found that it was sufficient
for our simulations.

## References

  * [1] Rajesh Balagam and Oleg A. Igoshin.  Mechanism for Collective Cell Alignment in Myxococcus xanthus Bacteria.  PLOS Computational Biology, 11(8):e1004474, August 2015.
  * [2] Rajesh Balagam, Douglas B. Litwin, Fabian Czerwinski, Mingzhai Sun, Heidi B. Kaplan, Joshua W. Shaevitz, and Oleg A. Igoshin.  Myxococcus xanthus Gliding Motors Are Elastically Coupled to the Substrate as Predicted by the Focal Adhesion Model of Gliding Motility.  PLOS Computational Biology, 10(5):e1003619, May 2014.
  * [3] Aparna Baskaran and M. Cristina Marchetti.  Enhanced Diffusion and Ordering of Self-Propelled Rods.  Phys. Rev. Lett., 101(26):268101, December 2008.
  * [4] Aparna Baskaran and M. Cristina Marchetti.  Hydrodynamics of self-propelled hard rods.  Phys. Rev. E, 77(1):011920, January 2008.
  * [5] Aparna Baskaran and M. Cristina Marchetti.  Statistical mechanics and hydrodynamics of bacterial suspensions.  PNAS, 106(37):15567–15572, September 2009.
  * [6] E. Bertin, M. Droz, and G. Grégoire.  Boltzmann and hydrodynamic description for self-propelled particles.  Phys Rev E Stat Nonlin Soft Matter Phys, 74, 2006.
  * [7] N. N. Bogoliubov.  Kinetic equations.  Journal of Physics USSR, 10:265–274, 1946.
  * [8] François Bolley, José A. Cañizo, and José A. Carrillo.  Mean-field limit for the stochastic Vicsek model.  Applied Mathematics Letters, 25(3):339–343, March 2012.
  * [9] C. Cercignani.  The Boltzmann Equation and Its Applications.  Springer-Verlag, New York, 1988.
  * [10] Yen-Liang Chou and Thomas Ihle.  Active matter beyond mean-field: Ring-kinetic theory for self-propelled particles.  Physical Review E, 91(2):022103, February 2015.  Publisher: American Physical Society.
  * [11] Pierre Degond, Angelika Manhart, and Hui Yu.  A continuum model for nematic alignment of self-propelled particles.  Discrete & Continuous Dynamical Systems - B, 22(4):1295, 2017.
  * [12] Pierre Degond, Angelika Manhart, and Hui Yu.  An age-structured continuum model for myxobacteria.  Math. Models Methods Appl. Sci., 28(09):1737–1770, August 2018\.
  * [13] Pierre Degond and Sébastien Motsch.  Continuum limit of self-driven particles with orientation interaction.  Math. Models Methods Appl. Sci., 18(supp01):1193–1215, August 2008\.
  * [14] Pierre Degond and Tong Yang.  Diffusion in a continuum model of self-propelled particles with alignment interaction.  Math. Models Methods Appl. Sci., 20(supp01):1459–1490, September 2010.
  * [15] F. Ginelli, F. Peruani, M. Bär, and H. Chaté.  Large-scale collective properties of self-propelled rods.  Phys Rev Lett, 104, 2010.
  * [16] Harold Grad.  On the kinetic theory of rarefied gases.  Communications on pure and applied mathematics, 2:331–407, 1949\.
  * [17] Timo Hanke, Christoph A. Weber, and Erwin Frey.  Understanding collective dynamics of soft active colloids by binary scattering.  Physical Review E, 88(5):052309, November 2013.  Publisher: American Physical Society.
  * [18] Sabine Hittmeir, Laura Kanzler, Angelika Manhart, and Christian Schmeiser.  Kinetic modelling of colonies of myxobacteria.  Kinetic and Related Models, 14(1):1–24, 2021.
  * [19] T. Ihle.  Towards a quantitative kinetic theory of polar active matter.  The European Physical Journal Special Topics, 223(7):1293–1314, June 2014.
  * [20] Thomas Ihle.  Kinetic theory of flocking: Derivation of hydrodynamic equations.  Physical Review E, 83(3):030901, March 2011.  Publisher: American Physical Society.
  * [21] Julian Jeggle, Joakim Stenhammar, and Raphael Wittkowski.  Pair-distribution function of active brownian spheres in two spatial dimensions: Simulation results and analytic representation.  The Journal of Chemical Physics, 152(19):194903, 2020.
  * [22] Dale Kaiser.  Coupling cell movement to multicellular development in myxobacteria.  Nat Rev Microbiol, 1(1):45–54, October 2003.
  * [23] Rüdiger Kürsten and Thomas Ihle.  Quantitative kinetic theory of flocking with three-particle closure.  Physical Review E, 104(3):034604, 2020.
  * [24] Rüdiger Kürsten, Sven Stroteich, Martín Zumaya Hernández, and Thomas Ihle.  Multiple Particle Correlation Analysis of Many-Particle Systems: Formalism and Application to Active Matter.  Physical Review Letters, 124(8):088002, 2020.
  * [25] R. LeVeque.  Finite volume methods for hyperbolic problems.  Cambridge University Press, 2002.
  * [26] Vijay Narayan, Sriram Ramaswamy, and Narayanan Menon.  Long-Lived Giant Number Fluctuations in a Swarming Granular Nematic.  Science, 317(5834):105–108, July 2007.
  * [27] F. J. Ndlec, T. Surrey, A. C. Maggs, and S. Leibler.  Self-organization of microtubules and motors.  Nature, 389(6648):305–308, September 1997.
  * [28] Misha Perepelitsa, Ilya Timofeyev, Patrick Murphy, and Oleg A. Igoshin.  Mean-field model for nematic alignment of self-propelled rods.  Physical Review E, 106(3):034613, September 2022.  Publisher: American Physical Society.
  * [29] F. Peruani, L. Schimansky-Geier, and M. Bär.  Cluster dynamics and cluster size distributions in systems of self-propelled particles.  Eur. Phys. J. Spec. Top., 191(1):173–185, December 2010.
  * [30] Fernando Peruani, Andreas Deutsch, and Markus Bär.  Nonequilibrium clustering of self-propelled rods.  Phys. Rev. E, 74(3):030904, September 2006.
  * [31] Fernando Peruani, Andreas Deutsch, and Markus Bär.  A mean-field theory for self-propelled particles interacting by velocity alignment mechanisms.  Eur. Phys. J. Spec. Top., 157:111–122, 2008.
  * [32] Fernando Peruani, Jörn Starruß, Vladimir Jakovljevic, Lotte Søgaard-Andersen, Andreas Deutsch, and Markus Bär.  Collective Motion and Nonequilibrium Cluster Formation in Colonies of Gliding Bacteria.  Physical Review Letters, 108(9):098102, February 2012.  Publisher: American Physical Society.
  * [33] A. Peshkov, E. Bertin, F. Ginelli, and H. Chaté.  Boltzmann-Ginzburg-Landau approach for continuous descriptions of generic Vicsek-like models.  The European Physical Journal Special Topics, 223(7):1315–1344, June 2014.
  * [34] Len Pismen.  Active Matter Within and Around Us: From Self-Propelled Particles to Flocks and Living Forms.  The Frontiers Collection. Springer International Publishing, Cham, 2021.
  * [35] Ingmar H. Riedel, Karsten Kruse, and Jonathon Howard.  A Self-Organized Vortex Array of Hydrodynamically Entrained Sperm Cells.  Science, 309(5732):300–303, July 2005.
  * [36] Volker Schaller, Christoph Weber, Christine Semmrich, Erwin Frey, and Andreas R. Bausch.  Polar patterns of driven filaments.  Nature, 467(7311):73–77, September 2010.
  * [37] Yutaka Sumino, Ken H. Nagai, Yuji Shitaka, Dan Tanaka, Kenichi Yoshikawa, Hugues Chaté, and Kazuhiro Oiwa.  Large-scale vortex lattice emerging from collectively moving microtubules.  Nature, 483(7390):448–452, March 2012.
  * [38] Ryo Suzuki, Christoph A. Weber, Erwin Frey, and Andreas R. Bausch.  Polar pattern formation in driven filament systems requires non-binary particle collisions.  Nature Physics, 11(10):839–843, October 2015.  Number: 10 Publisher: Nature Publishing Group.
  * [39] Shashi Thutupalli, Mingzhai Sun, Filiz Bunyak, Kannappan Palaniappan, and Joshua W. Shaevitz.  Directional reversals enable Myxococcus xanthus cells to produce collective one-dimensional streams during fruiting-body formation.  J R Soc Interface, 12(109):20150049, August 2015.
  * [40] Florian Thüroff, Christoph A. Weber, and Erwin Frey.  Critical Assessment of the Boltzmann Approach to Active Systems.  Physical Review Letters, 111(19):190601, November 2013.  Publisher: American Physical Society.
  * [41] Christoph A. Weber, Volker Schaller, Andreas R. Bausch, and Erwin Frey.  Nucleation-induced transition to collective motion in active systems.  Physical Review E, 86(3):030901, September 2012.  Publisher: American Physical Society.
  * [42] Christoph A. Weber, Florian Thüroff, and Erwin Frey.  Role of particle conservation in self-propelled particle systems.  New Journal of Physics, 15(4):045014, April 2013.  Publisher: IOP Publishing.
  * [43] Yilin Wu, Yi Jiang, Dale Kaiser, and Mark Alber.  Social Interactions in Myxobacterial Swarming.  PLoS Comput Biol, 3(12):e253, December 2007.


Marginal Likelihood Estimation via Arrogance Sampling

By Benedict Escoto

###### Abstract

This paper describes a method for estimating the marginal likelihood or Bayes
factors of Bayesian models using non-parametric importance sampling
(“arrogance sampling”). This method can also be used to compute the
normalizing constant of probability distributions. Because the required inputs
are samples from the distribution to be normalized and the scaled density at
those samples, this method may be a convenient replacement for the harmonic
mean estimator. The method has been implemented in the open source R package
margLikArrogance.

##  1 Introduction

When a Bayesian evaluates two competing models or theories, $T_{1}$ and
$T_{2}$, having observed a vector of observations $\boldsymbol{x}$, Bayes’
Theorem determines the posterior ratio of the models’ probabilities:

$\frac{p(T_{1}|\boldsymbol{x})}{p(T_{2}|\boldsymbol{x})}=\frac{p(\boldsymbol{x}|T_{1})}{p(\boldsymbol{x}|T_{2})}\frac{p(T_{1})}{p(T_{2})}.$
(1)

The quantity $\frac{p(\boldsymbol{x}|T_{1})}{p(\boldsymbol{x}|T_{2})}$ is
called a _Bayes factor_ and the quantities $p(\boldsymbol{x}|T_{1})$ and
$p(\boldsymbol{x}|T_{2})$ are called the theories’ _marginal likelihoods_.

The types of Bayesian models considered in this paper have a fixed finite
number of parameters, each with their own probability function. If
$\boldsymbol{\theta}$ are parameters for a model $T$, then

$p(\boldsymbol{x}|T)=\int
p(\boldsymbol{x}|\boldsymbol{\theta},T)p(\boldsymbol{\theta}|T)\,d\boldsymbol{\theta}=\int
p(\boldsymbol{x}\wedge\boldsymbol{\theta}|T)\,d\boldsymbol{\theta}$ (2)

Unfortunately, this integral is difficult to compute in practice. The purpose
of this paper is to describe one method for estimating it.

Evaluating integral (2) is sometimes called the problem of computing
normalizing constants. The following formula shows how $p(\boldsymbol{x}|T)$
is a normalizing constant.

$p(\boldsymbol{\theta}|\boldsymbol{x},T)=\frac{p(\boldsymbol{\theta}\wedge\boldsymbol{x}|T)}{p(\boldsymbol{x}|T)}$
(3)

Thus the marginal likelihood $p(\boldsymbol{x}|T)$ is also the normalizing
constant of the posterior parameter distribution
$p(\boldsymbol{\theta}|\boldsymbol{x},T)$ assuming we are given the density
$p(\boldsymbol{\theta}\wedge\boldsymbol{x}|T)$ which is often easy to compute
in Bayesian models. Furthermore, Bayesian statisticians typically produce
samples from the posterior parameter distribution
$p(\boldsymbol{\theta}|\boldsymbol{x},T)$ even when not concerned with theory
choice. In these case, computing the marginal likelihood is equivalent to
computing the normalizing constant of a distribution from which samples and
the scaled density at these samples are available. The method described in
this paper takes this approach.

##  2 Review of Literature

Given how basic (1) is, it is perhaps surprising that there is no easy and
definitive way of applying it, even for simple models. Furthermore, as the
dimensionality and complexity of probability distributions increase, the
difficulty of approximation also increases. The following three techniques for
computing bayes factors or marginal likelihoods are important but will not be
mentioned further here.

  1. 1.

Analytic asymptotic approximations such as Laplace’s method, see for instance
Kass and Raftery (1995),

  2. 2.

Bridge sampling/path sampling/thermodynamic integration (Gelman and Meng,
1998), and

  3. 3.

Chib’s MCMC approximation (Chib, 1995; Chib and Jeliazkov, 2005).

Kass and Raftery (1995) is a popular overview of the earlier literature on
Bayes factor computation. All these methods can be very successful in the
right circumstances, and can often handle problems too complex for the method
described here. However, the method of this paper may still be useful due to
its convenience.

The rest of section 2 describes three approaches that are relevant to this
paper.

###  2.1 Importance Sampling

Importance sampling is a technique for reducing the variance of monte carlo
integration. This section will note some general facts; see Owen and Zhou
(1998) for more information.

Suppose we are trying to compute the (possibly multidimensional) integral $I$
of a well-behaved function $f(\boldsymbol{\theta})$. Then

$I=\int
f(\boldsymbol{\theta})\,d\boldsymbol{\theta}=\int\frac{f(\boldsymbol{\theta})}{g(\boldsymbol{\theta})}g(\boldsymbol{\theta})\,d(\boldsymbol{\theta})$

so if $g(\boldsymbol{\theta})$ is a probability density function and
$\boldsymbol{\theta}_{i}$ are independent samples from it, then

$I=\mbox{E}_{g}[f(\boldsymbol{\theta})/g(\boldsymbol{\theta})]\approx\frac{1}{n}\sum_{i=1}^{n}\frac{f(\boldsymbol{\theta}_{i})}{g(\boldsymbol{\theta}_{i})}=I_{n}.$
(4)

$I_{n}$ is an unbiased approximation to $I$ and by the central limit theorem
will tend to a normal distribution. It has variance

$\mbox{Var}[I_{n}]=\frac{1}{n}\int\left(\frac{f(\boldsymbol{\theta})}{g(\boldsymbol{\theta})}-I\right)^{2}g(\boldsymbol{\theta})\,d\boldsymbol{\theta}=\frac{1}{n}\int\frac{(f(\boldsymbol{\theta})-Ig(\boldsymbol{\theta}))^{2}}{g(\boldsymbol{\theta})}\,d\boldsymbol{\theta}$
(5)

Sometimes $f$ is called the _target_ and $g$ is called the _proposal_
distribution.

Assuming that $f$ is non-negative, then minimum variance (of $0!$) is achieved
when $g=f/I$—in other words when $g$ is just the normalized version of $f$.
This cannot be done in practice because normalizing $f$ requires knowing the
quantity $I$ that we wanted to approximate; however (5) is still important
because it means that the more similar the proposal is to the target, the
better our estimator $I_{n}$ becomes. In particular, $f$ must go to 0 faster
than $g$ or the estimator will have infinite variance.

To summarize this section:

  1. 1.

Importance sampling is a monte carlo integration technique which evaluates the
target using samples from a proposal distribution.

  2. 2.

The estimator is unbiased, normally distributed, and its variance (if not 0 or
infinity) decreases as $O(n^{-1})$ (using big-$O$ notation).

  3. 3.

The closer the proposal is to the target, the better the estimator. The
proposal also needs to have longer tails than the target.

###  2.2 Nonparametric Importance Sampling

A difficulty with importance sampling is that it is often difficult to choose
a proposal distribution $g$. Not enough is known about $f$ to choose an
optimal distribution, and if a bad distribution is chosen the result can have
large or even infinite variance. One approach to the selection of proposal $g$
is to use non-parametric techniques to build $g$ from samples of $f$. I call
this class of techniques self-importance sampling, or arrogance sampling for
short, because they attempt to sample $f$ from itself without using any
external information. (And also isn’t it a bit arrogant to try to evaluate a
complex, multidimensional integral using only the values at a few points?) The
method of this paper falls into this class and particularly deserves the name
because the target and proposal (when they are both non-zero) have exactly the
same values up to a multiplicative constant.

Two papers which apply nonparametric importance sampling to the problem of
marginal likelihood computation (or computation of normalizing constants) are
Zhang (1996) and Neddermeyer (2009). Although both authors apply their methods
to more general situations, here I will use the framework suggested by (3) and
assume that we can compute $p(\boldsymbol{\theta}\wedge\boldsymbol{x}|T)$ for
arbitrary $\boldsymbol{\theta}$ and also that we can sample from the posterior
parameter distribution $p(\boldsymbol{\theta}|\boldsymbol{x},T)$. The goal is
to estimate the normalizing constant, the marginal likelihood
$p(\boldsymbol{x}|T)$.

Zhang’s approach is to build the proposal $g$ using traditional kernel density
estimation. $m$ samples are first drawn from
$p(\boldsymbol{\theta}|\boldsymbol{x},T)$ and used to construct $g$. Then $n$
samples are drawn from $g$ and used to evaluate $p(\boldsymbol{x}|T)$ as in
traditional importance sampling. This approach is quite intuitive because
kernel estimation is a popular way of approximating an unknown function. Zhang
proves that the variance of his estimator decreases as
$O(m^{\frac{-4}{4+d}}n^{-1})$ where $d$ is the dimensionality of
$\boldsymbol{\theta}$, compared to $O(n^{-1})$ for standard (parametric)
importance sampling.

There were, however, a few issues with Zhang’s method:

  1. 1.

A kernel density estimate is equal to 0 at points far from the points the
kernel estimator was built on. This is a problem because importance sampling
requires the proposal to have longer tails than the target. This fact forces
Zhang to make the restrictive assumption that
$p(\boldsymbol{\theta}|\boldsymbol{x},T)$ has compact support.

  2. 2.

It is hard to compute the optimal kernel bandwidth. Zhang recommends using a
plug-in estimator because the function
$p(\boldsymbol{\theta}\wedge\boldsymbol{x}|T)$ is available, which is unusual
for kernel estimation problems. Still, bandwidth selection appears to require
significant additional analysis.

  3. 3.

Finally, although the variance may decrease as $O(m^{\frac{-4}{4+d}}n^{-1})$
as $m$ increases, the difficulty of computing $g(\boldsymbol{\theta})$ also
increases with $m$, because it requires searching through the $m$ basis points
to find all the points close to $\boldsymbol{\theta}$. In multiple dimensions,
this problem is not trivial and may outweigh the $O(m^{\frac{-4}{4+d}})$
speedup (in the worst case, practical evaluation of $g(\boldsymbol{\theta})$
at a single point may be $O(m)$). See Zlochin and Baram (2002) for some
discussion of these issues.

Neddermeyer (2009) uses a similar approach to Zhang and also achieves a
variance of $O(m^{\frac{-4}{4+d}}n^{-1})$. It improves on Zhang’s approach in
two ways relevant to this paper:

  1. 1.

The support of $p(\boldsymbol{\theta}|\boldsymbol{x},T)$ is not required to be
compact.

  2. 2.

Instead of using kernel density estimators, linear blend frequency polynomials
(LBFPs) are used instead. LBFPs are basically histograms whose density is
interpolated between adjacent bins. As a result, the computation of
$g(\boldsymbol{\theta})$ requires only finding which bin $\boldsymbol{\theta}$
is in, and looking up the histogram value at that and adjacent bins ($2^{d}$
bins in total).

As we will see in section 3, the arrogance sampling described in this paper is
similar to the methods of Zhang and Neddermeyer.

###  2.3 Harmonic Mean Estimator

The harmonic mean estimator is a simple and notorious method for calculating
marginal likelihoods. It is a kind of importance sampling, except the proposal
$g$ is actually the distribution
$p(\boldsymbol{\theta}|\boldsymbol{x},T)=p(\boldsymbol{\theta}\wedge\boldsymbol{x}|T)/p(\boldsymbol{x}|T)$
to be normalized and the target $f$ is the known distribution
$p(\boldsymbol{\theta}|T)$. Then if $\boldsymbol{\theta}_{i}$ are samples from
$p(\boldsymbol{\theta}|x,T)$, we apparently have

$1\approx\frac{1}{n}\sum_{i=1}^{n}\frac{p(\boldsymbol{\theta}_{i}|T)}{p(\boldsymbol{\theta}_{i}|\boldsymbol{x},T)}=\frac{1}{n}\sum_{i=1}^{n}\frac{p(\boldsymbol{\theta}_{i}|T)}{p(\boldsymbol{x}|\boldsymbol{\theta}_{i},T)p(\boldsymbol{\theta}_{i}|T)/p(\boldsymbol{x}|T)}=\frac{1}{n}\sum_{i=1}^{n}\frac{1}{p(\boldsymbol{x}|\boldsymbol{\theta}_{i},T)/p(\boldsymbol{x}|T)}$

hence

$p(\boldsymbol{x}|T)\stackrel{{\scriptstyle?}}{{\approx}}\left(\frac{1}{n}\sum_{i=1}^{n}\frac{1}{p(\boldsymbol{x}|\boldsymbol{\theta}_{i},T)}\right)^{-1}$
(6)

Two advantages of the harmonic mean estimator are that it is simple to compute
and only depends on samples from $p(\boldsymbol{\theta}|x,T)$ and the
likelihood $p(\boldsymbol{x}|\boldsymbol{\theta},T)$ at those samples. The
main drawback of the harmonic mean estimator is that it doesn’t work—as
mentioned earlier the importance sampling proposal distribution needs to have
longer tails than the target. In this case the target
$p(\boldsymbol{\theta}|T)$ typically has longer tails than the proposal
$p(\boldsymbol{\theta}|\boldsymbol{x},T)$ and thus (6) has infinite variance.
Despite not working, the harmonic mean estimator continues to be popular
(Neal, 2008).

##  3 Description of Technique

This paper’s arrogance sampling technique is a simple method that applies the
nonparametric importance techniques of Zhang and Neddermeyer in an attempt to
develop a method almost as convenient as the harmonic mean estimator.

The only required inputs are samples $\boldsymbol{\theta}_{i}$ from
$p(\boldsymbol{\theta}|\boldsymbol{x},T)$ and the values
$p(\boldsymbol{\theta}_{i}\wedge\boldsymbol{x}|T)=p(x|\boldsymbol{\theta}_{i},T)p(\boldsymbol{\theta}_{i}|T)$.
This is similar to the harmonic mean estimator, but perhaps slightly less
convenient because $p(\boldsymbol{\theta}_{i}\wedge\boldsymbol{x}|T)$ is
required instead of $p(\boldsymbol{x}|\boldsymbol{\theta}_{i},T)$.

There are two basic steps:

  1. 1.

Take $m$ samples from $p(\boldsymbol{\theta}|\boldsymbol{x},T)$ and using
modified histogram density estimation, construct probability density function
$f(\boldsymbol{\theta})$.

  2. 2.

With $n$ more samples from $p(\boldsymbol{\theta}|\boldsymbol{x},T)$, estimate
$1/p(\boldsymbol{x}|T)$ via importance sampling with target $f$ and proposal
$p(\boldsymbol{\theta}|\boldsymbol{x},T)$.

These steps are described in more detail below.

###  3.1 Construction of the Histogram

Of the $N$ total samples $\boldsymbol{\theta}_{i}$ from
$p(\boldsymbol{\theta}|\boldsymbol{x},T)$, the first $m$ will be used to make
a histogram. The optimal choice of $m$ will be discussed below, but in
practice this seems difficult to determine. An arbitrary rule of
$\mbox{min}(0.2N,2\sqrt{N})$ can be used in practice.

With a traditional histogram, the only available information is the location
of the sampled points. In this case we also know the (scaled) heights
$p(\boldsymbol{\theta}\wedge\boldsymbol{x}|T)$ at each sampled point. We can
use this extra information to improve the fit.

Our “arrogant” histogram $f$ is constructed the same as a regular histogram,
except the bin heights are not determined by the number of points in each bin,
but rather by the minimum density over all points in the bin. If a bin
contains no sampled points, then $f(\boldsymbol{\theta})=0$ for
$\boldsymbol{\theta}$ in that bin. Then $f$ is normalized so that $\int
f(\boldsymbol{\theta})\,d\boldsymbol{\theta}=1$.

To determine our bin width, we can simply and somewhat arbitrarily set our bin
width $h$ so that the histogram is positive for 50% of the sampled points from
the distribution $p(\boldsymbol{\theta}|\boldsymbol{x},T)$. To approximate
$h$, we can use a small number of samples (say, 40) from
$p(\boldsymbol{\theta}|\boldsymbol{x},T)$ and set $h$ so that
$f(\boldsymbol{\theta})>0$ for exactly half of these samples.

Figure 1 compares the traditional and new histograms for a one dimensional
normal distribution based on 50 samples. The green rug lines indicate the $50$
sampled points which are the same for all. The arrogant histogram’s bin width
is chosen as above. The traditional histogram’s optimal bin width was
determined by Scott’s rule to minimize mean squared error. As the figure
shows, the modified histogram is much smoother for a given bin width, so a
smaller bin width can be used. On the other hand, $f$ will either equal 0 or
have about twice the original density at each point, while the traditional
histogram’s density is numerically close to the original density.

Figure 1: Histogram Comparison

###  3.2 Importance Sampling

The remaining $n=N-m-40$ sampled points can be used for importance sampling.
Using equation (4) with histogram $f$ as our target and
$p(\boldsymbol{\theta}|\boldsymbol{x},T)$ as the proposal, we have

$1\approx
I_{n}=\frac{1}{n}\sum_{i=1}^{n}\frac{f(\boldsymbol{\theta}_{i})}{p(\boldsymbol{\theta}_{i}|\boldsymbol{x},T)}=\frac{1}{n}\sum_{i=1}^{n}\frac{f(\boldsymbol{\theta}_{i})}{p(\boldsymbol{\theta}_{i}\wedge\boldsymbol{x}|T)/p(\boldsymbol{x}|T)}$

hence

$p(\boldsymbol{x}|T)\approx
p(\boldsymbol{x}|T)/I_{n}=\left(\frac{1}{n}\sum_{i=1}^{n}\frac{f(\boldsymbol{\theta}_{i})}{p(\boldsymbol{\theta}_{i}\wedge\boldsymbol{x}|T)}\right)^{-1}=A_{n}$
(7)

To underscore the self-important/arrogant nature of this approximation
$A_{n}$, we can rewrite (7) as

$p(\boldsymbol{x}|T)\approx
H\left(\frac{1}{n}\sum_{i=1}^{n}\frac{\mbox{min}\\{p(\boldsymbol{\theta}_{j}\wedge\boldsymbol{x}|T):\boldsymbol{\theta}_{j}\mbox{
and }\boldsymbol{\theta}_{j}\mbox{ are in the same
bin}\\}}{p(\boldsymbol{\theta}_{i}\wedge\boldsymbol{x}|T)}\right)^{-1}$

where $H$ is the histogram normalizing constant. This equation shows that all
the values in the numerator and the denominator of our importance sampling are
from the same distribution $p(\boldsymbol{\theta}\wedge\boldsymbol{x}|T)$.

Note that the histogram $f$ is the target of the importance sampling and
$p(\boldsymbol{\theta}\wedge\boldsymbol{x}|T)$ is the proposal. This is
backwards from the usual scheme where the unknown distribution is the target
and the known distribution is the proposal. Instead here the unknown
distribution is the proposal, as in the harmonic mean estimator (see Robert
and Wraith (2009) for another example of this.)

As in section 2.1, our approximation of $p(\boldsymbol{x}|T)^{-1}$ tends to a
normal distribution as $n\to\infty$ by the central limit theorem. This fact
can be used to estimate a confidence interval around $p(\boldsymbol{x}|T)$.

##  4 Validity of Method

This section will investigate the performance of the method. First, note that
this method is just an implementation of importance sampling, so $A_{n}^{-1}$
should converge to $p(\boldsymbol{x}|T)^{-1}$ with finite variance as long as
the proposal density $p(\boldsymbol{\theta}|\boldsymbol{x},T)$ exists and is
finite and positive on the compact region where the target histogram density
is positive.

To calculate the speed of convergence we will use equation (5) where $f$ is
the histogram,
$g(\boldsymbol{\theta})=p(\boldsymbol{\theta}|\boldsymbol{x},T)$, and $I=1$
because the histogram has been normalized. Unless otherwise noted, we will
assume below that $g:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is finite, twice
differentiable and positive, and that $\int\frac{\lVert\nabla\cdot
g(\boldsymbol{\theta})\rVert^{2}}{g(\boldsymbol{\theta})}d\boldsymbol{\theta}$
is finite.

###  4.1 Histogram Bin Width

One important issue will be how quickly the $d$-dimensional histogram’s
selected bin width $h$ goes to 0 as the number of samples
$m\rightarrow\infty$. This section will only offer an intuitive argument. For
any $m$, the histogram will enclose about the same probability ($\frac{1}{2}$)
and will have about the same average density in a fixed region. Each bin has
volume $h^{d}$, so if $l$ is the number of bins then $lh^{d}=O(1)$ and
$h\propto l^{-d}$.

Furthermore, the distribution of the sampled points converges to the actual
distribution $g(\boldsymbol{\theta})$. If $m>O(l)$, an unbounded number of
sampled points would end up in each bin. If $m<O(l)$, then some bins would
have no points in them. Neither of these is possible because exactly one
sampled point is necessary to establish each bin. Thus $m\propto l$ and
$h\propto m^{-d}$.

###  4.2 Conditional Variance

Before estimating the convergence rate of $A_{n}$ we will prove something
about the conditional variance of importance sampling. Let
$A=\\{\boldsymbol{\theta}:f(\boldsymbol{\theta})>0\\}$, $\mathbf{1}_{A}$ be
the characteristic function of $A$, and
$q=\int_{A}g(\boldsymbol{\theta})\,d\boldsymbol{\theta}$. Define

$g_{A}(\boldsymbol{\theta})=\left\\{\begin{array}[]{cl}g(\boldsymbol{\theta})/q&\mbox{
if }\boldsymbol{\theta}\in A\\\ 0&\mbox{ otherwise }\end{array}\right.$

Then $g_{A}$ is the density of $g$ conditional on $f>0$. Define
$\mbox{Var}_{A}$ and $\mbox{E}_{A}$ to mean the variance and expectation
conditional on $f(\boldsymbol{\theta})>0$. Thus

$\displaystyle\mbox{Var}(f(\boldsymbol{\theta})/g(\boldsymbol{\theta}))$
$\displaystyle=$
$\displaystyle\mbox{Var}(\mbox{E}(f(\boldsymbol{\theta})/g(\boldsymbol{\theta})|\mathbf{1}_{A}))+\mbox{E}(\mbox{Var}(f(\boldsymbol{\theta})/g(\boldsymbol{\theta})|\mathbf{1}_{A}))$
$\displaystyle=$
$\displaystyle\mbox{Var}\left(\begin{array}[]{cl}\mbox{E}_{A}(f(\boldsymbol{\theta})/g(\boldsymbol{\theta}))&\mbox{
if }\boldsymbol{\theta}\in A\\\ 0&\mbox{ otherwise}\\\ \end{array}\right)$
$\displaystyle+\,\mbox{E}\left(\begin{array}[]{cl}\mbox{Var}_{A}(f(\boldsymbol{\theta})/g(\boldsymbol{\theta}))&\mbox{
if }\boldsymbol{\theta}\in A\\\ 0&\mbox{ otherwise}\\\ \end{array}\right)$
$\displaystyle=$ $\displaystyle\mbox{Var}\left(\begin{array}[]{cl}1/q&\mbox{
if }\boldsymbol{\theta}\in A\\\ 0&\mbox{ otherwise}\\\
\end{array}\right)+q\mbox{Var}_{A}(f(\boldsymbol{\theta})/g(\boldsymbol{\theta}))$
$\displaystyle=$
$\displaystyle(1/q)^{2}q(1-q)+\frac{1}{q}\mbox{Var}_{A}(f(\boldsymbol{\theta})/qg(\boldsymbol{\theta}))$
$\displaystyle=$
$\displaystyle\frac{1-q}{q}+\frac{1}{q}\mbox{Var}_{A}(f(\boldsymbol{\theta})/g_{A}(\boldsymbol{\theta}))$

We will assume below that $q=\frac{1}{2}$, so that

$\mbox{Var}(f(\boldsymbol{\theta})/g(\boldsymbol{\theta}))=1+2\mbox{Var}_{A}(f(\boldsymbol{\theta})/g_{A}(\boldsymbol{\theta}))$
(11)

###  4.3 Importance Sampling Convergence

With $f$, $g$, and $A$ as defined above, $f$ and $g_{A}$ have the same domain.
Assuming errors in estimating $q$ and normalization errors are of a lesser
order of magnitude, we can treat the histogram heights as being sampled from
$g_{A}$. Suppose the histogram has $l$ bins $\\{B_{j}\\}$, each with width $h$
and based around the points $g_{A}(\boldsymbol{\theta}_{j})$. Then by equation
(5),

$\displaystyle\mbox{Var}_{A}(f(\boldsymbol{\theta})/g_{A}(\boldsymbol{\theta}))$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{l}\int_{B_{j}}\frac{(f(\boldsymbol{\theta})-g_{A}(\boldsymbol{\theta}))^{2}}{g_{A}(\boldsymbol{\theta})}d\boldsymbol{\theta}$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{l}\int_{B_{j}}\frac{(g_{A}(\boldsymbol{\theta})+\nabla
g_{A}(\boldsymbol{\theta})\cdot(\boldsymbol{\theta}_{j}-\boldsymbol{\theta})+O((\boldsymbol{\theta}_{j}-\boldsymbol{\theta})^{2})-g_{A}(\boldsymbol{\theta}))^{2}}{g_{A}(\boldsymbol{\theta})}d\boldsymbol{\theta}$
$\displaystyle=$ $\displaystyle\sum_{j=1}^{l}\int_{B_{j}}\frac{(\nabla
g_{A}(\boldsymbol{\theta})\cdot(\boldsymbol{\theta}_{j}-\boldsymbol{\theta}))^{2}+O((\boldsymbol{\theta}_{j}-\boldsymbol{\theta})^{3})}{g_{A}(\boldsymbol{\theta})}d\boldsymbol{\theta}$
$\displaystyle\leq$
$\displaystyle\sum_{j=1}^{l}\int_{B_{j}}\frac{\lVert\nabla\cdot
g_{A}(\boldsymbol{\theta})\rVert^{2}h^{2}}{g_{A}(\boldsymbol{\theta})}d\boldsymbol{\theta}$
$\displaystyle=$ $\displaystyle h^{2}\int\frac{\lVert\nabla\cdot
g_{A}(\boldsymbol{\theta})\rVert^{2}}{g_{A}(\boldsymbol{\theta})}d\boldsymbol{\theta}$

Because $h\propto m^{-d}$ where $d$ is the number of dimensions, and $m$ is
the number of samples used to make the histogram,

$\mbox{Var}_{A}(f(\boldsymbol{\theta})/g_{A}(\boldsymbol{\theta}))\leq
Cm^{-2/d}\\\ $

where $C\propto\int\frac{\lVert\nabla\cdot
g_{A}(\boldsymbol{\theta})\rVert^{2}}{g_{A}(\boldsymbol{\theta})}d\boldsymbol{\theta}$.
Putting this together with (11), we get

$\mbox{Var}(I_{n})=\mbox{Var}(p(\boldsymbol{x}|T)/A_{n})=n^{-1}(1+O(Cm^{-2/d}))$
(12)

##  5 Implementation Issues

###  5.1 Speed of Convergence

The variance of $n^{-1}(1+O(Cm^{-2/d}))$ given by (12) is asymptotically equal
to $n^{-1}$, which is the typical importance sampling rate. In practice
however, the asymptotic results cannot distinguish useful from impractical
estimators. If $Cm^{-2/d}$ is small and
$\mbox{Var}(p(\boldsymbol{x}|T)/A_{n})\approx n^{-1}$, then
$p(\boldsymbol{x}|T)$ can be approximated in only 1000 samples to about
$6\%=\frac{1.96}{\sqrt{1000}}$ with 95% confidence. For many theory choice
purposes, this is quite sufficient. Thus in typical problem cases the factor
of $Cm^{-2/d}$ will be very significant. If $Cm^{-2/d}\gg 1$, then the
convergence rate may in practice be similar to $n^{-1}m^{-2/d}$. Compare this
to the rate of $n^{-1}m^{-4/(4+d)}$ for the methods proposed by Zhang and
Neddermeyer.

This method also uses simple histograms, instead of a more sophisticated
density estimation method (Zhang uses kernel estimation, Neddermeyer uses
linear blend frequency polynomials). Although simple histograms converge
slower for large $d$ as shown above, they are much faster to compute for large
$d$.

Neddermeyer’s LBFP algorithm is quite efficient compared to Zhang’s, but its
running time is $O(2^{d}d^{2}n^{\frac{d+5}{d+4}})$. $d$ is a constant for any
fixed problem, but if, say, $d=10$, then the dimensionality constant
multiplies the running time by $2^{10}10^{2}\approx 10^{5}$.

By contrast, this paper’s method takes only $O(dm\mbox{log}(m))$ time to
construct the initial histogram, and an additional $O(dn\mbox{log}(m))$ time
to do the importance sampling. The main reason for the difference is that
querying a simple histogram can be done in $\mbox{log}(m)$ time by computing
the bin coordinates and looking up the bin’s height in a tree structure.
However, querying a LBFP requires blending all nearby bins and is thus
exponential in $d$.

###  5.2 When $g=0$

Our discussion assumed that
$g(\boldsymbol{\theta})=p(\boldsymbol{\theta}|\boldsymbol{x},T)$ was always
positive. If $g$ goes to 0 where the histogram is positive, the variance of
$A_{n}^{-1}$ will be infinite. However, this paper’s method can still be used
if $g(\boldsymbol{\theta})$ is 0 over some well-defined area.

For instance, suppose one dimension $\theta_{k}$ of $p(\boldsymbol{\theta}|T)$
is defined by a gamma distribution, so that $p(\theta_{k}|T)=0$ if and only if
$\theta_{k}\leq 0$. Then we can ensure the variance is not infinite by
checking that the histogram is only defined where $\theta_{k}>\epsilon>0$ for
some fixed $\epsilon$.

The margLikArrogance package contains a simple mechanism to do this. The user
may specify a range along each dimension of $\boldsymbol{\theta}$ where it is
known that $g>0$. If the histogram is non-zero outside of this range, the
method aborts with an error.

Note that the variance of the estimator increases with
$\int\frac{\lVert\nabla\cdot
g_{A}(\boldsymbol{\theta})\rVert^{2}}{g_{A}(\boldsymbol{\theta})}d\boldsymbol{\theta}$.
In practice the estimator will work well only when $g$ doesn’t go to 0 too
quickly where the histogram is positive. In these cases the histogram will be
defined well away from any region where $g=0$ and infinite variance won’t be
an issue even if $g=0$ somewhere.

###  5.3 Bin Shape

Cubic histogram bins were used above—their widths were fixed at $h$ in each
dimension. Although the asymptotic results aren’t affected by the shape of
each bin, for usable convergence rates the bins’ dimensions need to compatible
with the shape of the high probability region of
$p(\boldsymbol{\theta}|\boldsymbol{x},T)$. Unfortunately, it is difficult to
determine the best bin shapes.

The margLikArrogance package contains a simple workaround: by default the
distribution is first scaled so that the sampled standard deviation along each
dimension is constant. This is equivalent to setting each bin’s width by
dimension in proportion to that dimension’s standard deviation. If this simple
rule of thumb is insufficient, the user can scale the sampled values of
$p(\boldsymbol{\theta}|\boldsymbol{x},T)$ manually (and make the corresponding
adjustment to the estimate $A_{n}$).

##  6 Conclusion

This paper has described an “arrogance sampling” technique for computing the
marginal likelihood or Bayes factor of a Bayesian model. It involves using
samples from the model’s posterior parameter distribution along with the
scaled values of the distribution’s density at those points. These samples are
divided into two main groups: $m$ samples are used to build a histogram; $n$
are used to importance sample the histogram using the posterior parameter
distribution as the proposal.

This method is simple to implement and runs quickly in
$O(d(m+n)\mbox{log}(m))$ time. Its asymptotic convergence rate,
$n^{-1}(1+O(Cm^{-2/d}))$, is not remarkable, but in practice convergence is
fast for many problems. Because the required inputs are similar to those of
the harmonic mean estimator, it may be a convenient replacement for it.

##  7 References

  1. 1.

S. Chib. “Marginal Likelihood from the Gibbs Output” _Journal of the American
Statistical Association_. Vol 90, No 432. (1995)

  2. 2.

S. Chib and I. Jeliazkov. “Accept-reject Metropolis-Hastings sampling and
marginal likelihood estimation” _Statistica Neerlandica_. Vol 59, No 1. (2005)

  3. 3.

A. Gelman and X. Meng. “Simulating Normalizing Constants: From Importance
Sampling to Bridge Sampling to Path Sampling” _Statistical Science_. Vol 13,
No 2. (1998)

  4. 4.

R. Kass and A. Raftery. “Bayes Factors” _Journal of the American Statistical
Association_. Vol 90, No 430. (1995)

  5. 5.

R. Neal. “The Harmonic Mean of the Likelihood: Worst Monte Carlo Method Ever”.
Blog post, http://radfordneal.wordpress.com/2008/08/17/the-harmonic-mean-of-
the-likelihood-worst-monte-carlo-method-ever/. (2008)

  6. 6.

J. Neddermeyer. “Computationally Efficient Nonparametric Importance Sampling”
_Journal of the American Statistical Association_. Vol 104, No 486. (2009)
arXiv:0805.3591v2

  7. 7.

A. Owen and Y. Zhou. “Safe and effective importance sampling” _Journal of the
American Statistical Association_. Vol 95, No 449. (2000)

  8. 8.

C. Robert and D. Wraith. “Computational methods for Bayesian model choice”
arXiv:0907.5123v1

  9. 9.

P. Zhang. “Nonparametric Importance Sampling” _Journal of the American
Statistical Association_. Vol 91, No 435. (1996)

  10. 10.

M. Zlochin and Y. Baram. “Efficient Nonparametric Importance Sampling for
Bayesian Inference” _Proceedings of the 2002 International Joint Conference on
Neural Networks_ 2498–2502. (2002)