Ant-Inspired Density Estimation via Random Walks
Cameron Musco
cnmusco@mit.edu

Hsin-Hao Su
hsinhao@mit.edu

Nancy Lynch
lynch@csail.mit.edu

arXiv:1603.02981v1 [cs.DC] 9 Mar 2016

Massachusetts Institute of Technology
March 10, 2016
Abstract
Many ant species employ distributed population density estimation in applications ranging from
quorum sensing [Pra05], to task allocation [Gor99], to appraisal of enemy colony strength [Ada90]. It has
been shown that ants estimate density by tracking encounter rates – the higher the population density,
the more often the ants bump into each other [Pra05, GPT93].
We study distributed density estimation from a theoretical perspective. We show that a group of
anonymous agents randomly walking on a grid are able to estimatetheir density d to within a multiplicasteps by measuring their encounter
tive factor 1 ± ǫ with probability 1 − δ in just Õ log(1/δ)dǫlog(1/dǫ)
2
rates with other agents. Despite dependencies inherent in the fact that nearby agents may collide repeatedly (and, worse, cannot recognize when this happens), this bound nearly matches what is required
to estimate d by independently sampling grid locations.
From a biological perspective, our work helps shed light on how ants and other social insects can
obtain relatively accurate density estimates via encounter rates. From a technical perspective, our
analysis provides new tools for understanding complex dependencies in the collision probabilities of
multiple random walks. We bound the strength of these dependencies using local mixing properties
of the underlying graph. Our results extend beyond the grid to more general graphs and we discuss
applications to social network size estimation, density estimation by robot swarms, and random walkedbased sampling of sensor networks.

1

Introduction

The ability to sense local population density is an important tool used by many ant species. When a colony
must relocate to a new nest, scouts search for potential nest sites, assess their quality, and recruit other scouts
to high quality locations. A high enough density of scouts at a potential new nest (a quorum threshold )
triggers those ants to decide on the site and transport the rest of the colony there [Pra05]. When neighboring
colonies compete for territory, a high relative density of a colony’s ants in a contested area will cause those
ants to attack enemies in the area, while a low relative density will cause the colony to retreat [Ada90].
Varying densities of ants successfully performing certain tasks such as foraging or brood care can trigger
other ants to switch tasks, maintaining proper worker allocation within in the colony [Gor99, SHG06].
It has been shown that ants estimate density in a distributed manner, by measuring encounter rates
[Pra05, GPT93]. As ants randomly walk around an area, if they bump into a larger number of other ants,
this indicates a higher population density. By tracking encounters with specific types of ants, e.g. successful
foragers or enemies, ants can estimate more specific densities. This strategy allows each ant to obtain an
accurate density estimate and requires very little communication – ants must simply detect when they collide
and do not need to perform any higher level data aggregation.

1.1

Density Estimation on the Grid

We study distributed density estimation from a theoretical perspective. We model a colony of ants as a set
of anonymous agents randomly distributed on a two-dimensional grid. Computation proceeds in rounds,
with each agent stepping in a random direction in each round. A collision occurs when two agents reach
the same position in the same round and encounter rate is measured as the number of collisions an agent is
involved in during a sequence of rounds divided by the number of rounds. Aside from collision detection,
the agents have no other means of communication.
The intuition that encounter rate tracks density is clear. It is easy to show that, for a set of randomly
walking agents, the expected encounter rate measured by each agent is exactly the density d – the number
of agents divided by the grid size (see Lemma 3). However, it is unclear if encounter rate actually gives a
good density estimate – i.e., if it concentrates around its expectation.
Consider agents positioned not on the grid, but on a complete graph. In each round, each agent steps to
a uniformly random position and in expectation, the number of other agents they collide with in this step
is d. Since each agent chooses its new location uniformly at random in each step, collisions are essentially
independent between rounds. The agents are effectively taking
Bernoulli samples with success

 independent
rounds
obtain a (1 ± ǫ) multiplicative
probability d, and by a standard Chernoff bound, within O log(1/δ)
dǫ2
approximation to d with probability 1 − δ.
On the grid graph, the picture is significantly more complex. If two agents are initially located near each
other on the grid, they are more likely to collide via random walking. After a first collision, due to their
proximity, they are likely to collide repeatedly in future rounds. The agents cannot recognize repeat collisions
since they are anonymous and even if they could, it is unclear that it would help. On average, compared to
the complete graph, agents collide with fewer individuals and collide multiple times with those individuals
that they do encounter, causing an increase in encounter rate variance and making density estimation more
difficult.
Mathematically speaking, on a graph with a fast mixing time [Lov93], like the complete graph, each
agent’s location is only weakly correlated with its previous locations. This ensures that collisions are also
weakly correlated between rounds and encounter rate serves as a very accurate estimate of density. The grid
graph on the other hand is slow mixing – agent positions and hence collisions are highly correlated between
rounds. This correlation increases encounter rate variance.

1.2

Our Contributions

Surprisingly, despite this increased variance, encounter
rate-based density estimation
on the grid is nearly as


log(1/δ) log log(1/δ) log(1/dǫ)
rounds, each agent’s encounter
accurate as on the complete graph. After just O
dǫ2
rate is a (1 ± ǫ) approximation to d with probability 1 − δ (see Theorem 2).
1

Technically, to bound accuracy on the grid, we obtain moment bounds on the number of times that two
randomly walking agents repeatedly collide over a set of rounds. These bounds also apply to the number of
equalizations (returns to starting location) of a single walk. While expected random walk hitting times, return
times, and collision rates are well understood [Lov93, ES09], higher moment bounds and high probability
results are much less common. We hope our bounds are of general use in the theoretical study of random
walks and random-walk based algorithms.
Our moment bounds show that, while the grid graph is slow mixing, it has sufficiently strong local mixing
to make random walk-based density estimation accurate. Random walks tend to spread quickly over a local
area and not repeatedly cover the same nodes. Significant work has focused on showing that random walk
sampling is nearly as good as independent sampling for fast mixing expander graphs [Gil98, CLLM12]. We
are the first to extend this type of analysis to slowly mixing graphs, showing that strong local mixing is
sufficient in many applications.
Beyond the grid, we show how to generate moment bounds from a bound on the probability that two
random walks re-collide (or analogously, that a single random walk equalizes) after a certain number of steps,
and demonstrate application of this technique to d-dimensional grids, regular expanders, and hypercubes.
We discuss applications of our results to social network size estimation via random walk [KLS11], obtaining
significant improvements over known work for networks with relatively slow global mixing times, but strong
local mixing. We also discuss connections to density estimation by robot swarms and random walk-based
sensor network sampling [AB04, LB07].

1.3

Road Map

In Section 2 we overview our theoretical model for distributed density estimation on the grid.
In Section 3, as a warm up, we give a simple density estimation algorithm that does not employ random
walks, is easy to analyze, but is not biologically plausible.
In Section 4 we give our main technical results on random walk-based density estimation.
In Section 5 we show how to extend our bounds to a number of graphs other than the grid.
In Section 6 we discuss applications of our results to social network size estimation, robot swarm, and sensor
network algorithms.

2

Theoretical Model for Density Estimation

√
√
We consider a two-dimensional torus with A nodes (dimensions A × A) populated with identical anonymous agents. We assume that A is large – larger than the area agents traverse over the runtimes of our
algorithms. We feel that this torus model successfully captures the dynamics of density estimation on a
surface, while avoiding complicating factors of boundary behavior.
Initially each agent is placed independently at a uniform random node in the torus. Computation proceeds
in discrete, synchronous rounds. In each round an agent may either remain in its current location or step
to any of its four neighboring grid squares. Formally, each agent has an ordered pair position which it may
update in each round by adding step s ∈ {(0, 1), (0, −1), (1, 0), (−1, 0), (0, 0)}.
A randomly walking agent chooses s uniformly at random from {(0, 1), (0, −1), (1, 0), (−1, 0)} in each
round. Of course, in reality ants do not move via pure random walk – observed encounter rates seem to
actually be lower than predicted by a pure random walk model [GPT93, NTD05]. However, we feel that our
model sufficiently captures the highly random movement of ants while remaining tractable to analysis and
applicable to ant-inspired random walk-based algorithms (Section 6).
Aside from the ability to move in each round, agents can sense the number of agents other than themselves
at their position at the end of each round, formally through calling count(position). We say that two agents
collide in round r if they have the same position at the end of the round. Outside of collision counting,
agents have no means of communication. They are anonymous (cannot uniquely identify each other) and all
execute identical density estimation routines.

2

Density Estimation Problem Let (n + 1) be the number of agents and define population density as
def

d = n/A. Each agent’s goal is to hestimate d to (1 ± ǫ) accuracy
with probability 1 − δ for ǫ, δ ∈ (0, 1) –
i
i.e., to return an estimate d˜ with P d˜ ∈ [(1 − ǫ)d, (1 + ǫ)d] ≥ 1 − δ. As a technicality, with n + 1 agents we

define d = n/A instead of d = (n + 1)/A for convenience of calculation. In the natural case, when n is large,
the distinction is minor.

3

Density Estimation via Simulation of Independent Sampling

As discussed, the challenge in analyzing random walk-based density estimation arises from increased variance
due to repeated collisions of nearby agents. Here we show that, if not restricted to random walking, agents
can avoid collision correlations by splitting into ‘stationary’ and ‘mobile’ groups and only counting collisions
between members of different groups. This allows them to essentially simulate independent sampling of gird
locations to estimate density. This algorithm is not ‘natural’ in a biological sense, however it is easy to
analyze, and demonstrates the feasibility of density estimation by anonymous agents on the grid. We give
pseudocode in Algorithm 1.
Algorithm 1 Independent Sampling-Based Density Estimation
input: runtime t
Set c := 0 and with probability 1/2, state := walking, else state := stationary.
for r = 1, ..., t do
if state := walking then
position := position + (0, 1)
⊲ Deterministic walk step.
end if
c := c + count(position)
⊲ Update collision count.
end for
return d˜ = 2c
t

3.1

Independent Sampling Accuracy Bound

We now present our main accuracy bound for the independent sampling algorithm.
√
Theorem 1 (Independent Sampling Accuracy Bound). After running for t rounds, assuming t < A and
d ≤ 1, 
Algorithm 
1 returns d˜ such that, for any δ > 0, with probability ≥ 1 − δ, d˜ ∈ [(1 − ǫ)d, (1 + ǫ)d] for
q


log(1/δ)
log(1/δ)
, d˜ is a (1 ± ǫ) multiplicative
.
In
other
words,
for
any
ǫ,
δ
∈
(0,
1)
if
t
=
Θ
ǫ =O
2
td
dǫ

estimate of d with probability ≥ 1 − δ.

Proof. Our analysis is from the perspective of an agent with state = walking. By symmetry, the distribution
of d˜ is identical for walking and stationary agents, so considering this case is sufficient.
Initially, assume that no two walking agents start in the same location. Given this assumption, we know
that a walking agent never collides with another walking agent – by assumption they all start in different
positions and update these positions identically in each round. In the written implementation agents always
step up, however any fixed pattern (e.g. a spiral) suffices.
In t steps, a walking agent visits t unique squares. Each of the n other agents is located in this set of
t
squares and stationary with probability 2A
. Further, each of these events is entirely independent from the
rest, as the agents are positioned and choose their state independently. So, for a walking agent, c is just a
t
t
. Clearly, E c = n · 2A
= td
sample of n independent random coin flips, each with success probability 2A
2 so
2
E d˜ = E c t = d. Further, by a Chernoff bound, for any ǫ ∈ (0, 1), the probability that d˜ is not a (1 ± ǫ)
multiplicative estimate of d is:
h
i
2
2
δ = P |d˜ − d| ≥ ǫd = P [|c − E c| ≥ ǫ E c] ≤ 2e−ǫ E c/3 ≤ 2e−ǫ td/6 .
3

This gives: log(1/δ) ≥ ǫ2 td/6 so ǫ = O

q

log(1/δ)
td


, yielding the result.

We now remove the assumption that no two walking agents start in the same location. We slightly
modify the algorithm – each agent sets c := c (mod t) before returning d˜ = 2c
t . If an agent starts alone and
is involved in < t collisions, this operation has no effect – the above bound holds.
If a walking agent is involved in < t ‘true collisions’ but starts in the same position as w ≥ 1 other
walking agents, the agents move in lockstep throughout the algorithm and are involved in w · t ‘spurious
collisions’ (w in each round). Setting c := c (mod t) exactly corrects for these spurious collisions and since
c now only includes collisions with stationary agents, the bound above holds.
Finally, if an agent is involved in ≥ t true collisions, this modification cannot worsen their estimate. If
c ≥ t and the agent does not set c := c (mod t), they compute d˜ ≥ 2tt ≥ 2. For ǫ < 1, the agent fails since
d ≤ 1. So setting c := c (mod t) can only increase probability of success.

4

Density Estimation via Random Walk Collision Rates

In Algorithm 1, each pair of agents can only collide once, at a specific location – the starting position of the
stationary agent in the pair. Collisions are independent and it is easy to show that the number of collisions
(and hence the density estimate) concentrates around its expectation. However, as discussed, independence
of collisions is unnecessary! Algorithm 2 describes a simple random walk-based approach that gives a nearly
matching bound.
Algorithm 2 Random Walk Encounter Rate-Based Density Estimation
input: runtime t
c := 0
for r = 1, ..., t do
position := position + rand{(0, 1), (0, −1), (1, 0), (−1, 0)}
c := c + count(position)
end for
return d˜ = ct

4.1

⊲ Random walk step.
⊲ Update collision count.

Random Walk-Based Density Estimation Analysis

Our main result follows; its proof appears at the end of Section 4 after some preliminary lemmas.
Theorem 2 (Random Walk Sampling Accuracy Bound). After running for t rounds, assuming t ≤ A,
Algorithm
2 returns d˜ such that, for any δ > 0, with probability ≥ 1 − δ, d˜ ∈ [(1 − ǫ)d, (1 + ǫ)d] for
q


log(1/dǫ)
log(1/δ) log(t)
, d˜ is a (1 ± ǫ)
. In other words, for any ǫ, δ ∈ (0, 1) if t = Θ log(1/δ) log log(1/δ)
ǫ=
td
dǫ2
multiplicative estimate of d with probability ≥ 1 − δ.
Throughout our analysis, we take the viewpoint of a single agent executing Algorithm 2, referred to as
‘agent a’. To start, we show that the encounter rate d˜ is an unbiased estimator of d:
Lemma 3 (Unbiased Estimator). E d˜ = d.
Proof. We can decompose c as the sum of collisions with different agents over different rounds. Specifically,
give the n other agents arbitrary ids 1, 2, ..., n and let cj (r) equal 1 if agent a collides with agent j in round
Pn Pt
r, and 0 otherwise. By linearity of expectation: E c = j=1 r=1 E cj (r).
Since each agent is initially at a uniform random location and after any number of steps, is still at uniform
random location, for all j, r, E cj (r) = 1/A. Thus, E c = nt/A = dt and E d˜ = E c/t = d.
With Lemma 3 in place, we now must show that the encounter rate is close to its expectation with high
probability and hence provides a good estimate of density.

4

4.2

Bounding the Effects of Repeat Collisions

P
Let cj = tr=1 cj (r) be the total number of collisions with agent j. Due to the initial uniform distribution
of the agents, the cj ’s are all independent and identically distributed.
Each cj is the sum of highly correlated random variables – due to the slow mixing of the grid, if two
agents collide at round r, they are much more likely to collide in successive rounds. However, by bounding
the strength of this correlation, we are able to give strong bounds
P on the moments of the distribution of
each cj , showing that it is sub-exponential. It follows that d˜ = 1t nj=1 cj , is also sub-exponential and hence
concentrates strongly around its expectation, the true density d.
We first bound the probability of a re-collision in round r + m, assuming a collision in round r:
Lemma 4 (Re-collision Probability
Consider two agents a1 and a2 randomly walking on a two√ Bound).
√
dimensional torus of dimensions A × A. If a1 and
in round r, for any m ≥ 0, the
 a2 collide again

1
1
probability that a1 and a2 collide in round r + m is Θ m+1 + O A .

Proof. From round r to round r + m, a1 and a2 take 2m random steps in total. Let Mx be the total number
of steps they take in the x direction and My be the total number in the y direction. Mx + My = 2m.
We start by computing the probability that the agents collide in round r + m conditioned on the values
of Mx and My . All steps are chosen independently, so we can consider movement in the x and y directions
separately. Specifically, let C be the event that the a1 and a2 collide in round r + m, Cx be the event that
they have the same x position, and Cy be the event that they have the same y position. We have:
P [C|Mx = mx , My = my ] = P [Cx |Mx = mx ] · P [Cy |My = my ] .

(1)

We first consider P [Cx |Mx = mx ]. All bounds will hold symmetrically for the y dimension. We split our
analysis into two cases. Let Cx1 be the event that the two agents have the same x position after round r + m
and have identical displacements from their starting locations. Let Cx2 be the event that the two agents have
the same x position after round r + m but do not have identical displacements. This requires that the agents
‘wrap’ around the torus, ending at the same position despite moving different amounts in the x direction.
We have P[Cx |Mx = mx ] = P[Cx1 |Mx = mx ] + P[Cx2 |Mx = mx ].
P[Cx1 |Mx = mx ] is identical to the probability that a single random walk takes mx steps and has 0 overall
displacement – i.e. takes an equal number of clockwise and counterclockwise steps. It is given by:
 mx
  mx

1
mx !
1
mx
= mx 2 ·
.
(2)
P[Cx1 |Mx = mx ] =
2
( 2 !)
2
mx /2
1
Above
mx is even
√ we assume

 – otherwise Cx cannot occur. By Stirling’s approximation for any n > 0,
1
n n
n! = 2πn e
1 + O n . Plugging this into 2:
 
mx 
√
  mx


 mx
1 + O m1x
2πmx mex
1
1
m
!
1
x
1
.
=
=Θ √
P[Cx |Mx = mx ] = mx 2 ·

mx 
2 ·

( 2 !)
2
2
mx + 1
1 + O mx1/2
πmx mxe/2

(We use
of mx in the denominator so that the bound holds in the case when mx = 0.)

 mx + 1 instead
P Cx2 |Mx = mx is the probability that two agents have the same x position after round r + m but have
different total displacements.
It is identical to the probability that a single mx step random walk has overall
√
displacement ±c A for some integer c ≥ 1 (and so ‘wraps around’ the torus, ending at its starting location).
Roughly, we bound the probability of this event by the probability that the random walk ends
 at any other
√
location on the torus. There are A such locations, so the probability is bounded by O √1A . We have:
j

m

k

√x

  mx X
A 

 2
mx
1
√
·
P Cx |Mx = mx = 2 ·
mx −c A .
2
2
c=1

(3)

where the extra factor of 2 comes
from the fact that the displacement may be either clockwise or counter√
A
clockwise. (Note that if mx −c
is not an integer we just define the binomial coefficient to equal 0.)
2
5

√
For i ∈ [1, ..., A − 1], let Dxi be the event that a single random walk is i steps clockwise from its starting
location after taking Mx steps. We have:
 mx
1
i
·
P[Dx |Mx = mx ] =
2

j

mx −i
√
A

X

c=−

j

√
For any i ∈ [1, ..., A − 1], and any c ≥ 1,


k

mx +i
√
A

k



mx √

mx +i+c A
2



 mx
1
·
≥
2

−1
X

c=−

j

j

k

m

mx +i
√
A

k



mx √

mx +i+c A
2



√x
 mx X

A 
mx √
1
·
≥
mx +i−c A .
2
2
c=1

√
mx +i−c A
2

is closer to m2x than
 

mx√
mx √
>
mx −c A
mx +i−c A
2

√
mx −c A
2

(4)

is, so
(5)

2

√


mx +i−c A
is an integer. This allows us to lower bound P[Dxi |Mx = mx ] using P Cx2 |Mx = mx .
2
√
A
Ei,c equal 1 if mx +i−c
is an integer and 0 otherwise. Since Cx2 and each Dxi are disjoint events:
2
√
A−1
X



 2
P Dxi |Mx = mx ≤ 1
P Cx |Mx = mx +
i=1
j mx k

√
√

 mx X
A 
A−1


mx √ 
1
X
(applying (4))
P Cx2 |Mx = mx +
·

A  ≤1
m
+i−c
x
2
2
c=1
i=1

as long as
Let

j
k


m
√x
 mx X

 √X
A
A−1

 2
mx√
1

·
P Cx |Mx = mx +
·
Ei,c  ≤ 1
mx −c A
2
2
c=1
i=1
√

 2
P Cx |Mx = mx · Θ( A) ≤ 1.

(by (5) and switching summations)

√ 
√
P√A−1
A
A for all c since mx +i−c
The last step follows from combining (3) with the fact that i=1 Ei,c = Θ
2


√


is integral for half the possible i ∈ [1, ..., A − 1]. Rearranging, we have P Cx2 |Mx = mx = O √1A .


 
Combining our bounds for Cx1 and Cx2 , P [Cx |Mx = mx ] = Θ √m1 +1 + O √1A . Identical bounds hold
x

for the y direction and by (1) we have:
P [C|Mx = mx , My = my ] = Θ

1
p
(mx + 1)(my + 1)

!

+O

1

1
p
+p
A(my + 1)
A(mx + 1)

!

+O

 
1
.
A

(6)

Our final step is to remove the conditioning on Mx and My . Since direction is chosen independently and
uniformly at random for each step, E Mx = E My = m. By a standard Chernoff bound:


2
1
.
P[Mx ≤ m/2] ≤ 2e−(1/2) ·m/2 = O
m+1
(Again using m + 1instead
 of m to cover the m = 0 case). An identical bound holds for My , and so, except
with probability O

1
m+1

P [C] = Θ



both are ≥ m/2. Plugging into (6) this gives us:

1
m+1



+O

1

p
A(m + 1)

!

6

 


 
1
1
1
+O
=Θ
+O
.
A
m+1
A

We note that the techniques of Lemma 4 also apply to bounding the probability that a single random
walk returns to its origin (equalizes) after m steps.
Corollary 5 (Equalization
√
√ Probability Bound). Consider agent a1 randomly walking on a two-dimensional
torus of dimensions A× A. If a1 is located at position
after round r, for any even m ≥ 0, the probability
 p

1
1
+O A
.
that a1 is again at position p after round r + m is Θ m+1

Proof. The analysis of Lemma 4 treats the two walks of a1 and a2 as a single walk with 2m total steps. An
identical analysis where 2m is replaced by m yields the corollary.

Roughly, assuming as in Theorem

 2 that t ≤ A, by Lemma 4, in t rounds, a expects to re-collide with any
Pt−1
1
agent it encounters m=0 Θ m+1 = Θ(log t) times. By Lemma 3, a expects to be involved in dt = nt/A


total collisions. So accounting for re-collisions, it expects to collide with Θ A nt
log t unique individuals. This
is formalized in Lemma 6.


t
Lemma 6 (First Collision Probability). Assuming t ≤ A, for all j ∈ [1, ..., n], P [cj ≥ 1] = Θ A log
t .
Proof. Using the fact that cj is identically distributed for all j,

n
X
n
n
1
ci = · E cj = · P [cj ≥ 1] · E[cj |cj ≥ 1]
E d˜ = d = · E
t
t
t
i=1
n
n
= · P [cj ≥ 1] · E[cj |cj ≥ 1].
A
t

Rearranging gives:
P [cj ≥ 1] =

t
.
A · E[cj |cj ≥ 1]

(7)


4 and linearity of expectation. Since t ≤ A, the O A1 term in
To compute E[cj |cj ≥ 1], we use

 Lemma
1
Lemma 4 is absorbed into the Θ m+1
. Let r ≤ t be the first round that the two agents collide. We have:
E[cj |cj ≥ 1] =

t−r
X

m=0

Θ



1
m+1



= Θ (log(t − r)) .

(8)

After any round the agents are located at uniformly and independently chosen positions, so collide with
probability exactly 1/A. So, the probability of the first collision between the agents being in a given round
can only decrease as the round number increases. So, at least 1/2 of the time that cj ≥ 1, there is a
collision in thefirst t/2
 rounds. So, overall, by (8), E[cj |cj ≥ 1] = Θ (log(t − t/2)) = Θ (log t). Using (7),

P [cj ≥ 1] = Θ

t
A·log t

, completing the proof.

We now give our main technical lemma – a strong moment bound on the distribution of cj . Intuitively,
not only does an agent expect to collide at most O(log t) times with any other agent it encounters, but this
bound extends to the higher moments of the collision distribution, and so holds with high probability. In
this sense, the grid has strong local mixing – random walks spread quickly over a local area and do not cover
the same nodes too many times.
def

Lemma 7. (Collision
Moment Bound)
For all j ∈ [1, ..., n], let c̄j = cj − E cj . For all k ≥ 2, assuming


 k
k−1
t
t ≤ A, E c̄j = O A · k! log
t .

Proof. We expand E[c̄kj ] = P[cj ≥ 1] · E[c̄kj |cj ≥ 1] + P[cj = 0] · E[c̄kj |cj = 0], and so by Lemma 6:


 




t
E c̄kj = O
· E c̄kj |cj ≥ 1 + E c̄kj |cj = 0 .
A log t
7



k
E c̄kj |cj = 0 = (E cj ) = (t/A)k ≤

k−1
t
A k! log





t for all k ≥ 2. Further, E c̄kj |cj ≥ 1 ≤ E ckj |cj ≥ 1 , since




E cj = At ≤ 1. So to prove the lemma, it just remains to show that E ckj |cj ≥ 1 = O k! logk t .
Conditioning on cj ≥ 1, we know the agents have an initial collision in some round t′ ≤ t. We split cj
P
P ′ +t−1
over rounds: cj = tr=t′ cj (r) ≤ tr=t
cj (r). To simplify notation we relabel round t′ round 1 and so
′
round t′ + t − 1 becomes round t. Expanding ckj out fully using the summation:
" t
#
t
t
t
t
t
X X
X
X
X
X
 k
E cj = E
...
cj (r1 )cj (r2 )...cj (rk ) =
...
E [cj (r1 )cj (r2 )...cj (rk )] .
r1 =1 r2 =1

rk =1

r1 =1 r2 =1

rk =1

E [cr1 (j)cr2 (j)...crk (j)] is just the probability that the two agents
 collide in each of rounds r1 , r2 ,...rk . Assume

w.l.o.g. that r1 ≤ r2 ≤ ... ≤ rk . By Lemma 4 this is: O r1 (r2 −r1 +1)(r3 −r21+1)...(rk −rk−1 +1) . So we can
rewrite, by linearity of expectation:


t
t
t
X
X
X
 k
1
.
O
...
E cj = k!
r1 (r2 − r1 + 1)(r3 − r2 + 1)...(rk − rk−1 + 1)
r =r
r =1 r =r
1

2

k

1

k−1

The k! comes from the fact that in this sum we only have ordered k-tuples and so need to multiple by k! to
account for the fact that the original sum is over unordered k-tuples. We can bound:
t
X

rk =rk−1

1
1
1
= 1 + + ... + = O(log t)
rk − rk−1 + 1
2
t

so rearranging the sum and simplifying gives:
t
X
 k
1
E cj = k!
r
r =1 1
1

t
X

r2 =r1

t
X
1
...
r − r1 r =r
+1 2
k

k−1

1
r
−
rk−1
+1 k

t
t
t
X
X
1 X
1
= k!
...
r
r2 − r1 + 1 r =r
r =1 1 r =r
1

2

k−1

1

k−2

1
· O(log t).
rk−2 − rk−1 + 1

P
We repeat this simplification for each level of summation replacing tri =ri−1
 k
ating through the k levels gives E cj = O(k! logk t) giving the lemma.

1
ri −ri−1 +1

with O(log t). Iter-

As with Lemma 4, the techniques used in Lemmas 6 and 7 can be applied to a single walk. We give two
bounds that may be of independent interest.

Corollary 8√(Random
Walk Visits Moment Bound). Consider an agent a1 randomly walking on a two√
dimensional A × A torus that is initially located at a uniformly random location and takes t ≤ A steps.
Let cj be the number of times that a1 visits node j. For all j ∈ [1, ...A] and all k ≥ 2,


 
t
E c̄kj = O
· k! logk−1 t .
A

Proof. This follows from noting that the expected number of visits to a given
node
 is t/A and so Lemma 6

t
can be used in conjunction with Corollary 5 to show that P[cj ≥ 1] = Θ A log t . We can then just follow
the proof of Lemma 7, using Corollary 5 where needed to obtain the result.

Corollary
√ 9 (Equalization Moment Bound). Consider an agent a1 randomly walking on a two-dimensional
√
A × A torus. If a1 takes t ≤ A steps and c is the number of times it returns to its starting position (the
number of equalizations), for all k ≥ 2,


 
E c̄k = O k! logk t .

Proof. This follows directly from the proof of the moment bound given in Lemma 7 for the number of
collisions between two agents that are assumed to collide at least once: E[ckj |cj ≥ 1] = O(k! logk t). By
replacing the application of Lemma 4 with Corollary 5 we obtain the result.
8

4.3

Concentration of Encounter Rate-Based Density Estimate

P
Armed with the moment bound of Lemma 7 we can finally show that nj=1 c̄j concentrates strongly about
Pn
Pn
its expectation. Since j=1 c̄j is just a mean-centered and scaled version of d˜ = 1t j=1 cj , this is enough
P
to prove the accuracy of encounter rate-based density estimation. We start by showing that nj=1 c̄j is a
sub-exponential random variable.
Pn
Pn
Corollary 10 ( j=1 c̄j is sub-exponential). Assuming t ≤ A, j=1 c̄j is sub-exponential with parameters
i
h Pn
σ 2 λ2
b = Θ(log t) and σ 2 = Θ(td log t). Specifically, for any λ with |λ| < 1b E eλ j=1 c̄j ≤ e 2 .

t
Proof. By Lemma 7, for σ 2 = Θ( t log
A ) and b = Θ(log t), c̄j satisfies the Bernstein condition:

  1
E c̄kj ≤ k!σ 2 bk−2 .
2

t
This implies that c̄j is sub-exponential with parameters σ 2 = Θ( t log
A ) and b = Θ(log t) (see [Wai15], Chapter
2). Since each c̄j is independent, this gives us, for all λ with |λ| ≤ 1b :
n
i Y
h Pn
t log t
2
λ2
E eλc̄j ≤ en· 2 ·Θ( A ) = eΘ(λ td log t) .
E eλ j=1 c̄j =
j=1

This completes the proof by the definition of a sub-exponential random variable.
We finally apply a standard sub-exponential tail bound [Wai15] to prove our main result.
Lemma 11 (Sub-exponential tail bound). Suppose that X is sub-exponential with parameters (σ 2 , b). Then,
for any ∆ ≤

σ2
b ,

∆2

P[|X − E X| ≥ ∆] ≤ 2e− 2σ2 .

Pn
Proof of Theorem 2. Since c̄j is just a mean-centered version of cj , j=1 c̄j deviates from its mean exactly
Pn
Pn
the same amount as j=1 cj . Further, d˜ is just equal to 1t j=1 cj , so the probability that it falls within an
Pn
ǫ multiplicative factor of its mean is the same as the probability that j=1 cj falls within an ǫ multiplicative
factor of its mean. By Corollary 10 and Lemma 11:







 2 
n
n
n
n
X
X
X
X
Θ − ǫlogtdt
cj − td ≥ ǫtd ≤ 2e
cj  = P 
cj − E 
cj  ≥ ǫ E 
δ = P
.
j=1

ǫ2 td
log t

5

j=1

= Θ (log(1/δ)) and so ǫ = Θ

q

j=1

j=1

log(1/δ) log t
td



, yielding the theorem.

Extensions to Other Topologies

We now extend our results to a broader set of graph topologies, demonstrating the generality of the local
mixing analysis discussed above. We illustrate divergence between local and global mixing properties, which
can have significant effects on random walk-based algorithms.

5.1

From Repeat Collision Bounds to Estimation Accuracy

Our proofs are largely independent of graph structure, using just a re-collision probability bound (Lemma 4)
and the regularity (uniform node degrees) of the grid, so agents remain uniformly distributed on the nodes
in each round. Hence, extending our results to other regular graphs primarily involves obtaining re-collision
probability bounds for these graphs.
We consider agents on a graph with A nodes that execute analogously to Algorithm 2, stepping to a
random neighbor in each round. Again, we focus on the multi-agent case but similar bounds (resembling
Corollaries 8 and 9) hold for single random walk. We start with a general lemma, giving density estimation
accuracy in terms of re-collision probability.
9

Lemma 12 (Re-collision Probability to Density Estimation Accuracy). Consider a regular graph with A
nodes such that, if two randomly walking agents a1 and a2 collide in round r, for any 0 ≤ m ≤ t, the
probability that they collide again in round r + m is Θ (β(m)) for some non-increasing function β(m). Let
def Pt
B(t) = m=0 β(m). After running for t ≤ A steps,Algorithm 2 returns
d˜ such that, for any δ > 0, with

q
log(1/δ)B(t)
.
probability ≥ 1 − δ, d˜ ∈ [(1 − ǫ)d, (1 + ǫ)d] for ǫ = O
td
Note that in the special case of the grid, by Lemma 4, we can set β(m) = 1/(m + 1) and hence B(t) =
Θ(log t), yielding Theorem 2.
Proof. E d˜ = d (Lemma 3) still holds as the regularity of the graph ensures that agents remain uniformly distributed on the nodes in every round (the stable distribution of any regular graph is the uniform distribution).
Lemma 6 is also analogous except that (8) becomes:
!
t−r
X
E[cj |cj ≥ 1] = Θ
β(m)
m=0

and using the fact that at least 1/2 the time
P that cj ≥1, there is a collision in the first t/2 rounds and that
t/2
β(m) is non-increasing, E[cj |cj ≥ 1] = Θ
m=0 β(m) = Θ (B(t)) . This gives:
P[cj ≥ 1] = Θ



t
A · B(t)



.


Following the moment calculations in Lemma 7, E[ckj |cj ≥ 1] = O k!B(t)k and hence:
E[c̄kj ]

=O




t
k−1
.
· k!B(t)
A



As in Corollary 10, this gives that c̄j is sub-exponential with parameters b = Θ(B(t)) and σ 2 = Θ tB(t)
A


P
so nj=1 c̄k is sub-exponential with parameters b = Θ(B(t)) and σ 2 = Θ n · tB(t)
= Θ (tdB(t)). Plugging
A

into Lemma 11 gives:

ǫ2 td
B(t)

= Θ(log(1/δ)). Rearranging yields the result.

Applying the above bound requires a constant factor approximation to the re-collision probability – the
probability is Θ(β(m)). Sometimes however, it is much easier to give just an upper bound – so the probability
is O(β(m)). In this case a slightly weaker bound holds:
Lemma 13 (Re-collision Probability Upper Bound to Density Estimation Accuracy). Consider a regular
graph with A nodes such that, if two randomly walking agents a1 and a2 collide in round r, for any 0 ≤ m ≤ t,
the probability that they collide again in round r + m is O (β(m)) for some non-increasing function β(m).
def Pt
˜
Let B(t) = m=0 β(m). After running for t ≤ A steps,
q Algorithm 2 returns d such that, for any δ > 0, with
log(1/δ)·B(t)2
probability ≥ 1 − δ d˜ ∈ [(1 − ǫ)d, (1 + ǫ)d] for ǫ = O
.
td
Proof. The proof is identical to that of Lemma 12 except that, we can only show P[cj ≥ 1] = Θ
Therefore, our moment bound becomes:


t
k
k
E[c̄j ] = O
· k!B(t) .
A
2

t
A


.

This gives that c̄kj is sub-exponential with parameters b = Θ(B(t)) and σ 2 = Θ( tB(t)
A ). Following Lemma
12 we therefore have

ǫ2 td
B(t)2

= Θ(log(1/δ)). Rearranging yields the proof.

10

5.2

k-Dimensional Tori

We first consider k-dimensional tori for general k. As k increases, local mixing becomes stronger, fewer
re-collisions occur, and density estimation becomes easier. In fact, for k ≥ 3, although the torus still mixes
slowly, density estimation is as accurate as on the complete graph! We first give the 1-dimensional case:
5.2.1

Density Estimation on the Ring

Lemma 14 (Re-collision Probability Bound – Ring). If two randomly walking agents a1 and a2 are located
on a 1-dimensional torus (a ring) with A nodes, and collide
in round r, for any m ≥ 0, the probability that


1
a1 and a2 collide again in round r + m for k ≥ 1 is Θ √m+1 + O A1 .

Proof. We have already shown this re-collision
of Lemma 4. It is identical to P[Cx |Mx = m]
 in the proof
 bound

1
1
on an A × A grid, which is bounded by Θ √m+1 + O A .




1
1
For m ≤ A, the O A
term. We can estimate the sum of repeat collision
is absorbed into the Θ √m+1
probabilities by using the fact that
√
√
√
√
1
2
2( x − x − 1)
√ ≤√
√
√
√
= √
= 2( x − x − 1).
√
x
x+ x−1
( x + x − 1)( x − x − 1)
So:
t
X

√
√
√
√
√ 
√
√
1
√
≤ 1+2
2 − 1) + ( 3 − 2) + ... + ( t + 1 − t) = 2 t + 1 − 1.
m+1
m=0
Similarly,

√
√
2
1
√
√
≥ √
= 2( x − x − 1)
x−1
x+ x−1

and so:
t
X


√
√
√
√
√
√
√
1
√
2 − 1) + ( 3 − 2) + ... + ( t + 2 − t + 1) = 2 t + 2 − 2.
≥2
m+1
m=0
√
Pt
1
= Θ( t). Plugging into Lemma 12, on a ring, random walk-based density
So, overall m=0 √m+1
q

q
2 

√ 
log(1/δ) t
log(1/δ)
log(1/δ)
√
estimation gives: ǫ = O
= O
rounds are
. Rearranging, t = Θ
td
ǫ2 d
td

necessary to obtain a 1 ± ǫ approximation
√ with probability ≥ 1 − δ. Local mixing on the ring is much worse
than on the torus– we expect to see Θ( t) rather than Θ(log t) repeat collisions with every agent interacted
with. Hence, density estimation is much more difficult, requiring t to be quadratic rather than linear in 1/d.
5.2.2

Density Estimation on Higher Dimensional Tori

We now cover k ≥ 3. While global mixing time is on the order of A2/k and so is slow if k << A, local
mixing is so strong that our accuracy bounds actually match those of independent sampling! Throughout
this section, we assume that k is a small constant and hide it in our asymptotic notation.
Lemma 15 (Re-collision Probability Bound – High-Dimensional Torus). If two randomly walking agents a1
and a2 are located on a k-dimensional torus with A nodes, and collide
in round
r, for any constant k ≥ 3,



1
m ≥ 0, the probability that a1 and a2 collide in round r + m is Θ (m+1)k/2 + O A1 .

Proof. We closely follow the proof of Lemma 4. In total, a1 and a2 take 2m steps: Mi in each dimension for
i ∈ [1, ..., k]. Let Ci be the event that the agents have the same position in the ith dimension in round r + m.
By the analysis of Lemma 4,




1
1
+O
.
P[Ci |Mi = mi ] = Θ √
A1/k
mi + 1
11

So,
 

 





1
1
1
1
P[C|M1 = m1 , ..., Mk = mk ] = Θ √
· .... · Θ √
.
+O
+O
A1/k
A1/k
m1 + 1
mk + 1

(9)

In expectation, Mi = 2m/k. So by a Chernoff bound,
P[Mi ≤ m/k] ≤ 2e

−(1/2)2 ·2m/3k

=O



1
(m + 1)k/2



again assuming k is a small
 constant.
 Union bounding over all k dimensions, we have Mi ≥ m/k for all i
1
except with probability O (m+1)
and hence by (9):
k/2
P[C] = O



1
(m + 1)k/2



"

+ Θ

1
p
m/k + 1

!

+O



1
A1/k

#k

=Θ



1
(m + 1)k/2



 
1
+O
A

giving the lemma (again, asymptotic notation hides multiplicative factors in k since it is constant).


Pt
1
1
< 1+
We can plug the above bound into Lemma 13. For t ≤ A and k ≥ 3,
m=0 (m+1)k/2 + A
√

P∞
log(1/δ)
1
. Rearranging, we require
m=0 (m+1)k/2 = O(1). So we can set B(t) = 1 and have ǫ = O
td


. This matches the performance of independent sampling up to constants.
t = Θ log(1/δ)
ǫ2 d

5.3

Regular Expanders

When a graph does mix well globally, it also mixes well locally. The number of repeat collisions is low
and accurate density estimation is possible. The most obvious example is the complete graph, on which
random-walk based density estimation is equivalent to density estimation via independent sampling. We
generalize this intuition to any regular expander .
Lemma 16 (Re-collision Probability Bound – Regular Expander). Let G be a k-regular expander with A
nodes and adjacency matrix M. Let W = k1 · M be its random walk matrix, with eigenvalues λ1 ≥ λ2 ≥ ... ≥
λA . Let λ = max{|λ2 |, |λA |}. If two randomly walking agents a1 and a2 collide in round r, for any m ≥ 0,
the probability that they collide again in round r + m is at most λm + 2/A.
Proof. Suppose that a1 and a2 collide at node i in round r. The probability they re-collide at round r + m
m
is ||Wm ei ||22 , since for each j, Wi,j
= (Wm ei )j is the probability an agent is at node j after round r + m
given that it is at node i after round r. We bound this norm using the following lemma on how rapidly an
expander random walk converges to its stable distribution:
Lemma 17 ([Lov93]). Let G be a k-regular expander with A nodes, adjacency matrix M, and random walk
matrix W = k1 · M. Let λ1 ≥ λ2 ≥ . . . ≥ λA be the eigenvalues of W and λ = max{|λ2 |, |λA |}. For each
1 ≤ j ≤ n,
1
(Wm · ei )j −
≤ λm .
A
Now we can bound ||Wm ei ||22 by:
||W

m

ei ||22

=

A
X

(W

m

ei )2j

=

A 
X
1
j=1

j=1

12

A

+ χj

2

def

where χj = (Wm · ei )j −
A · (1/A) = 0. Therefore,

1
A

so that χj ∈ [−1/A, λm ] due to Lemma 17. We have

||Wm ei ||22 =

A 
X
1
j=1

A

+ χj

P

j

χj =

P

j (W

m

ei )j −

2

!
 2
A
A
X
1
2χj
1 X 2
2
≤
+
+ χj = +
χ .
A
A
A j=1 j
j=1
m
χ2j is maximized when
χP
the indices
j = λ is maximized. Let S ⊂ [1, A] beP
P j with
Pthe number of possible
m
m
m
χ
=
0.
Therefore,
|S|
·
λ
≤
−
j with χj = λ . Since j χj = 0, we have j∈S λ + j ∈S
j
j ∈A
/ χj ≤
/
1
1 − |S|/A and so |S| ≤ λm +1/A
. Therefore,

P

j

A
X
j=1

χ2j ≤
≤

X

λ2m +

j∈S

X

χ2j

j ∈S
/

2m

A − |S|
λ
≤ λm + 1/A.
+
λm + 1/A
A2

Thus, ||Wm ei ||22 ≤ λm + 2/A, giving the lemma.
Pt
1
+ 2t/A. Assuming t = O(A),
We now apply Lemma 13, with B(t) = m=0 β(m) ≤ 1−λ
s
!
!
r
log(1/δ)(1/(1 − λ) + 2t/A)2
log(1/δ)
=O
.
ǫ=O
td
td(1 − λ)2


2
Rearranging, t = Θ ǫ2log(1/δ)
d(1−λ)2 , matching independent sampling up to a factor of O(1/(1 − λ) ).

5.4

k-Dimensional Hypercube

Finally, we give bounds for a k-dimensional hypercube. Such a graph has A = 2k vertices mapped to the
elements of {±1}k , with an edge between any two vertices that differ by hamming distance 1. The hypercube
is relatively fast mixing. Its adjacency matrix eigenvalues are [−k, −k + 2, −k + 4, ..., k − 4, k − 2, k]. Since it
is bipartite, we can effectively ignore
eigenvalues and apply Lemma 16 with λ = Θ(1 − 2/k) =
 the negative
log(1/δ) log2 (A)
. However, it is possible to remove the dependence on A
Θ(1 − 1/ log A). This yields t = Θ
ǫ2 d
via a more refined analysis – while the global mixing time of the graph increases as A grows, local mixing
becomes stronger!
Lemma 18 (Re-collision Probability Bound – k-Dimensional Hypercube). If two randomly walking agents
a1 and a2 are located on a k-dimensional hypercube with A = 2k vertices and collide
 inround r, for any
m ≥ 0, the probability that a1 and a2 collide in round r + m is at most (7/10)m + O √1A .

Proof. A node of the hypercube can be represented as a k-bit string and each random walk step seen as
choosing one of the bits uniformly at random and flipping it. If a1 and a2 collide, for each of the bit, the total
number of times a1 and a2 chose that bit must be even. The total number of possible ways for re-collision
to occur at round r + m is exactly the number of ways 2m flips can be placed into k buckets, where each
bucket has even number of elements. This quantity is:
X
(2m)!
.
a
!
· . . . · ak !
1
a +...+a =2m
1

(ai

k

mod 2)≡0

This value is equal to the coefficient of x2m in the exponential generating function
k
k
 x

k  
(2m)! X k x(2i−k)
x4
e + e−x
x2
=
e
.
+
+ ...
= (2m)!
(2m)! 1 +
2!
4!
2
2k i=0 i
13

By differentiating 2m times, we find that the coefficient of x2m is:

k  
k  
X
1 X k
k
k
2m
/2 · (2i − k)2m .
(2i − k) =
i
2k i=0 i
i=0
This summation is exactly E[X 2m ], where X is a sum of k i.i.d. random variables each equal to 1 with
probability 1/2 and −1 otherwise. For any c ∈ (0, 1], we can split the expectation:
E[X 2m ] = E[X 2m ||X| ≥ ck] · P[|X| ≥ ck] + E[X 2m ||X| ≤ ck] · P[|X| ≤ ck]
≤ k 2m · P[|X| ≥ ck] + (ck)2m .

To bound the return probability bound, we this count by the the total number of possible paths taken
by a1 and a2 in m steps, k 2m , giving an upper bound of:
P[|X| ≥ ck] + c2m .
By a Hoeffding bound, P[|X| ≥ ck] ≤ 2e−c
So our final probability bound is:

2

k/2

. If we set c =

p
√
√
ln A/k = ln 2 then P[|X| ≥ ck] ≤ 1/ A.

√
1
1
P[|X| ≥ ck] + c2m ≤ √ + ( ln 2)2m < √ + (7/10)m
A
A
yielding the lemma. Note that, by adjusting c, it is possible trade off the terms in the above bound, giving
stronger inverse dependence on A at the expense of slower exponential decay in m.
√
√
Pt
10
+
t/
Plugging
into
Lemma
13,
we
have
B(t)
=
β(m)
≤
A.
If
we
assume
t
=
O(
A), this gives
m=0
3

q


log(1/δ)
log(1/δ)
and so t = Θ
, matching independent sampling.
ǫ=O
td
ǫ2 d

6

Applications

We conclude by discussing algorithmic applications of our ant-inspired density estimation algorithm (Algorithm 2), variations on this algorithm, and the analysis techniques we develop.

6.1

Social Network Size Estimation

Random walk-based density estimation is closely related to work on estimating the size of social networks
and other massive graphs using random walks [KLS11, KBM12, LL12, LW14]. In these applications, one
does not have access to the full graph (so cannot exactly count the nodes), but can simulate random walks
by following links between nodes [MMG+ 07, GKBM09]. One approach is to run a single random walk and
count repeat node visits [LL12, KBM12]. Alternatively, [KLS11] proposes running multiple random walks
and counting their collisions. This can be significantly more efficient since the dominant cost is typically in
link queries to the network. With multiple random walks, this cost can be trivially distributed to multiple
servers simulating walks independently.
Walks are first run for a ‘burn-in period’ so that their locations are distributed approximately by the
stable distribution of the network. The walks are then halted, and the number of collisions in this final
round are counted. The collision count gives an estimate of the walks’ density. Since the number of walks is
known, this yields an estimate for network size.
We show that ant-inspired algorithms can give runtime improvements over this method. After burn-in,
instead of halting the walks, it is possible to run the random walks for multiple rounds, recording encounter
rates as in Algorithm 2. This allows the use of fewer random walks, decreasing total burn-in cost, and leading
to faster runtimes when mixing time is relatively slow, as is common in social network graphs [MYK10].

14

6.1.1

Random Walk-Based Algorithm for Network Size Estimation

Consider an undirected, connected, non-bipartite graph G = (V, E). Let S be the set of vertices of G that
are ‘known’. Initially, S = {v} where v is a seed vertex. We can access G by looking up the neighborhood
Γ(vi ) of any vertex vi ∈ S and adding Γ(vi ) to S.
To compute the number of nodes |V | in the network, we could scan S, looking up the neighbors of each
vertex and adding them to the set. After querying all nodes in S we will have S = V and will know the
network size. However, the number of queries required equals |V |. The goal is to estimate network size using
a significantly more efficient random-walk based approach.
A number of challenges are introduced by this application. While we can simulate many random walks
on G, we can no longer assume these random walks start at randomly chosen nodes, as we do not have the
ability to uniformly sample nodes from the network. Instead, we must allow the random walks to run for
a burn-in phase of length proportional to the mixing time of G. After this phase, the walks are distributed
approximately according to the stable distribution of G.
In general, G is not regular. In the stable distribution, a random walk is located at a vertex with
probability proportional to its degree. Hence, collisions tend to occur more at higher degree vertices. To
correct for this bias, we count a collision at vertex vi with weight 1/ deg(vi ). As we will see, with this
modification, we must adjust our final estimate by the average degree of the graph, which we must estimate.
Our results depend on a natural generalization of re-collision probability. For any i, j ∈ |V |, let p(vi , vj , m)
be the probability of an m step random walk starting at vi ending at vj . Define:
def

β(m) =

maxi,j p(vi , vj , m)
deg(vj )

Intuitively, this is the maximum m step collision probability, weightedP
by degree since higher degree vertices
t
are more likely to be visited in the stable distribution. Let B(t) = m=1 β(m). Note that this weighted
B(t) is upper bounded by the unweighted B(t) used in Lemmas 12 and 13.
For simplicity, we initially ignore burn-in and assume that our walks start distributed exactly by the
def
deg(vi )
i)
and initial
stable distribution of G. A walk starts at vertex vi with probability pi = Pdeg(v
deg(vi ) = 2|E|
i

locations are independent. We also assume knowledge of the average degree deg = 2|E|/|V |. We later
rigorously analyze burn-in and show to estimate deg, completing our analysis.
Algorithm 3 Random Walk-Based Network Size Estimation
input: step count t, average degree deg, n random starting locations [w1 , ..., wn ] distributed according to
the network’s stable distribution
[c1 , ..., cn ] := [0, 0, ..., 0]
for r = 1, ..., t do
∀j, set wj := randomElement(Γ(wj ))
⊲ Γ(wj ) denotes the neighborhood of wj .
count(wj )
⊲ count(wj ) returns number of other walkers currently at wj .
∀j, set cj := cj + deg(wj )
end for P
deg
c
C := ·n(n−1)tj
return Ã = 1/C

Theorem 19. If n2 t = Θ



Ã ∈ [(1 − ǫ)|V |, (1 + ǫ)|V |].

B(t)deg
ǫ2 δ


· |V | , then with probability at least 1 − δ, Algorithm 3 returns

Throughout this section, we work directly with the weighted total collision count C, showing that it is
close to its expectation with high probability and hence giving the accuracy bound for Ã. As in the density
estimation case, we start by showing that C is correct in expectation.
Lemma 20. E C = 1/|V |.

15

Proof. Let cj (r) be the number of collisions, weighted by inverse vertex degree, walk j expects to be involved
i)
in at round r. In each round all walks are at vertex vi with probability pi = deg(v
2|E| , so:
E cj (r) =


|V | 
|V |
X
n−1
deg(vi ) (n − 1) deg(vi )
n−1X
1
deg(vi ) =
=
·
·
.
2
2|E|
2|E|
deg(v
)
4|E|
2|E|
i
i=1
i=1

By linearity of expectation: E cj =

t(n−1)
2|E| ,

E

P

cj =

tn(n−1)
2|E|

and hence, E C =

deg
2|E|

= 1/|V |.

We now need to show concentration of C about its expectation. Let ci,j be the weighted collision count
between walks wi and wj where i 6= j. It is possible to closely follow the moment bound proof of Lemma 7
and show that ci,j is sub-exponential. However, unlike in the case of regular graphs, we will not be able to
claim that the different ci,j ’s are independent. Hence, we will not be able to use the same sub-exponential
tail bounds employed in Section 4.3.
Instead, we bound the second moment (the variance) of each ci,j and obtain our concentration results via
Chebyshev’s inequality. This leads to a linear rather than logarithmic dependence on the failure probability
1/δ. However, we note that we can simply perform log(1/δ) estimates each with failure probability 1/3 and
return the median, which will be correct with probability 1 − δ.
def

Lemma 21. (Degree
Collision Variance Bound) For all i, j ∈ [1, ..., n] with i 6= j, let c̄i,j = ci,j −
 Weighted

 2 
tB(t)
E ci,j . E c̄i,j = O |E| .
Proof. We can split E c̄2i,j over rounds as:


E

t
X
r=1

!2 
t
t−1 X
t
X
X


ci,j (r)  =
E c̄i,j (r)2 + 2
E [c̄i,j (r)c̄i,j (r′ )]
r=1

r=1 r ′ =r+1




|V |
|V |
2
′ 2
X
X
p(vi , vj , r − r ) 
1
1
 deg(vi ) ·

+2
=t
·
·
2
2
2
(2|E|)
deg(v
)
(2|E|)
deg(v
)
deg(vj )
i
i
r=1 r ′ =r+1
i=1
i=1
j=1



|V |
|V |
t−1
X
X
X
t|V |
 deg(vi ) · β(m)

p(vi , vj , m)
+
2t
≤
4|E|2
(2|E|)2
m=1
i=1
j=1
|V | 
X
deg(vi )2



t−1 X
t
X



t−1
X
β(m)
t
tB(t)
≤
.
+ 2t
=O
4|E|
2|E|
|E|
m=1

Lemma 22 (Total Collision Variance Bound). Let C =
Proof.

Pn

j=1 c̄j

=

P




E 

i,j∈[1,...,n],i6=j c̄i,j .

X

i,j∈[1,...,n],i6=j

P
deg j c̄j
n(n−1)t .


 
E C̄ 2 = O n12 t ·

B(t)|E|
|V |2



.

We closely follow the variance calculation in [KLS11]:

2 

c̄i,j   =

X

i,j∈[1,...,n],i6=j




X

i′ ,j ′ ∈[1,...,n],i6=j



c̄i,j · c̄i′ ,j ′ 

 
 
 
 2 
n
n
n
2
=2
(E c̄i,j ) + 2 · 3!
E c̄i,j + 4!
E c̄i,j c̄i,k .
4
2
3

The first term corresponds to the cases when i = i′ and j = j ′ . The second corresponds
to i 6= i′ and j 6= j ′ ,

in which case c̄i,j and c̄i′ ,j ′ are independent and identically distributed. The 4! n4 multiplier is the number
of ways to choose an ordered set
 of four distinct indices. The last term corresponds to all cases when either
i = i′ or j = j ′ . There are 3! n3 ways to choose an ordered set of three distinct indices, multiplied by two to
16

account for the repeated index being in either the first or second position. Using E c̄i,j = 0 and the bound
on E[c̄2i,j ] from Lemma 21:

2 
 2

 
X
n tB(t)
n


c̄i,j   = O
E 
+ 0 + 2 · 3!
E c̄i,j c̄i,k .
(10)
|E|
3
i,j∈[1,...,n],i6=j

When j 6= k, c̄i,j and c̄i,k are independent and identically distributed conditioned on the path that walk wi
traverses. Let Ψi be the t step path chosen by wi .
X
P [Ψi = ψi ] · E [c̄i,j |Ψi = ψi ] · E [c̄i,k |Ψi = ψi ]
E [c̄i,j c̄i,k ] =
ψi

=

X
ψi

=

X
ψi

P [Ψi = ψi ] · E [c̄i,j |Ψi = ψi ]

2

P [Ψi = ψi ] · (E [ci,j |Ψi = ψi ] − E [ci,j ])2 .

(11)

Pt
t
i (r))
E [ci,j |Ψi = ψi ] = r=1 deg(ψ
· deg(ψ1i (r)) = 2|E|
= E [ci,j ]. That is, the expected number of collisions is
2|E|
identical for every path of wi . Plugging into (11), E [c̄i,j c̄i,k ] = 0.

2 
 2
P
tB(t)
and thus:
= O n |E|
So finally, plugging back into equation (10), E
i,j∈[1,...,n],i6=j c̄i,j
h 2i
E C =O

n2 tB(t)
·
|E|



deg
n(n − 1)t

2 !

=O



1 B(t)|E|
·
n2 t
|V |2



.

With this variance bound in place, we can now finally prove Theorem 19.
Proof of Theorem 19. Note that C̄ = C − E C. By Chebyshev’s inequality Lemma 22 gives:
P [|C − E C| ≥ ǫ E C] ≤
Rearranging gives us that, in order to have C̄ ∈
2

h

n t=Θ

1
ǫ 2 n2 t

(1−ǫ) 1+ǫ
|V | , |V |



i

· B(t)|E|.

with probability δ, we must have:


1
B(t)|E| .
ǫ2 δ

As long as we have ǫ < 1/4 then this immediately gives Ã ∈ [(1 − 2ǫ)|V |, (1 + 2ǫ)|V |], giving the lemma after
adjusting constants on ǫ.
6.1.2

Estimating The Average Degree

we need a (1 ± ǫ) apWe now show how to estimate the value of deg used in Algorithm 3. Specifically,
P
cj

j
1
proximation to deg
. If we then substitute this into the formula Ã = deg·n(n−1)t
, we still have a (1 ± O(ǫ))
approximation to the true network size. We use the algorithm and analysis of [KLS11], which gives a simple
approximation via inverse degree sampling.

Algorithm 4 Average Degree Estimation
input: n random starting locations [w1 , ..., wn ] distributed according to the network’s stable distribution.
1
∀j, set dj := deg(w
⊲ Sampling
j)
return D :=

P

dj
·n

17

Theorem 23 (Average Degree Estimation). If n = Θ
h
i
1+ǫ
probability at least 1 − δ, D ∈ 1−ǫ
.
,
deg deg



1
ǫ2 δ

·

deg
degmin


, Algorithm 4 returns D such that, with

Proof.


|V | 
n
X
1
|V |
deg(vi )
1
1
1X
.
E dj = · n ·
=
·
=
ED =
n j=1
n
2|E|
deg(vi )
2|E|
deg
i=1
Since each dj is independent, letting D̄ = D − E D,
|V |
  1
1
1
1
1 X deg(vi )
|V |
1
1
E D̄2 = E[d¯2j ] ≤ E[d2j ] =
.
≤ ·
=
·
n
n
n i=1 (2|E|) deg(vi )2
n 2|E| degmin
n degmin deg
h
i
deg
1
ǫ
Applying Chebyshev’s inequality: P D − deg
≤ deg
and rearranging, to succeed with proba≤ ǫ2 n deg
min
bility at least 1 − δ it suffices to set:


1
deg
.
n=Θ 2 ·
ǫ δ degmin

6.1.3

Handling Burn-In Error

Finally, we remove our assumption that walks start distributed exactly according to the network’s stable
distribution, rigorously bounding the length of burn-in required before running Algorithm 3.
n
Let D∗ ∈ R|V | be a vector representing the true stable distribution of n random walks on G and
n
Dt ∈ R|V | be a vector representing the distribution of the walks after running for t burn-in steps. Specifically,
each walk w1 , ..., wn is initialized at a single seed vertex v. For t rounds we then update the location of each
walk independently by moving to a randomly chosen neighbor. Both vectors are probability distributions:
they have all entries in [0, 1] and kD∗ k1 = kDk1 = 1.
Let ∆ = D∗ − Dt and assume that k∆k1 ≤ δ. We can consider two equivalent algorithms: draw an
initial set of locations W = w1 , ..., wn from D∗ , run Algorithm 3, and then artificially fail with probability
max{0, ∆(W )}. Alternatively, draw W = w1 , ..., wn from Dt , run Algorithm 3, and then artificially fail with
probability max{0, −∆(W )}. These algorithms are clearly equivalent. The first obtains a good estimator
with probability 1 − 2δ - probability δ that Algorithm 3 fails when initialized via the stable distribution D∗
by Theorem 19 plus an artificial failure probability of ≤ k∆k1 ≤ δ. The second then clearly also fails with
probability 2δ. This can only be higher than if the we did not perform the artificial failure after running
Algorithm 3. Therefore, running Algorithm 3 with a set of random walks initially distributed according to
Dt yields success probability ≥ 1 − 2δ.
How long must the burn-in period be to ensure kD∗ − Dt k1 ≤ δ? Let W be the random walk matrix of
G. Let λ1 ≥ λ2 ≥ . . . ≥ λA be the eigenvalues of W and λ = max{|λ2 |, |λ|V | |}. Let Ct ∈ R|V | denote the
location distribution for a single random walk after burn-in and C ∗ ∈ R|V | denote the stable distribution of
a single random walk. If we have, for all i, |Ct (vi ) − C ∗ (vi )| ≤ δ/n · C ∗ (vi ) then we have for any W :
|Dt (W ) − D∗ (W )| =
≤

n
Y

i=1
n
Y

i=1

Ct (wi ) −

n
Y

i=1

C ∗ (wi )

(C ∗ (wi ) + δ/n · C ∗ (wi )) −

n
Y

i=1

C ∗ (wi )

n
n  
X
X
n
δ i ≤ 2δ · D∗ (W )
(δ/n)i ≤ D∗ (W )
< D∗ (W )
i
i=1
i=1

as long as δ < 1/2. This multiplicative bound gives kD∗ − Dt k1 ≤ 2δ.
mixing
time


 bounds
 By standard
log(|E|/δ)
log(n|E|/δ)
δ
∗
∗
([Lov93], Theorem 5.1), |Ct (vi ) − C (vi )| ≤ n|E| · C (vi ) for all i after O
=O
burn-in
1−λ
1−λ

steps (since n < |E| or else we could have scanned the full graph.)
18

6.1.4

Overall Runtime and Comparision to Previous Work


Let M = O log(|E|/δ)
denote the burn-in time required before running Algorithm 3. In order to obtain a
1−λ
(1 ± ǫ) estimate of network size with probability 1 − δ we must run n random walks for M + t steps, making
n(M + t) link queries, where by Theorems 19 and 23:



s

deg
|V | · B(t)deg 
n = Θ max
,
.
(12)
2
 degmin ǫ δ

t · ǫ2 δ

Typically, the second term dominates since deg << |V |. Hence, by increasing t, we are able to use fewer
random walks, significantly decreasing the number of link queries required if M is large.
In the special case with t = 1 we obtain a somewhat simpler bound. Instead of using the more general
analysis of Lemma 21 we can directly calculate:
|V |

  X deg(vi )2
1
1
.
=
E c̄2i,j =
2 deg(v )2
(2|E|)
deg
·
2|E|
i
i=1
2

Plugging this into Lemma 22 gives E[C ] =

n2
deg·2|E|

·



deg
n(n−1)

2

= O



|E|



. So, to have, C ∈ [(1 −


 
deg
ǫ) E C, (1 + ǫ) E C] with probability 1 − δ by Chebyshev’s inequality, we need: n = Θ |E|
ǫ2 δ = Θ ǫ2 δ · |V | .
So overall require:



s

deg
|V | · deg 
n = Θ max
,
.
(13)
2
 degmin ǫ δ
ǫ2 δ 
n2 |V |2
2

[KLS11] also halts random walks after burn-in (uses t = 1), but uses a different estimator, tracking
average degree, average inverse degree, and unweighted collisions. They require:
)!
(s
P 3
P
pi
1/pi
1
P ,
.
(14)
n = Θ max
,
P
ǫ2 δ p2i ǫ2 δ ( p2i )2 ǫ2 δ|V |2

i)
Where pi = deg(v
2|E| . Their third term comes from deg estimation and can be upper bounded by our first,
although, by giving a tighter analysis in Theorem 23 we match this term. Their first term can be rewritten:
s
s
s
s
1
2|E|
2|E|
|V | · deg
= P|V |
.
·
P|V | 2 =
P|V |
2
2
2
2
ǫ2 δ
ǫ δ
p
ǫ δ
deg(vi ) /(2|E|)
deg(vi )

i=1

i

i=1

i=1

P|V |
This will always be somewhat smaller than our second term as 2|E| ≤ i=1 deg(vi )2 . Their second term is
a bit harder to compare, however, can be easily upper bounded by:
P 3
deg3max 2|E|
pi
1
deg3max
≤
· 2 =
P
2
4
3 · ǫ2 δ|V | .
2
2
ǫ
δ
2
ǫ δ ( pi )
deg |V |
deg

Assuming degmax /deg is not too large, this term will be small. However, a few very high degree nodes
in an otherwise sparse graph can make it very large.
In sum, while our bounds are not directly comparable to those of [KLS11], in the t = 1 case, assuming
reasonably bounded degrees, they are of the same order of magnitude. Further, 12 gives an important
tradeoff for graphs with relatively slow mixing time – we can increase the number of steps in our random
walks, decreasing the total number of walks.
As an illustrative example, consider a k-dimensional torus graph for k ≥ 3 (for k = 2 mixing time is Θ(|V |)
so we might as well census the full graph). For any δ > 0, the mixing time required for kD∗ − DM k1 = O(δ)
19

(see Section 6.1.3) is M = Θ(log(|V |/δ)|V |2/k ). All nodes have degree 2k, and using the bounds above, to
obtain a (1 ± ǫ) estimate of |V |, the algorithm of [KLS11] requires:
(
)!
r
1
|V |
n = Θ max 2 ,
,
ǫ δ
ǫ2 δ


|/δ)
(k+4)/(2k)
√
Assuming |V | is large so the second term dominates, they require M · n = Θ log(|V
·
|V
|
ǫ d
link queries to obtain a size estimate. In contrast, we require:



s
)!
(
r


deg
|V
|
·
B(t)deg
|V |
1


.
n = Θ max
,
= Θ max 2 ,

 degmin ǫ2 δ
t · ǫ2 δ
ǫ δ
t · ǫ2 δ

since by Lemma 15, B(t) = O(1/k) and deg = degmin = k. If we set t = Θ(M ), and again assume that the
second term dominates since |V | is large, the total number of link queries we need is:
!
p
log(|V |/δ)
(k+1)/2k
√
· |V |
.
n(M + t) = O
ǫ d
This beats [KLS11] by decreasing the polynomial in |V | and the logarithmic burn-in term. Ignoring error
dependences, if k = 3, [KLS11] requires Θ(n7/6 ) queries which is more expensive than fully censusing the
graph, whereas we require O(n2/3 ) queries, which is sublinear in the graph size.
We leave it as an open question to compare our bounds with those of [KLS11] on more natural classes of
graphs, and to determine either experimentally or theoretically, typical values of B(t) on these graphs. As
used above, any of the bounds obtained in Section 5 apply since the degree weighted B(t) is upper bounded
by the unweighted B(t) used in our regular graph analysis. However, bounds for real networks or popular
random graph models used to study these networks would be very interesting.

6.2

Distributed Density Estimation by Robot Swarms

Algorithm 2 can be directly applied as a simple and robust density estimation algorithm for robot swarms.
Additionally, the algorithm can be used to estimate the frequency of certain properties within the swarm. Let
d be the overall population density and dP be the density of agents with some property P . Let fP = dP /d
be the relative frequency of P .
Assuming that agents with property P are distributed uniformly in population and that agents can
detect this property (through direct communication or some other signal), then they can separately track
encounters with these
agents. They can compute
an estimate d˜ of d and d˜P of dPh. ByTheorem
i


2, after

log(1/δ) log log(1/δ) log(1/dǫ)
1−ǫ
1+ǫ
˜
˜
steps,
with
probability
1
−
2δ,
d
/
d
∈
f
,
running for t = Θ
P
P
d P ǫ2
1+ǫ
1−ǫ fP =

[(1 − O(ǫ))fP , (1 + O(ǫ))fP ] for small ǫ.
In a biological setting, properties may include if an ant has recently completed a successful food foraging
trip [Gor99], or if an ant is a nestmate or enemy [Ada90]. In a robotics setting, properties may include
whether a robot is part of a certain task group, whether a robot has completed a certain task, or whether a
robot has detected a certain event or environmental property.

6.3

Random Walk-Based Sensor Network Sampling

Finally, we believe our moment bounds for a single random walk (Corollaries 8 and 9) can be applied to
random walk-based distributed algorithms for sensor network sampling. We leave obtaining rigorous bounds
in this domain to future work.
Random walk-based sensor network sampling [LB07, AB04] is a technique in which a query message (a
‘token’) is initially sent by a base station to some sensor. The token is relayed randomly between sensors,
which are connected via a grid network, and its value is updated appropriately at each step to give an answer
to the query. This scheme is robust and efficient - it easily adapts to node failures and does not require
setting up or storing spanning tree communication structures.
20

However, if attempting to estimate some quantity, such as the percentage of sensors that have recorded
a specific condition, as in density estimation, unless an effort is made to record which sensors have been previously visited, additional variance is added due to repeat sensor visits. Recording previous visits introduces
computational burden – either the token message size must increase or nodes themselves must remember
which tokens they have seen. We are hopeful that our moment bounds can be used to show that this is
unnecessary – due to strong local mixing, the number of repeat sensor visits will be low, and increased
variance due to random walking will be limited.

6.4

Generalization to Data Aggregation

We note that, estimating the percentage of sensors in a network or the density of robots in a swarm with
a property that is uniformly distributed is a special case of a more general data aggregation problem: each
agent or sensor holds a value vi drawn independently from some distribution D. The goal is to estimate
some statistic of D - e.g. its expectation. In the case of density estimation vi is simply an indicator random
variable which is 1 with probability d and 0 otherwise.
Ideally, to estimate the expectation of D, one would independently sample the agents and compute
the sample mean. As in the case of density estimation, if instead, one simply takes an average of values
encountered via a random walk, additional variance is added due to repeated encounters with many agents.
However, it is likely that the effect of this variance can be bounded, as we have done in the special case of
density estimation. We believe extending our results to more general data aggregation problems via random
walk is an interesting future direction.

7

Acknowledgements

We thank Amartya Shankha Biswas, Christopher Musco, and Mira Radeva for helpful discussions. This
work was partially supported by NSF Grants BIO-1455983 and CCF-1461559, the NSF Center for Science
of Information grant CCF-0939370, and AFOSR grant FA9550-13-1-0042.

References
[AB04]

Chen Avin and Carlos Brito. Efficient and robust query processing in dynamic environments using random walk techniques. In Proceedings of the 3rd International Symposium on Information
Processing in Sensor Networks, pages 277–286. ACM, 2004.

[Ada90]

Eldridge S Adams. Boundary disputes in the territorial ant Azteca trigona: effects of asymmetries
in colony size. Animal Behaviour, 39(2):321–328, 1990.

[CLLM12] Kai-Min Chung, Henry Lam, Zhenming Liu, and Michael Mitzenmacher. Chernoff-hoeffding
bounds for markov chains: Generalized and simplified. arXiv preprint arXiv:1201.0559, 2012.
[ES09]

Robert Elsässer and Thomas Sauerwald. Tight bounds for the cover time of multiple random
walks. In Automata, Languages and Programming, pages 415–426. Springer, 2009.

[Gil98]

David Gillman. A chernoff bound for random walks on expander graphs. SIAM Journal on
Computing, 27(4):1203–1220, 1998.

[GKBM09] Minas Gjoka, Maciej Kurant, Carter T Butts, and Athina Markopoulou. A walk in facebook:
Uniform sampling of users in online social networks. arXiv preprint arXiv:0906.0060, 2009.
[Gor99]

Deborah M Gordon. Interaction patterns and task allocation in ant colonies. In Information
Processing in Social Insects, pages 51–67. Springer, 1999.

[GPT93]

Deborah M Gordon, Richard E Paul, and Karen Thorpe. What is the function of encounter
patterns in ant colonies? Animal Behaviour, 45(6):1083–1100, 1993.

21

[KBM12]

Maciej Kurant, Carter T Butts, and Athina Markopoulou. Graph size estimation. arXiv preprint
arXiv:1210.0460, 2012.

[KLS11]

Liran Katzir, Edo Liberty, and Oren Somekh. Estimating sizes of social networks via biased
sampling. In Proceedings of the 20th International Conference on World Wide Web, pages 597–
606. ACM, 2011.

[LB07]

Luisa Lima and Joao Barros. Random walks on sensor networks. In 5th International Symposium
on Modeling and Optimization in Mobile, Ad Hoc and Wireless Networks and Workshops, 2007,
pages 1–5. IEEE, 2007.

[LL12]

Jianguo Lu and Dingding Li. Sampling online social networks by random walk. In Proceedings
of the First ACM International Workshop on Hot Topics on Interdisciplinary Social Networks
Research, pages 33–40. ACM, 2012.

[Lov93]

László Lovász. Random walks on graphs: A survey. Combinatorics, Paul Erdos is Eighty,
2(1):1–46, 1993.

[LW14]

Jianguo Lu and Hao Wang. Variance reduction in large graph sampling. Information Processing
& Management, 50(3):476–491, 2014.

[MMG+ 07] Alan Mislove, Massimiliano Marcon, Krishna P Gummadi, Peter Druschel, and Bobby Bhattacharjee. Measurement and analysis of online social networks. In Proceedings of the 7th ACM
SIGCOMM Conference on Internet Measurement, pages 29–42. ACM, 2007.
[MYK10]

Abedelaziz Mohaisen, Aaram Yun, and Yongdae Kim. Measuring the mixing time of social
graphs. In Proceedings of the 10th ACM SIGCOMM Conference on Internet Measurement,
pages 383–389. ACM, 2010.

[NTD05]

Stamatios C Nicolis, Guy Theraulaz, and Jean-Louis Deneubourg. The effect of aggregates on
interaction rate in ant colonies. Animal Behaviour, 69(3):535–540, 2005.

[Pra05]

Stephen C Pratt. Quorum sensing by encounter rates in the ant Temnothorax albipennis. Behavioral Ecology, 16(2):488–496, 2005.

[SHG06]

Robert J Schafer, Susan Holmes, and Deborah M Gordon. Forager activation and food availability in harvester ants. Animal Behaviour, 71(4):815–822, 2006.

[Wai15]

Martin J Wainwright. High-dimensional statistics: A non-asymptotic viewpoint, draft.
http://www.stat.berkeley.edu/~mjwain/stat210b/Chap2_TailBounds_Jan22_2015.pdf,
2015.

22

