arXiv:1405.2322v2 [] 4 Oct 2016

Parameter estimation of a two-colored urn model
class
Line Le Goff

Philippe Soulier

January 29, 2018
Abstract
Though widely used in applications, reinforced random walk on graphs
have never been the subject of a valid statistical inference. We develop in
this paper a statistical framework for a general two-colored urn model. The
probability to draw a ball at each step depends on the number of balls of
each color and on a multidimensional parameter θ through a function f ,
called a choice function. We introduce two estimators of θ: the maximum
likelihood estimator and a weighted least squares estimator which is less
efficient, but is closer to the calibration techniques used in the applied literature. In general, the model is an inhomogeneous Markov chain and
this property makes the estimation of the parameter impossible on a single path, even if it were infinite. Therefore we assume that we observe
i.i.d. experiments, each of a predetermined finite length. This is coherent
with the usual experimental set-ups. We apply the statistical framework to
a real life experiment: the selection of a path among pre-existing channels
by an ant colony. We performed experiments, which consisted of letting
ants pass through the branches of a fork. We consider the particular urn
model proposed by J.-L. Deneubourg in 1990 to describe this phenomenon.
We simulate this model for several parameter values in order to assess the
accuracy of the MLE and the WLSE. Then we estimate the parameter
from the experimental data and evaluate confident regions with Bootstrap
algorithms. The findings of this paper do not contradict the biological literature, but give statistical significance to the values of the parameter found
therein.

1

1

Introduction

Urn models have been studied for nearly one century. In 1931, G. Pólya provided
the first probabilistic result on the game consisting in drawing a ball from an urn
initially containing one red ball and one black ball (see Pólya, 1931). At each time
step, a ball is drawn and put back in the urn with an additional ball of the same
color. The probability to draw a red ball is the proportion of red balls in the urn.
G. Pólya proved that, as the number of draws tends to infinity, the proportion of
red balls tends to a random variable following the uniform distribution on [0, 1].
The Pólya urn is easily generalizable to a large class of two-colored urn models
characterized by a choice function f : Θ × N × N → [0, 1] which itself depends on
a parameter θ ∈ Θ ⊂ Rd , d ≤ 1. Let Rn and Bn be the numbers of red and black
balls in the urn after n draws. Note that, by construction, Rn + Bn = n. The
probability that the (n + 1)-th ball is red is given by
probR
n+1 = f (θ, Rn , Bn ) = f (θ, Rn , n − Rn ) .

(1)

Consequently, the probability to draw a black ball at time n + 1 is 1 − probR
n+1 .
The goal of this paper is to propose a valid statistical methodology to estimate
the parameter θ. We suppose that we observe N independent paths, each consisting of a sequence of n colors drawn by the model defined in (1). The statistical
theory is developed as n is fixed and N tends to infinity. We choose this framework since in some models, it is not possible to obtain consistent estimators of the
parameter with only one path, even if its length n increases to infinity. Moreover
data from real experiments comprises a set of finite paths. We define estimators
for θ that we prove to be consistent and asymptotically normal under some usual
regularity assumptions on the model (1). We study more precisely two particular
cases: the maximum likelihood estimator (MLE) and the weighted least squares
estimators (WLSE).
We have applied these statistical tools to the problem of path formation by
an ant colony. One of the fundamental factors affecting an organism’s survival is
its ability to optimally and dynamically exploit its environment. For example, in
order to take advantage of the best sites of resources, housing or reproduction,
these areas must be discovered and exploited at the earliest opportunity. Many
species of ants rise to this challenge by developing a network of paths, which
connects different strategic sites such as nests and food sources. These paths
consist of pheromones, attractive chemical substances. We focus on a specific
aspect of this phenomenon: the selection of a path among pre-existing channels.
When exploring their environment, ants often face bifurcations and must bypass
obstacles. As shown by experimental studies, the laying of pheromones by ants
passing successively through a bifurcation results in two possible outcomes: either
one branch is eventually selected and the other abandoned, or both branches end
2

up being uniformly chosen (see Deneubourg, Aron, Goss, and Pasteels, 1990).
The analysis of the spontaneous path formation by a colony of ants is made
difficult by the absence of any means to measure precisely the quantity of, or even
detect, the pheromones laid by the ants.
It is commonly assumed that when approaching a bifurcation, ants choose a
branch and lay a certain constant amount of pheromone without ever turning back.
Consequently, the quantity of pheromone laid on a branch is proportional to the
number of ants which passed through it. Thus this phenomenon can be described
by a urn model as proposed by J.-L. Deneubourg et al. in 1990 (see Deneubourg
et al., 1990). More precisely they define a choice function with a two-dimension
parameter (α, c) ∈ (0, 1)2 such that the probability for an ant to choose the right
branch after Rn passages through the right branch and n passages in total is given
by
(c + Rn )α
.
(2)
P robR
n+1 =
(c + Rn )α + (c + n − Rn )α
The probability to choose the left branch is consequently 1 − P robR
n+1 . The parameter α makes this model non linear with respect to the proportion of passages
through one branch. It models the sensitivity of the ant to the concentration of
pheromone. The parameter c is the intrinsic attractiveness of each branch and
can also be interpreted as the inverse of the attractiveness (or the strength) of the
pheromone deposit laid by each ant.
Several probabilistic studies provide the asymptotic behavior of Rn /n in terms
of α and c (see Pólya, 1931; Tarrès, 2011; Davis, 1990). The influence of the two
parameters is on different time scales, but they can contribute to the same effect
(selection of one branch or unifomization of the traffic on the two branches) or have
antagonistic effects. The model is thus characterized by four phases, according to
the values of α and c: slow or fast uniformization; slow or fast selection. The phase
most commonly considered in the literature is slow selection. This corresponds
in our model to α and c larger than 1: selection of one branch will eventually
happen, though slowly because of weak pheromone deposits (see Deneubourg
et al., 1990; Vittori, Talbot, Gautrais, Fourcassié, Araujo, and Theraulaz, 2006;
Garnier, Guérécheau, Combe, Fourcassié, and Theraulaz, 2009). However, the
model may account for other possibilities, such as fast uniformization, which occurs when α < 1 and c > 1. One purpose of this paper is to investigate more
thoroughly these possibilities which have been more or less overlooked in the previous literature.
Ethological studies have already provided values for the parameter (α, c), but
the methods used mainly consisted of calibration without control of the statistical
validity of these methods and results (see Deneubourg et al., 1990; Vittori et al.,
2006; Garnier et al., 2009; Thienen, Metzler, Choe, and Witte, 2014) and in particular these methods do not produce confidence regions. However this type of
3

information supplies interesting and important elements to the behavioral discussion. However this type of information supply interesting and important elements
to the behavioral discussion. In this paper, we define the MLE and the WLSE for
the ant behavior model (2). We assess the quality of these estimators in a simulation experiment. We then use them on experimental data provided by a real
life experiment performed to this purpose with ants. We also compute confidence
region by a Bootstrap algorithm.
The model (2) can be applied, mutatis mutandis, to many other fields. For
instance we can consider a fork with two branches as a neuron having two axons.
During each period of time, the length of one of the two axons increases. The
longer an axon is, the higher the probability that it will further grow. Thus
the dynamic of this biological system may also be modeled by (2) (see Khanin
and Khanin, 2001). Furthermore, using the notion of a choice function (see
Section 2.3), the model (1) can be adapted and applied to many situations where
a binary choice occurs, or to even more complex situations such as networks with
several nodes (see Jeanson, Ratnieks, and Deneubourg, 2003; Pemantle, 2007), for
instance network exploration by an ant colony (see Aron, Deneubourg, Goss, and
Pasteels, 1990; Beckers, Deneubourg, and Goss, 1993; Nicolis and Deneubourg,
1999; Dussutour, Deneubourg, and Fourcassié, 2005; Nicolis and Dussutour, 2008;
Thienen et al., 2014; Arganda, Nicolis, Perochain, Péchabadens, Latil, and Dussutour, 2014).
The paper is organized as follows. The model (1) and the statistical framework is rigorously defined in Section 2. Section 2.1 is focused on the MLE and
Section 2.2 on the WLSE. Under some usual regularity conditions, we prove that
both estimators are consistent and asymptotically normal. Moreover the MLE
is asymptotically efficient. The numerical implementation of the MLE may be
difficult and unstable (and lengthy), therefore, a WLSE is considered and theoretically studied. This estimator does not match the theoretical performances of
the MLE, but is easier to compute and is popular among practitioners. Section 2.3
proposes an extension of the general urn model to a class of vertex reinforced random walk on graphs. Section 3 is an adaptation of the statistical framework to
the ethological problem. In Section 3.1, we first introduce the assumptions on
the ant behavior and then we define the model proposed by J.-L. Deneubourg.
The four phases of the model are described more precisely and are interpreted
from the point of view of ethology. Section 3.2 rewrites the MLE and the WLSE
for this particular case. In Section 3.3, we show that it is not possible to consistently estimate the parameter of the model on a single experiment of length n,
even if n tends to infinity. In order to assess the performance of the estimators,
a short simulation experiment is reported in Section 3.4. Section 4 reports the
study on the experimental data. The experimental protocol and the data produced are described in Sections 4.1 and 4.2. Section 4.3 supplies the estimation
4

results (computation of the estimators and their confidence regions). We provide
some concluding remarks in Section 5 and prove the theoretical statistical results
of this paper in Section 6.

2

Parameter estimation

Let us first write precisely the model studied and statistical framework used. We
assume that {Xk , k ≥ 1} is a sequence of Bernoulli random variables (representing
the colors of the balls drawn: 1 for red and 0 for black) and that there exists a
function f0 : N2 → [0, 1] such that, for all integers k,
P(Xk+1 = 1 | Fk ) = f0 (Zk , k − Zk ) .
where Z0 = 0 and for k ≥ 1, Zk = X1 + · · · + Xk and Fk is the sigma-field
generated by Z0 , X1 , . . . , Xk . The random walk {Zk , k ≥ 0} is an inhomogeneous
Markov chain. For n ≥ 1 and a sequence (e1 , . . . , en ) ∈ {0, 1}n , applying the
Markov property, we obtain
P(X1 = e1 , . . . , Xn = en ) =

n−1
Y

f0 (zk , k − zk )ek+1 {1 − f0 (zk , k − zk )}1−ek+1 , (3)

k=0

where zk = e1 + · · · + ek with z0 = 0, by convention.
We assume that we observe N experiments, each consisting in a path of
length n of the model (3). In other words, we have a set of N sequences of n consecutive draws. For j = 1, . . . , N and k = 1, . . . , n, let Xkj ∈ {0, 1} denote the
color of the k-th ball drawn in the j-th experiment. Let Z0j = 0 and Zkj =
Pk
j
i=1 Xi , k ≥ 1 be the total number of red balls drawn at time k during the j-th
j
experiment, so that Xkj = Zkj − Zk−1
. In all the paper, n will be fixed and our
asymptotic results will be obtained with N (the number of experiments) tending
to ∞.
Let Θ ⊂ Rd and f : Θ × N2 → (0, 1) be a function, called the choice function
of the urn. We assume that there exists θ0 ∈ Θ such that f0 (·, ·) = f (θ0 , ·, ·),
i.e. for k = 0, . . . , n − 1 and i = 0, . . . , k,
j
P(Xk+1
= 1 | Zkj = i) = f (θ0 , i, k − i) .

To proof the consistency and the asymptotic normality of the estimators introduced below, we need the following assumptions on the choice function f . For
any function g defined on Θ, we denote ġ and g̈ the gradient and Hessian matrix
with respect to θ, ∂s g the partial derivative with respect to the s-th component θs
of θ, 1 ≤ s ≤ d and A0 the transpose of the vector or matrix A.
5

Assumption 1.
(i) (Regularity) The set Θ is a compact with non empty interior. For 0 ≤
i ≤ k ≤ n − 1, f0 (i, k − i) > 0 and the function θ → f (θ, i, k − i) is twice
continuously differentiable on Θ.
(ii) (Identifiability) If f (θ1 , i, k − i) = f (θ2 , i, k − i) for all 0 ≤ i ≤ k ≤ n − 1,
then θ1 = θ2 ,
(iii) The n(n − 1)-dimensional vectors {∂s f (θ0 , i, k − i), 0 ≤ i ≤ k ≤ n − 1},
1 ≤ s ≤ d, are linearly independent in Rn(n−1) .
Assumption 1 ensures that θ0 is the unique maximizer of L and that the Fisher
information matrix
In (θ0 ) = −L̈(θ0 ) =

n−1 X
k
X
k=0 i=0

P(Zk = i)
f˙0 (i, k − i)f˙0 (i, k − i)0
¯
f0 (i, k − i)f0 (i, k − i)

is invertible, where we denote f˙0 (i, k − i) = f˙(θ0 , i, k − i) and f¯0 (i, k − i) =
1 − f0 (i, k − i). The explicit expression of the probabilities P(Zk = i), i, k ∈ N, is
given in Section 6.1, Equation (13).

2.1

Maximum likelihood estimation (MLE)

The structure of the model (3) allows to have an explicit expression of likelihood VN (θ). The independence of the N experiments yields the following multiplicative form
VN (θ) =

N n−1
Y
Y

j

j

f (θ, Zkj , k − Zkj )Xk+1 {1 − f (θ, Zkj , k − Zkj )}1−Xk+1 .

j=1 k=0

The log-likelihood function LN based on N paths, is thus given by
LN (θ) =

N X
n−1
X


j
Xk+1
log f (θ, Zkj , k − Zkj )

j=1 k=0
j
+ (1 − Xk+1
) log{1 − f (θ, Zkj , k − Zkj )} . (4)

Let θ̂N be the maximum likelihood estimator of θ0 , that is
θ̂N = arg max LN (θ) .
θ∈Θ

6

(5)

Define L(θ) = N −1 E[LN (θ)] (where the dependence in n is omitted). Then,
L(θ) =

n−1
X

E f0 (Zk , k − Zk ) log f (θ, Zk , k − Zk )
k=0


+ {1 − f0 (Zk , k − Zk )} log{1 − f (θ, Zk , k − Zk )} . (6)
Let N (m, Σ) denote the Gaussian distribution with mean m and covariance Σ.
Theorem 2. If Assumptions 1-(i) and 1-(ii) hold then the maximum likelihood
estimator θ̂N is a strongly consistent estimator of θ0 . If moreover√Assumption 1-(iii)
holds and θ0 is an interior point of Θ, then, as N tends to ∞, N (θ̂N − θ0 ) converges weakly towards N (0, In−1 (θ0 )).
The proof is in Section 6.4. It is the consequence of a more general result
stated and proved therein.

2.2

Weighted least squares estimation (WLSE)

Least squares estimators are very popular among practitioners. Moreover for some
urn models, the MLE may be numerically unstable hence difficult (and lengthy) to
compute so the WLSE constitutes a convenient alternative. In order to describe
this estimator, we introduce some notation. For 0 ≤ i ≤ k ≤ n − 1, define
PN
j
1
N
1 X
j=1 1{Zkj =i} Xk+1
N
, (7)
aN (i, k − i) =
1 j
and pN (i, k − i) =
N j=1 {Zk =i}
aN (i, k − i)
with the convention 00 = 0. The quantity aN (i, k − i) is the empirical probability
that i red balls have been drawn at time k and pN (i, k − i) is the empirical
conditional probability that a red ball is again chosen at time k + 1 given i red
balls were drawn at time k.
We further define qN (i, k−i) = 1−pN (i, k−i) and f¯(θ, i, k−i) = 1−f (θ, i, k−i).
Let {wN (i, k − i), 0 ≤ i ≤ k} be a sequence of weights and define the contrast
function
WN (θ) =

n−1 X
k
X

wN (i, k − i){pN (i, k − i) − f (θ, i, k − i)}2 .

(8)

k=0 i=0

The weighted least squares estimator minimizes WN , that is
W
θ̂N
= arg min WN (θ) .
θ∈Θ

7

(9)

Theorem 3. If Assumptions 1-(i) and 1-(ii) hold and if the weights wN converge
almost surely to a sequence of positive weights w0 , then the weighted least squares
W
estimator θ̂N
is a strongly consistent estimator of θ0 .
If moreover
1-(iii) holds and θ0 is an interior point of Θ, then as N
√ Assumption
W
tends to ∞, N (θ̂N − θ0 ) converges weakly towards N (0, Σn (θ0 )), where Σn (θ0 )
is a definite positive covariance matrix.
If moreover wN (i, k − i) = pN −1 qN −1 aN (i, k − i), for all 0 ≤ i ≤ k ≤ n − 1, the
W
estimator θ̂N
is asymptotically efficient, i.e. Σn (θ0 ) = In (θ0 ).
The proof is in Section 6.4, where a explicit expression of Σn (θ0 )−1 is supplied.

2.3

Generalization

It is possible to extend the statistical framework introduced in the previous
sections to a large class of reinforced random walks on graphs. For instance,
let G = (X
 , E) be a locally finite non oriented graph, with X the set of its vertices
and E ⊂ {i, j} : i, j ∈ X , i 6= j the set of its non oriented edges. We denote
by x ∼ y, if {x, y} ∈ E. Let NX be the set of integer vectors indexed on X . We
define X = (Xn )n≥0 a random walk on G, i.e. a sequence of vertices, such that,
for all n ≥ 0, Xn+1 ∼ Xn . The vector Zn ∈ NX is such that, for all vertex x ∈ X ,
Zn (x) is the number of times the walk X has visited x up to time n. Let Θ ⊂ Rd ,
d ≥ 1, we suppose that the walk X is vertex-reinforced and that there exists a
choice function f : Θ × X × NX → [0, 1] and a parameter θ ∈ Θ such that, for
all n ≥ 0 and x ∈ X ,
P(Xn+1 = x|X0 , · · · , Xn ) = f (θ, x, Zn )1{x∼Xn } .
For instance, the choice function can be similar to the one proposed by J.-L.
Deneubourg (see Deneubourg et al., 1990), for n ≥ 0 and x ∈ X ,
α
c + Zn (x)
P(Xn+1 = x|X0 , · · · , Xn ) = X
α 1{x∼Xn } ,
c + Zn (y)
y∼Xn

where Θ = (0, +∞)2 and θ = (α, c).
The statistical framework introduced for the general urn model is easily adaptable to this vertex reinforced random work. Under some adequate regular conditions, it would be not difficult to prove the consistency and the asymptotic
normality of the MLE and the WLSE by establishing a theorem similar of the
general result proved in Section 6.3.

8

3

Application to an ethological problem

In 1990, J.-L. Deneubourg et al. used a particular urn model to reproduce the
sequences of consecutive choices made by ants at a fork (see Deneubourg et al.,
1990). Let replace the urn filled with balls of two colors by a fork with two
branches. Drawing a ball and adding a ball of the same color in the urn is
equivalent to an ant choosing a branch and reinforcing it with pheromone by going
throw it. Then the probability to draw a red ball depending on the previous draws
is equal to the probability to choose the right branch depending on the previous
passages.
We first introduce the model proposed by J.-L. Deneubourg et al. in the statistical framework described in the previous section. Futher, we provide a description of the model behavior depending on the parameter value. We also supply an
ethological interpretation of the parameter. We then prove that it is impossible
to estimate the parameter on a single path. Finally we report the study of the
estimator accuracy that we perform on simulated data.

3.1

Formalization and results

Behavioral assumptions We first introduce the hypotheses assumed on the
ant behavior.
1. Each ant regularly deposits a constant amount of pheromone as it walks.
2. The ants are strictly identical which means that every ant has the same
reaction to the same amount of pheromone.
3. Pheromone trails do not evaporate during the experiment.
4. Each ant reaches the fork alone, chooses a branch and leaves the bifurcation
by crossing only once into the chosen branch without passing through or
reinforcing the non-chosen branch.
Under these assumptions, the quantity of pheromone laid on each branch is proportional to the number of passages through it. In path formation modeling, these
assumptions are commonly made. However, because of the inter-individual variability in ant behavior, the first two assumptions are unrealistic. For instance, the
pheromone perception noise implies that each ant could detect a different signal
from the same quantity of pheromone. These assumptions are an approximation
of the real ant behavior. The implicit hypothesis here is that the inter-individual
variability is small enough to consider that all ants are identical. For the third
assumption, we suppose that the persistence of the pheromone trails allows to
ignore the evaporation of the pheromone. Experimental protocols are designed to
9

make the four assumptions more acceptable by choosing the ant species adequately
and by placing them in an appropriate situation (see Section 4).
The model The random variable Xk , introduced previously, is the choice of
the k-th ant going through the fork (1 for right and 0 for left). Consequently, for
k ≥ 1, Zk is the number of passages through the right branch after k passages.
For θ = (α, c) ∈ (0, ∞)2 and all integers 0 ≤ i ≤ k, we define choice function
f (θ, i, k − i) =

(c + i)α
.
(c + i)α + (c + k − i)α

(10)

We assume that the probability that an ant chooses the right branch at time k + 1
given the first k choices is given by:
P(Xk+1

(c + Zk )α
= 1|Fk ) = f (θ, Zk , k − Zk ) =
,
(c + Zk )α + (c + k − Zk )α

(11)

where c > 0 is the intrinsic attractiveness of each branch and α > 0 is the
possible non-linearity of the choice. This process is an urn model which has been
exhaustively investigated in the probabilistic literature. We recall here its main
features.
Theorem 4.
(i) If α < 1, then
lim

n→∞

1
Zn
= a.s.
n
2

(ii) If α = 1, then Zn /n converges almost surely to a random limit with a
Beta(c, c) distribution with density x → Γ−2 (c)Γ(2c)xc−1 (1 − x)c−1 with respect to Lebesgue’s measure on [0, 1], and Γ is the Gamma function.
(iii) If α > 1, then eventually only one branch will be chosen, i.e.
∃x ∈ {0, 1} , ∃n0 ∈ N , ∀n ≥ n0 , Xn = x .
The case α < 1 is due to Tarrès (2011); the case α > 1 to Davis (1990) and
the case α = 1 to Pólya (1931) (see Freedman (1965) for an online access).
Ethological interpretation of the parameters and properties of the model
The parameter α characterizes the ant’s differential sensitivity to the pheromone.
When α > 1, the ants can detect better and better increasing amounts of pheromones laid on each branch and thus are more likely to choose the branch with the
10

most pheromones. Moreover, after a random but almost surely finite number of
passages, one branch will eventually be selected, i.e. all ants will afterwards choose
this branch (see Theorem 4-(iii)). In the opposite case, when α < 1, the ants are
less able to perceive the differences between the amounts of pheromones laid on
each branch as these amounts increase. This minimization effect is so strong that
the proportion of passages on each branch converges to 1/2 (see Theorem 4-(i)). It
is important to note that, when α 6= 1, the asymptotic behavior of the proportion
of passages through each branch depends only on α and not on c. Furthermore
the larger α is or the closer α is to zero, the faster these effects will happen.
The role of c is clearer when (11) is rewritten as follows:
P(Xk+1 = 1|Fk ) =

(1 + Zk /c)α
.
(1 + Zk c)α + (1 + (k − Zk )/c)α

The parameter c is the inverse of the reinforcement incrementation and thus can be
interpreted as the inverse of the attractiveness (or the strength) of the pheromones
laid at each passage. Consequently when α is neither very close to zero nor very
large, c has a strong short term influence. When c is small compared to 1, the
first passage strongly reinforces the first chosen branch. Thus during the first
few passages, a branch will be highly favored even if α < 1 (in which case the
branches will eventually be uniformly crossed). When c is large, the first passages
weakly reinforce the chosen branches. Then if α > 1 (in which case a branch will
eventually be selected), a large number of passages must be observed before the
clear emergence of a preference. Naturally, the larger c is or the closer c is to zero,
the longer these effects will be seen.
When α = 1, the asymptotic behavior of the passage proportion is determined
by c (see Theorem 4-(ii)). As c grows from zero to infinity, the limiting distribution
of Zn /n (as n → ∞) evolves continuously from two Dirac point masses at 0 and 1
to a single Dirac mass at 1/2. To illustrate this point, we show in Figure 1 the
density of the Beta distribution for c = 0.1 and c = 10. We make some further
comments.
• If c < 1, a strong asymmetry in the choices of the branches appears. One
branch is eventually chosen much more frequently than the other. Furthermore as c tends to 0, the Beta distribution tends to the distribution with
two point masses at 0 and 1. This limit case corresponds to the situation in
which a branch is selected, i.e. α > 1.
• If c = 1, the limiting distribution is uniform on [0, 1].
• If c > 1, Zn /n appears to be much more concentrated around 1/2. This is
similar to what is observed in the case α < 1.
11

Figure 1: Graph of the Beta distribution for parameters (c, c) with c = 0.1, c = 1
and c = 10
To summarize, the model possesses four phases: fast and slow selections and
fast and slow uniformizations. These phases are delimited by two phase transitions: a discontinuous one between α > 1 (branch selection) and α < 1 (branch
uniformization) with a critical state α = 1 and a smooth one between c < 1 (strong
pheromone deposits) and c > 1 (weak pheromone deposits). These properties are
summarized in the phase diagram in Figure 2. When c is small, it is very likely
that one branch will be favored during the first passages, thus the empirical distribution of the choices resembles the Beta distribution with a small c (the grey
solid line in Figure 1). As the number of experiments increases, the shape of the
empirical distribution will be closer and closer to its limit: a Dirac mass at 1/2
if α < 1 and two Dirac masses at 0 and 1 if α > 1. When c is large, the earlier passages do not show any preference between the branch. Again, when the number
of experiments increases, the asymptotic behavior is progressively revealed.
To date, the most commonly used behavioral state in the model is the slow selection of a branch (see Deneubourg et al., 1990; Beckers, Deneubourg, and Goss,
1992). But at least two other states are interesting. The fast uniformization
can describe the case where none of the branches are preferred. The slow uniformization could reproduce the saturation phenomenon. There exists a threshold
concentration of pheromone upon which ants can no longer detect the pheromone
concentration variations (see Pasteels, Deneubourg, and Goss, 1987). In experiments involving many ants, one can first observe the favorization of a branch.
But when this branch is saturated (its attractiveness stops increasing), ants go
more and more through the other branch, whose attractiveness still increases.
Eventually, the two branches are uniformly chosen.

12

Remark 5. In the context of an estimation procedure, the similar effects of α
and c (favorization/selection of a branch or not) induce an identifiability issue.
Indeed, we observe a finite number of choices and consequently we only see the
short term behavior. We have seen that the favorization of a branch in the first
passages could be due to a pair of parameter values (α, c) with c small compared
to 1 and α close to 1 or to a pair (α, c) with c close to 1 and α large compared to 1.
On the other hand, if the first passages are nearly uniform on the two branches,
it could be the result of c close to 1 and α small, or of c large and α close to 1.
Thus, we can expect that the estimation of the parameters will be difficult
when both parameters contribute to the same effect, e.g. α large and c small (fast
selection of one branch) or α small and c large (no selection); and also when
the parameters have competing effects: very small c and α < 1, or very large c
and α > 1. The statistical procedure that we introduce in this paper partially
circumvents this difficulty, since it focuses on the transition probabilities rather
than on the general shape of a curve, which is what calibration methods do. This
will be illustrated in Section 3.4.

Figure 2: Phase diagram of the model. The graphs in the shaded boxes show
the shape of the empirical distribution of Zn /n for small n (left) and its limiting
distribution (as n → ∞, right).

13

3.2

Parameter estimation

We start with the maximum likelihood estimator, that is θ̂N = (α̂N , ĉN ) defined
by (5). The Fisher information matrix has the following expression.
In (θ) =

n−1 X
k
X

P(Zk = i)f (α, c, i, k − i)f¯(α, c, i, k − i)J (α, c, i, k − i) ,

(12)

k=0 i=0

with θ = (α, c) and for 0 ≤ i ≤ k ≤ n − 1,

c+i
log
log2 c+k−i
 α(k−2i)
J (α, c, i, k − i) =
c+i
log c+k−i (c+i)(c+k−i)

α(k−2i)
c+i
c+k−i (c+i)(c+k−i)
α2 (k−2i)2
(c+i)2 (c+k−i)2



!
.

It is important to note that the Fisher information matrix is not diagonal. Thus
the estimation of each parameter has an effect on the estimation of the other.
Corollary 6. Let Θ be a compact subset of (0, ∞)2 which contains (α0 , c0 ). Then
the maximum√
likelihood estimator θ̂N is consistent and asymptotically normal and
efficient, i.e. N (θ̂N − θ0 ) converges weakly to N (0, In−1 (θ0 )).
As mentioned above, we also use weighted least squared estimators defined
by (9) for several different weight sequences wN , such that wN (i, k − i) converges
almost surely to w0 (i, k − i) > 0, for all 0 ≤ k ≤ n − 1 and 0 ≤ i ≤ k.
Corollary 7. Let Θ be a compact subset of (0, ∞)2 which contains (α0 , c0 ) and
assume that wN converges almost surely to positive weights w0 . Then the
√weighted
W
W
least squared estimator θ̂N
is consistent and asymptotically normal, i.e. N (θ̂N
− θ0 )
converges weakly to N (0, Σn (θ0 )), where Σn (θ0 ) is a positive definite covariance
−1
matrix. It is efficient, if wN (i, k − i) = p−1
N qN aN (i, k − i), for all 0 ≤ k ≤ n − 1
and 0 ≤ i ≤ k.
The proofs of these corollaries are in Section 6.5. They are an application of
Theorems 2 and 3. Since Assumption 1-(i) obviously holds, it remains only to
check (ii) and (iii) of Assumption 1.

3.3

Estimation on a single path

The main feature of the binary choice model for α0 > 1 is that only one branch
will be crossed eventually. It seems clear then that a statistical procedure based
on only one path (one sequence of choices) cannot be consistent, since no new
information will be obtained after one branch is eventually abandoned. This
intuition is true and more surprisingly, it is also true in the case α0 = 1. This
is translated in statistical terms in the following theorem. Let `n denote the loglikelihood based on a single path of length n and `˙n its gradient. The model is
regular, so the Fisher information is varθ (`˙n (θ)).
14

Theorem 8.
(i) If α0 = 1 and c0 > 0, then limn→∞ In (c0 ) < ∞.
(ii) If α0 < 1, n−1 `n (θ0 ) → − log 2.
(iii) If α0 > 1, then `n (θ0 ) converges almost surely to a random variable as n → ∞.
The proof is in Section 6.6. Statement (i) means that, when α0 = 1, the
Fisher information is bounded. This implies that the parameter c0 cannot be
estimated on a single path. This also implies that the length n of each path
should be taken as large as possible (theoretically infinite) in order to minimize
the asymptotic variance of the estimators. Statements (ii) and (iii) imply that the
maximum likelihood estimator is inconsistent, since the likelihood does not tend
to a constant.

3.4

Simulation experiment

In order to assess the quality of the estimators proposed, we have made a short
simulation study. For several pairs (α, c), we have simulated 1000 experiments
of N = 50 paths of length n = 100 (recall that n is the number of ants going
through the bifurcation). These are reasonable values in view of the practical
experiments with actual ants. We compare the performance of the maximum
likelihood estimator θ̂N (MLE) defined in (5) and of the weighted least squares
W
estimator θ̂N
(WLSE) defined in (9) with the weights wN (i, k − i) = aN (i, k − i)
defined in (7). The asymptotically efficient WLSE with the weights wN (i, k − i) =
aN (i, k − i)pN (i, k − i)−1 qN (i, k − i)−1 provides a severely biased estimation of α
and always estimates a very small value of c with a very small dispersion. This
is caused by the fact that the empirical pN and qN vanish frequently, so that the
weights are infinite. We will not report the study for this estimator.
The theoretical standard deviation We first evaluate numerically some values of the theoretical standard deviations of both estimators for several values of α
and c. We have chosen arbitrary values of α and c in the range 0.5, 2. We have
also chosen values of α and c which correspond to those found in the literature
cited and to those that we have estimated in the real life experiment described in
Section 4. These results are reported in Table 1 and in Figure 3 and their features
are summarized in the following points.
• As theoretically expected, the asymptotic variance of the MLE, which is the
Fisher information bound, is smaller than the variance of the WLSE, but
the ratio between the variances of the two estimators is never less than one
fourth. Moreover, their overall behavior is similar.
15

• The variance of the estimators of α is smaller when both parameters do not
contribute to the same effect. The worst variance is for α large and c small,
that is when the values of both parameters imply fast selection of a branch.
The variance tend to infinity when α tends to infinity.
• The variance of the estimators of c increases with c and tends to infinity
when α tends to 0 and to ∞.
• These effects are explained by the fact that the coefficients of the Fisher
information matrix tend to zero when α tends to zero, except the coefficient
corresponding to α. See Formula (12).

(a) TSD of α

(b) TSD of c

Figure 3: Theoretical standard deviation for N = 50 paths of length n = 100 of
the MLE of α (a) and of c (b), for α in (0, 2] and fixed values of c. Theoretical
standard deviation of the WLSE has the same shape, but with higher convergence
speed as α tend to 0 or ∞.

Performance of the estimators Recall that we have simulated 1000 experiments, each of N = 50 paths of length n = 100. Because of the length of the
computations, each MLE was computed only 500 times. Table 1 reports root
mean squared error (MSE) of both estimators based on the simulated data for the
same values of the parameters and their features are summarized in the following
points.
• For most values of the parameters, the MSE are close to the theoretical
standard deviation.

16

• The MSE increase significantly when c is large or when both parameters
contribute to the same effect.
• This increase is more noticeable for the WLSE than for the MLE.
• This increase is in part due to the skewness of these estimators. For some
values of the parameters, both estimators tend to overestimate the parameters.
• For the MLE, the MSE is much larger in the case of non selection than in
the case of selection where the empirical performance of the MLE nearly
matches the theoretical value.
• These effects are always stronger for the estimation of c than for the estimation of α.
This degraded performance for some specific or extreme values of the parameters is in part due to numerical issues.
• In the case where selection of a branch is fast, many of the empirical weights
used to compute the WLSE vanish, and the least squares method uses very
few points to fit the curve. The MLE is not affected by this problem.
• In the case where both parameters concur to non selection, the probability
of choosing one branch converges very fast to 1/2, and thus the experiment
brings very little information. This affects both the MLE and the WLSE,
and in addition, many of the empirical weights vanish so the WLSE is even
less efficient.
The degraded performance may also be caused by the identifiability problem
explained in Remark 5, i.e. the similar effects of the two parameters makes the
estimation more difficult.
Bootstrap confidence intervals Since the asymptotic variance depends on
the unknown parameters, we have computed the pivotal Bootstrap 95% confidence
intervals for the parameters based on one simulation of N = 50 paths of length
n = 100 and a Bootstrap sample size of 500 (see Wasserman (2004), Section 8.3,
for details on this method). We have compared these Bootstrap intervals with
the corresponding Monte-Carlo intervals, based on 500 simulations (see Table 2).
The match is nearly perfect for the MLE for α, but as before, the performance is
poorer for the estimation of c. The intervals for c are noticeably skewed to the
right but always contain the true value. For further comparison, we only show
here the results corresponding to the values of the parameters estimated in the
real life experiment reported below and those corresponding to values found in
the earlier literature.
17

Table 1:
Theoretical standard deviations (TSD) for N = 50 paths of
length n = 100 and square root of the mean squared errors (MSE) for 500 (for
the MLE) or 1000 (for the WLSE) simulated experiments of N = 50 paths of
length n = 100. All figures of this table must be mutiply by 0.01.
(α, c)
(0.5, 0.5)
(0.5, 1.0)
(0.5, 2.0)
(1.0, 0.5)
(1.0, 1.0)
(1.0, 2.0)
(1.5, 0.5)
(1.5, 1.0)
(1.5, 2.0)
(2.0, 0.5)
(2.0, 1.0)
(2.0, 2.0)
(2.0, 20.0)
(2.6, 60.0)
(1.1, 3.0)
(1.1, 7.0)

α
c
MLE
WLSE
MLE
WLSE
√
√
√
√
TSD
M SE TSD
M SE TSD
M SE TSD
M SE
5.02
5.60
6.25
7.60
25.4
54.1
31.7
58.8
6.45
7.39
8.56
10.2
59.8
137
79.4
169
8.80
16.9
12.9
356
146
651
214
29500
3.81
3.90
4.97
6.25
11.8
12.4
15.9
21.0
4.34
4.98
5.91
7.34
25.2
31.3
34.9
47.3
5.83
5.80
8.72
9.38
59.3
64.9
89.1
111
7.83
7.85
12.8
19.4
12.3
13.0
18.7
24.6
6.69
6.59
10.8
15.1
21.0
21.4
33.1
40.3
6.88
7.24
11.5
13.9
41.7
45.9
68.8
78.4
19.4
26.1
41.8
2690
16.4
19.8
28.6
1410
13.5
14.2
28.7
985
23.9
26.0
44.3
988
11.2
11.7
23.1
35.4
40.7
43.4
77.2
101
35.5
58.8
129
124
794
1420
2880
3240
166
864
1230
3590
5780
32500
43200 142000
7.20
7.80
11.9
11.3
89.9
111
148
201
13.8
17.4
29.5
28.3
297
433
633
758
All figures must be multiply by 0.01

Table 2: Monte-Carlo 95% confidence intervals for 500 simulated experiment of
N = 50 paths of length n = 100 and Bootstrap 95% confidence intervals for one
simulated experiment of 50 paths of length 100
IDC 95% for α
IDC 95% for c
(α, c)
Monte-Carlo Bootstrap Monte-Carlo Bootstrap
MLE (2.0, 20) (1.50, 3.38) (1.90, 5.10) (9.90, 55.5) (19.1, 94.4)
(2.6, 60) (1.23, 29.3) (1.12, 38.7) (15.9, 1053) (19.9, 1637)
(1.1, 3.0) (0.98, 1.29) (1.01, 1.23) (1.81, 6.22) (1.34, 3.99)
WLSE (2.0, 20) (1.26, 4.85) (1.06, 3.99) (5.92, 88.1) (2.76, 69.3)
(2.6, 60) (0.74, 88.3) (0.24, 87.7) (3.86, 3754) (0.17, 4132)
(1.1, 7.0) (0.75, 1.83) (0.55, 2.56) (1.68, 30.3) (1.48, 66.0)

18

4

Real life experiment with ants

In this section, we apply the previous estimators on data from a path selection
experiment by a colony of ants.

4.1

Experiment description

This experiment was done in the Research Center on Animal Cognition (UMR 5169)
of Paul Sabatier University Toulouse under the supervision of Guy Theraulaz,
Hugues Chaté and the first author. A small laboratory colony (approximately
200 workers) of Argentine ants Linepithema humile was starved for two days before the experiment. During the experiment, the colony had access to a fork
carved in a white PVC slab, partially covered by a Plexiglas plate (see Figure 4).
The angle between the branches was 60◦ . The fork galleries had a 0.5 cm square
section. The entrance of the maze was controlled by a door. Food was never
present during the experiment. The maze was initially free of any pheromone
trail.

Figure 4: The experimental set-up: a fork carved in a white PVC slab, partially
covered by a Plexiglass plate.
Each trial (N = 50) consisted in introducing separately each ant to the entrance of the fork (see Figure 4) one at a time. Once inside, an ant must choose
between the left or the right branch of the fork. As soon as the ant had made a
choice and stepped into one branch, it was removed from the set-up and another
ant was introduced. All the choices were recorded and a trial ended when 100 ants
had passed through the fork.
This experimental protocol was designed to strengthen the behavioral assumptions described in Section 3.1. Any return to the fork is forbidden so that we can
19

consider that each ant passed only one time. There was never more than one
ant in the set-up. This implies that each ant in the maze received no other
cue about the previous passages than the pheromone that was been laid. The
species Linepithema humile was in part chosen to justify the assumption of identical pheromone deposits. Indeed, these ant may deposit regularly the same type
of pheromone on their trajectory (see Van Vorhis Key and Baker, 1982; Aron,
Pasteels, and Deneubourg, 1989). All ants were prepared the same way before
the experiments to increase the credibility of the assumption stating that each
ant behaved by the same way. The length of the experiments was limited to stay
close to the half-life duration of the pheromone trails (see Jeanson et al., 2003).

4.2

Data representation

Figures 5(a) shows the 50 paths of length n = 100, that is, 50 choice sequences of
100 ants that went through the bifurcation. The paths are represented as random
walks with increment +1 when the right branch is chosen, and −1 when the left
one is chosen. In less than ten experiments, a branch seemed to be selected,
whereas in the others, selection of a branch was not obvious. Figure 5(b) shows
the histogram of the distribution of Z100 /100, that is the final proportion of the
choices of the right branch. There is no clear visual evidence that α > 1 as it is
claim in the literature (see Deneubourg et al., 1990; Vittori et al., 2006; Garnier
et al., 2009).

(a) The 50 paths of n = 100 ants choosing
either left (+1) or right (-1)

(b) Histogram of the final proportion of right
passages (Z100 /100)

Figure 5: Data representation

20

4.3

Parameter estimation

Several values of these parameters have been proposed in the applied literature.
Deneubourg et al. (1990) proposed α = 2, c = 20 and more recently Garnier
et al. (2009) suggested α = 2.6 and c = 60. It must be noted however that
these values are not obtained by a statistical method but by the calibration of
a curve to a plot. Therefore, these methods do not lead to confidence intervals.
Moreover, a calibration method has an inherent risk of over fitting, because of the
identifiability problem explained in Remark 5. As illustrated in Figure 2, if for
instance α and c are both small, then both branches will be asymptotically equally
chosen, but paths of finite length n might be misleading and the calibration will
suggest values of α and c corresponding to the selection of a branch. The statistical
procedure is based on the dynamics of the process and is thus less prone to this
type of error. Nevertheless, we will see that our results do not contradict those of
Deneubourg et al. (1990) and Garnier et al. (2009), but complement them.
Table 3: The MLE and the WLSE for the 50 paths of real ants and their
Boostrap 95% confidence intervals
α̂
Bootstrap 95% CI
ĉ
Bootstrap 95% CI
MLE 1.07
(0.80, 1.99)
3.26
(1.14, 23.0)
WLSE 1.10
(0.62, 3.81)
6.91
(0.94, 85.4)
Table 3 shows the results of the maximum likelihood estimation and the
weighted least squares estimation. Both estimates of α are close to 1.1 and the
estimates of c are between 3 and 7. The 95% confidence intervals are slightly
larger than the simulated ones (see Table 2). This increased variability may be
due to the extreme paths which seem to show a very fast selection of one branch
(see Figures 5(a) and 5(b)). This may suggest that the ants did not have the
same behavior and that the distribution of Z100 /100 could be a mixture of two
distributions.
For both methods, the 95% Bootstrap confidence intervals of α contain the
value 1. More precisely, as shown in Figure 6, approximately 1/3 of the bootstrap parameters gives weak pheromone deposits (c > 1) and a weak differential
sensitivity (α < 1), which means that branches are eventually uniformly crossed.
In almost all the others cases, we conclude for weak pheromone deposits (c > 1)
and a strong differential sensitivity (α > 1), which means that a branch will be
eventually, though slowly, selected. In only a few cases do the estimators give
strong pheromone deposits (c < 1), but a weak differential sensitivity (α < 1),
which means that a branch is chosen more than the other at the beginning of the
experiment, but branches are eventually uniformly crossed. Finally, there are no
values which imply both strong pheromone deposits (c < 1) and a strong diffe21

rential sensitivity (α > 1). Therefore, we can conclude that pheromone deposits
are weak with a good confidence but we cannot confidently decide for α.
The values obtained by Deneubourg et al. (1990) (α = 2, c = 20) and more
recently by Garnier et al. (2009) (α = 2.6, c = 60) are both in the confidence
intervals for the WLSE found in Table 3. But the values of α suggested by
these authors are out of the 95% confidence interval for the MLE. Thus these
parameters, which decide for a slow branch selection, are no more likely than a
parameter set which would yield non selection of a path.
Figure 6 illustrates the fact that the two estimators are strongly positively
correlated. There seems to be two cutoff values for c: if ĉ∗ > 8, then α̂∗ > 1, and
if ĉ∗ < 1.5, then α̂∗ < 1. The above mentioned values reported by Deneubourg
et al. (1990) and Garnier et al. (2009) exhibit these features: they both have c > 8
and α > 1 and α increase with c.
If we fix the value of c and estimate only α, then the 95% Bootstrap confidence
intervals for α are smaller. Figure 7 shows the estimated values of α and the
confidence intervals as functions of the fixed value of c. We see that if c is greater
than 6 for the MLE (or than 12 for the WLSE), then the confidence intervals of α
lie entirely above 1. Furthermore if c is less than 0.8 for the MLE (or than 2 for the
WLSE), then the confidence intervals of α lie entirely under 1. This shows that if
the deposits are weak enough, i.e. c > 12, we can conclude that a slow selection
of a branch will occur with probability 1. On the other hand, if the deposits
are strong enough, i.e. c < 0.8, we can conclude that branches will eventually be
uniformly crossed with probability 1.

(a) MLE

(b) WLSE

Figure 6: Log-log scatterplots of the estimates (α̂∗ , ĉ∗ ) for the 500 Bootstrap
samples for the MLE (a) and the WLSE (b) for the 50 paths of real ants.

22

(a) MLE

(b) WLSE

Figure 7: Graph of the estimates α̂ and their Bootstrap 95% confidence interval
for the MLE (a) and the WLSE (b) for the 50 paths of real ants as a function of
fixed value of c.

5

Concluding remarks

In the literature no parameter estimation methods for reinforced random walks
can be found. To partially fill this void, this article proposes a statistical framework to estimate the parameter of a general two-colored urn model. We define
the maximum likelihood estimator (MLE) and the weighted least squares estimators (WLSE) for the parameter of this model and prove their consistency and
their asymptotically normality under some usual regularity conditions. The proof
lies on a general result for a large class of estimators called minimum contrast
estimators. The MLE is asymptotically efficient, but can be difficult (lengthy)
to compute, which can be an issue specially while using Bootstrap algorithms.
The WLSE is a suitable alternative. Moreover this estimator is popular among
practitioners.
We apply this statistical tools to the problem of path selection by an ant
colony. To this purpose, we performed experiments with actual ants to collect
data. The experiment consisted of introducing one hundred ants into a Y shaped
device, one at a time, and observing their successive choices. We also consider
the particular urn model introduced by Deneubourg et al. (1990) to describe this
phenomenon. This urn has two parameters, α and c, which have distinct biological
interpretations, but contribute to the same effect: either selection of a branch or
uniformization of the choices. The parameter c influences the short term behavior,
whereas α determines the asymptotic behavior. Consequently the model exhibits
four phases which are illustrated by Figure 2. The case most commonly considered
in the literature is the case of slow selection, which corresponds to α > 1 and c > 1:
the ants will eventually always choose the same branch, but this selection will take
23

a long time. For instance, Deneubourg et al. (1990) provides the values α = 2
and c = 20. However other phases can be relevant to describe the ant behavior.
For instance the fast uniformization, corresponding to α < 1 and c > 1, can model
the less likely but not negligible case in which ants do not select a branch.
After assessing the accuracy of the MLE and the WLSE on simulated data,
we estimate the value of α and c with the two estimators. We also evaluate
confidence regions by Bootstrap proceeding. The estimated values of α and c
ranged between 1.1 and 3 and between 3 and 7, respectively. This tends to imply
that slow selection of a branch will occur. However, the Bootstrap sample gives a
confidence level of 65% for the hypothesis of slow selection, while the hypothesis
of fast uniformization has a confidence of 35%.
This low level of confidence for the commonly assumed slow selection phase
might be explained by technical reasons. The number of experiments (50) is relatively small; increasing the number of replicas will reduce the confidence regions.
Moreover the competition between the parameters for the same effect induces an
identifiability issue. For instance the apparent preference of a branch may be due
to α > 1 or c small with respect to 1. Therefore, the model, which is biologically
relevant, is statistically difficult to estimate. Indeed, for an ethological study, discriminating the ant pheromone sensitivity from the pheromone deposit strength
is meaningful. But for a statistical procedure, the similarity of effect of the two
parameters scales down the estimation performance.
However, the uncertainty may not come from an inefficiency of the statistical
procedure, but from shortcomings of the ethological hypotheses. Indeed, the estimated confidence intervals computed from the experimental data are larger than
the ones computed from the simulated data (for similar parameter values). Moreover the assumption that the inter-individual variability is negligible is strong.
For instance, it may be necessary to consider that the pheromone deposit varies
at each passage, i.e. that c is random.
These ethological considerations will be further discussed in a forthcoming
paper which will analyze more elaborated experimental designs. The ants will be
observed while freely evolving in a network with several nodes. In addition of a
data analysis, we will model the experiment with a reinforced random walk on
a finite graph for which we have provided probabilistic results (see Le Goff and
Raimond, 2015). The statistical methodology introduced in this paper will be
extended to a larger class of reinforced random walks.

24

6
6.1

Proofs
Distribution of Zk , for k ∈ N

In order to compute the distribution of Zk , we introduce some notation. Let Sk
be the set of sequences of length k + 1 of integers i0 , . . . , ik such that i0 = 0 and
ij − ij−1 ∈ {0, 1} for j = 1, . . . , k. For i ≤ k let Sk (i) = {(i0 , . . . , ik ) ∈ Sk | ik = i}.
Then we have
P(Zk = i) =

X

k−1
Y

f0 (iq , q − iq )iq+1 −iq (1 − f0 (iq , q − iq ))1−iq+1 +iq . (13)

(i0 ,...,ik )∈Sk (i) q=0

6.2

A central limit theorem for the empirical conditional
probabilities

For 0 ≤ i ≤ k ≤ n − 1, recall the definition of aN (i, k − i) and pN (i, k − i) in (7)
and that f¯0 (i, k − i) = 1 − f0 (i, k − i).
√
Lemma 9. { N (pN (i, k − i) − f0 (i, k − i)), 0 ≤ i ≤ k ≤ n − 1} converges weakly
to a Gaussian vector with diagonal covariance matrix Γ0 with diagonal elements
γ0 (i, k − i) =

f0 (i, k − i)f¯0 (i, k − i)
.
P(Zk = i)

(14)

P
j
Proof. Define bN (i, k − i) = N −1 N
j=1 1{Zkj =i} Xk+1 , the empirical estimate of
b(i, k − i) = P(Zk = i, Xk+1 = 1) and a(i, k − i) = E[aN (i, k − i)] = P(Zk = i).
Write then
pN (i, k − i) − f0 (i, k − i) =

bN (i, k − i) − b(i, k − i)
aN (i, k − i)
b(i, k − i)
−
(aN (i, k − i) − a(i, k − i)) .
aN (i, k − i)a(i, k − i)

Since the paths (Z1j , . . . , Znj ), 1 ≤ j ≤ N are i.i.d., the multivariate central limit
holds for the sequence of 2n(n − 1) dimensional vectors {(bN (i, k − i) − b(i, k −
i), aN (i, k − i) − a(i, k − i)), 0 ≤ i ≤ k ≤ n − 1}. The proof is concluded by tedious
computations using the Markov property, which we omit.
Remark 10. We can prove that the covariance matrix Γ0 is diagonal by a statistical argument. If we consider the tautological model {f (i, k − i), 0 ≤ i ≤ k ≤

25

n − 1}, i.e. θ = f and f0 is the true value. Then the likelihood is
LN (f ) =

n−1 X
k
X

aN (i, k − i){pN (i, k − i) log f (i, k − i)

k=0 i=0

+ qN (i, k − i) log f¯(i, k − i)} ,
where f¯(i, k−i) = 1−f (i, k−i). Thus we see that {pN (i, k − i), 0 ≤ i ≤ k ≤ n − 1}
is the maximum
likelihood estimator of f0 . This model is a regular statistical
√
model, thus N (pN − f0 ) converges weakly to the Gaussian distribution with covariance matrix In−1 (f0 ), where In (f ) is the Fisher information matrix of the
model. It is easily seen that In (f0 ) is the n(n − 1) dimensional diagonal matrix
with diagonal elements given by (14).

6.3

A general result for minimum contrast estimators

Theorems 2 and 3 are a consequence of the general result we prove in this section.
More precisely we demonstrate the consistency and the asymptotic normality of
a general estimator of which the MLE and the WLSE are particular cases.
For 0 ≤ i ≤ k ≤ n − 1, recall the definition of aN (i, k − i), pN (i, k − i) in (7)
and that f¯(θ, i, k − i) = 1 − f (θ, i, k − i). Let wN (i, k − i), 0 ≤ i ≤ k ≤ n − 1
be a sequence of random weights and let G be function defined on [0, 1] × (0, 1).
Define the empirical contrast function by
WN (θ) =

n−1 X
k
X

wN (i, k − i)G(pN (i, k − i), f (θ, i, k − i)) .

k=0 i=0

For instance, choosing G(p, q) = −p log q − (1 − p) log(1 − q) and wN (i, k − i) =
aN (i, k − i) yields
WN (θ) = −

n−1 X
k
X


aN (i, k − i) pN (i, k − i) log f (θ, i, k − i)

k=0 i=0

+ qN (i, k − i) log f¯(θ, i, k − i)
= − N −1 LN (θ) ,
so that minimizing WN is equivalent to maximizing the likelihood LN , defined
in (4). Choosing G(p, q) = (p − q)2 yields the weighted least squares contrast
function WN , defined in (8). We now define the minimum contrast estimator of θ0
by
W
θ̂N
= arg min WN (θ) .
θ∈Θ

26

W
In order to prove the consistency and asymptotic normality of θ̂N
, we make
the following assumptions on G and on the weights wN (i, k − i). Let ∂2 G and ∂22 G
denote the first and second derivatives of G with respect to its second argument.

Assumption 11. The function G is non negative, twice continuously differentiable
on [0, 1]×(0, 1) with G(p, q)−G(p, p) > 0 if p 6= q, ∂2 G(p, p) = 0 and ∂22 G(p, p) > 0.
Assumption 12. For all 0 ≤ i ≤ k ≤ n − 1, wN (i, k − i) converge almost surely
to w0 (i, k − i) and w0 (i, k − i) > 0.
W
Theorem 13. If Assumptions 1-(i), 1-(ii), 11 and 12 hold, then θ̂N
is consistent.
If moreover θ0 is an interior point of Θ and Assumption 1-(iii) holds, then
√
W
N (θ̂N
− θ0 ) converges weakly to a Gaussian distribution with zero mean.

The exact expression of the variance is given in the proof.
Proof. Under Assumption 11, the strong law of large numbers shows that WN (θ)
converges almost surely to
W(θ) =

k
n−1 X
X

w0 (i, k − i)G(f0 (i, k − i), f (θ, i, k − i)) .

k=0 i=0

Assumptions 1-(ii) and 11 ensure that θ0 is the unique minimum of W. Indeed,
G(p, q) > 0 if p 6= q and G(p, p) = 0. Thus, W is minimized by any value of θ
such that f (θ, i, k − i) = f (θ0 , i, k − i). By Assumption 1-(ii), this implies θ = θ0 .
Moreover the convergence of WN to W is uniform, since Θ is compact and the
function f is twice continuously differential with respect to θ, its first variable.
W
. For the sake of completeness, we give a brief
This yields the consistency of θ̂N
W
minimizes WN , we have
proof. Since θ0 minimizes W and θ̂N
W
0 ≤ W(θ̂N
) − W(θ0 )
W
W
W
) − WN (θ0 ) + WN (θ0 ) − W(θ0 )
= W(θ̂N
) − WN (θ̂N
) + WN (θ̂N
W
W
) − WN (θ̂N
) + WN (θ0 ) − W(θ0 ) ≤ 2 sup |WN (θ) − W(θ)| .
≤ W(θ̂N
θ∈Θ

Since θ0 is the unique minimizer of W, for  > 0, we can find δ such that if θ ∈ Θ
and kθ − θ0 k > , then W(θ) − W(θ0 ) ≥ δ. Thus
P(kθ̂N − θ0 k > ) ≤ P(W(θ̂N ) − W(θ0 ) ≥ δ)


≤ P 2 sup |WN (θ) − W(θ)| ≥ δ → 0 .
θ∈Θ

27

The central limit theorem is a consequence of the consistency and Lemma 9. A
first order Taylor extension of ẆN (θ) at θ0 yields
W
W
0 = ẆN (θ̂N
) = ẆN (θ0 ) + ẄN (θ̃N )(θ̂N
− θ0 ) ,
W
where θ̃N ∈ [θ0 , θ̂N
]. Setting f˙0 (i, k − i) = f˙(θ0 , i, k − i), we have

ẆN (θ0 ) =

n−1 X
k
X

wN (i, k − i)∂2 G(pN (i, k − i), f0 (i, k − i))f˙0 (i, k − i) .

k=0 i=0
2
G be the mixed second derivative of G. Note that
Let ∂12

∂2 G(f0 (i, k − i), f0 (i, k − i)) = 0 .
Thus, by the delta-method (see Dacunha-Castelle and Duflo,√
1986, Theorem 3.3.11)
and since wN converges almost surely to w0 , we obtain that N ẆN (θ0 ) converges
weakly towards
k
n−1 X
X

2
w0 (i, k − i)∂12
G(f0 (i, k − i), f0 (i, k − i))Λ0 (i, k − i)f˙0 (i, k − i) ,

k=0 i=0

where Λ0 (i, k) are independent Gaussian random
√ variables with zero mean and
variance γ0 (i, k) defined in 14. Equivalently, N ẆN (θ0 ) converges weakly to a
Gaussian vector with zero mean and covariance matrix H(θ0 ) defined by
H(θ0 ) =

n−1 X
k
X

2
w02 (i, k − i){∂12
G(f0 (i, k − i), f0 (i, k − i))}2

k=0 i=0

× γ0 (i, k)f˙0 (i, k − i)(f˙0 (i, k − i))0 .
By the law of large numbers, ẄN (θ) converges almost surely to Ẅ(θ) and this convergence is also locally uniform. Thus, ẄN (θ̃N ) converges almost surely to Ẅ(θ0 ).
Using again the fact that ∂2 G(p, p) = 0, we obtain
Ẅ(θ0 ) =

n−1 X
k
X

w0 (i, k − i)∂22 G(f0 (i, k − i), f0 (i, k − i))

k=0 i=0

× f˙0 (i, k − i)(f˙0 (i, k − i))0 .
Denote for brevity g(i, k − i) = w0 (i, k − i)∂22 G(f0 (i, k − i), f0 (i, k − i)). Then, for
any u ∈ Rd , we have
!2
n−1 X
k
d
X
X
uẄ(θ0 )u0 =
g(i, k − i)
us ∂s f (θ0 , i, k − i)
.
(15)
s=1

k=0 i=0

28

By assumption 12, g(i, k − i) > 0 for all 0 ≤ i ≤ k ≤P
n − 1, thus (15) is zero only
if for all k = 0, . . . , n − 1 and i = 0, . . . , k, we have ds=1 us ∂s f (θ0 , i, k) = 0. By
Assumption 1-(iii), this is possible only if us = 0 for all s = 1, . . . , d. Thus Ẅ(θ0 )
is positive definite.
We can now conclude that for large enough N , ẄN (θ̃N ) is invertible and we
can write
√
√
W
N (θ̂N
− θ0 ) = −Ẅ−1
(
θ̃
)
N ẆN (θ0 ) .
N
N
The right hand side converges weakly to the Gaussian distribution with zero mean
and covariance matrix Ẅ−1 (θ0 )H(θ0 )Ẅ−1 (θ0 ).

6.4

Proofs of theorems 2 and 3

Theorems 2 and 3 are a consequence of Theorem 13.
Lemma 14. Assumption 12 holds for the weights wN (i, k − i) = aN (i, k − i) and
−1
wN (i, k − i) = aN (i, k − i)p−1
N (i, k − i)qN (i, k − i), 0 ≤ i ≤ k ≤ n − 1.
Proof. For all 0 ≤ i ≤ k ≤ n − 1, aN (i, k − i) converges almost surely to P(Zk = i)
−1
−1
−1 ¯
and aN (i, k − i)p−1
N (i, k − i)qN (i, k − i) to P(Zk = i)f0 (i, k − i) f0 (i, k − i) .
Moreover Assumption 1-(i) implies that f0 (i, k − i) > 0 and f¯0 (i, k − i) > 0 for all
0 ≤ i ≤ k ≤ n − 1. Using Formula (13), this in turn implies that P(Zk = i) > 0
for all 0 ≤ i ≤ k ≤ n − 1.
Proof of Theorem 2. As mentioned above, the maximum likelihood estimator minimizes the contrast function W written with the function G(p, q) = −p log q −
(1 − p)log(1 − q) and the weights aN (i, k − i). Thus the proof of Theorem 2 consists in checking Assumption 11 and 12 to apply Theorem 13. Lemma 14 implies
that Assumption 12 holds.
The function G considered here satisfies Assumption 11. Indeed, for p, q ∈ (0, 1),
define
K(p, q) = G(p, q) − G(p, p) = p log(p/q) + (1 − p) log((1 − p)/(1 − q) .
Remark that K(p, q) is the Kullback-Leibler distance between the Bernoulli measures with respective success probabilities p and q. Then it is well known that
K(p, q) > 0 except if p = q. Indeed, by Jensen’s inequality,
K(p, q) ≥ − log(pq/p + (1 − p)(1 − q)/(1 − p)) = log 1 = 0 ,
and by strict concavity of the log function, equality holds only if p = q. Moreover,
∂2 G(p, q) = −p/q+(1−p)/(1−q) so ∂2 G(p, p) = 0 and ∂22 G(p, p) = p−1 (1 − p)−1 > 0.

29

Proof of Theorem 3. Again, the proof consists in checking Assumption 11 and 12
to apply Theorem 13. The latter holds by virtue of Lemma 14 and Assumption 11 trivially holds for the function G(p, q) = (p − q)2 . If wN (i, k − i) =
−1
p−1
N (i, k − i)qN (i, k − i)aN (i, k − i), then
H(θ0 ) = 2Ẅ(θ0 ) = 4In (θ0 )
=4

n−1 X
k
X
k=0 i=0

P(Zk = i)
f˙0 (i, k − i)(f˙0 (i, k − i))0 . (16)
f0 (i, k − i)f¯0 (i, k − i)

So that Σn (θ0 ) = Ẅ−1 (θ0 )H(θ0 )Ẅ−1 (θ0 ) = In−1 (θ0 ).
If the weights are chosen as wN (i, k) = aN (i, k), then w0 (i, k − i) = P(Zk = i)
and
H(θ0 ) = 4
Ẅ (θ0 ) = 2

k
n−1 X
X
k=0 i=0
n−1
k
XX

P(Zk = i)f0 (i, k − i)f¯0 (i, k − i)f˙0 (i, k − i)(f˙0 (i, k − i))0 , (17)
P(Zk = i)f˙0 (i, k − i)(f˙0 (i, k − i))0 .

(18)

k=0 i=0

6.5

Proofs of Corollaries 6 and 7

Corollaries 6 and 7 are a consequence of Theorem 13. The assumptions on the
weights wN and on the functions G have been already verified in the previous
section. We have to prove the Assumption 1 on the choice function f defined
in (10). Hypothesis 1-(i) is obvious.
By elementary computations, we have, for 0 ≤ i ≤ k ≤ n − 1,

α 
α0
c+i
c0 + i
f (α, c, i, k − i) = f (α0 , c0 , i, k − i) ⇔
=
c+k−i
c0 + k − i
log(c0 + i) − log(c0 + k − i)
α
=
. (19)
⇔
α0
log(c + i) − log(c + k − i)
Plugging the pairs (i, k) = (0, 1) and (i, k) = (0, 2) into (19) yields
log(c0 ) − log(c0 + 1)
log(c0 ) − log(c0 + 2)
=
,
log(c) − log(c + 1)
log(c) − log(c + 2)
or equivalently
log(1 + 1/c0 )
log(1 + 1/c)
=
.
log(1 + 2/c0 )
log(1 + 2/c)
30

(20)

It is easily checked that the function x → log(1 + x)/ log(1/2x) is strictly increasing on (0, ∞). Thus (20) implies that c = c0 . Plugging this equality into (19)
yields α = α0 . This proves Assumption 1-(ii).
We now prove that if n ≥ 2, the vectors {∂α f (θ0 , i, k − i), 0 ≤ i ≤ k ≤ n − 1}
and {∂c f (θ0 , i, k − i), 0 ≤ i ≤ k ≤ n − 1} are linearly independent in Rn(n−1) .
For 0 ≤ i ≤ k ≤ n − 1, we have,


c+i
,
∂α f (α, c, i, k − i) = f (α, c, i, k − i)f (α, c, k − i, i) log
c+k−i
α(k − 2i)
∂c f (α, c, i, k − i) = f (α, c, i, k − i)f (α, c, k − i, i)
.
(c + i)(c + k − i)
Let (u, v) ∈ R2 and assume that for all i, j ≤ n − 1 such that i + j ≤ n − 1, it
holds that


c0 + i
α0 (j − i)
=0.
u log
+v
c0 + j
(c0 + i)(c0 + j)
Replacing (i, j) for instance successively by (0, 1) and (0, 2) yields



c
α0

0

+v
=0,
u log
c0 (c0 + 1)
 c0 + 1 
c0
2α0


=0.
+v
u log
c0 + 2
c0 (c0 + 2)
If (u, v) 6= (0, 0), this implies




c0 + 2
c0 + 2
c0 + 1
c0 + 1
log
+2
log
=0.
c0
c0
c0
c0
By strict convexity of the function x → x log x on (0, ∞), this is impossible.
Thus u = v = 0 and Assumption 1-(iii) holds.

6.6

Proof of Theorem 8

Proof of Theorem 8, case α0 = 1. In this case the model is Pólya’s urn, and we
have
 




n−1
X
1
1
4
1
E
+E
−
.
(21)
In (c) =
2c + k
c + Zk
c + k − Zk
2c + k
k=0
The distribution of Zk is given by
 
k c(c + 1) · · · (c + i − 1) × c(c + 1) · · · (c + k − i − 1)
P(Zk = i) =
.
2c(2c + 1) · · · (2c + k − 1)
i
31

Thus,
 X
k  
k c(c + 1) · · · (c + i − 1) × c(c + 1) · · · (c + k − i − 1) 1
1
=
.
E
c + Zk
i
2c(2c + 1) · · · (2c + k − 1)
c+i
i=0


For any c > 0, there exists constants C1 < C2 such that, for all integers h ≥ 1,
h
Y
C1 h ≤
(1 + c/i) ≤ C2 hc .
c

i=1

Therefore, there exists a constant C > 0 such that for all k ≥ 1,

−1







k−1
c−2
c−1
O(k ) if c > 1 ,
X i
i
1
= O(k −1 log k) if c = 1 ,
≤ Ck −2
1−
E

c + Zk
k
k

i=1
O(k −c ) if c < 1 .
In all three cases, we obtain that the first series in (21) is summable. By symmetry,
the sum of the second expectations is also finite.
Proof of Theorem 8, case α0 < 1. In this case, we know by Theorem 4 that Zn /n
converges almost surely to 1/2. This implies that f (θ, Zn , n−Zn ) converges almost
surely to 1/2 for all θ. By Cesaro’s Lemma, this implies that n−1 `n (θ) → − log 2
a.s.
Proof of Theorem 8, case α0 > 1. Let Ω1 be the event that color 1 is eventually selected, which happens with probability 1/2 by Theorem 4. Then, on Ω1 , Zn /n → 1
and if k > T∞ , then Xk+1 = 1 and Zk = k − Q∞ . Thus for large enough n, the
log-likelihood on one path becomes
`n (θ) =

T∞
X

Xk+1 log f (θ, Zk , k − Zk ) + (1 − Xk+1 ) log{1 − f (θ, Zk , k − Zk )}

k=0

+

n
X

log f (θ, k − Q∞ , Q∞ ) .

k=T∞ +1

As k → ∞, for any α > 0,
(c + Q∞ )α
log f (θ, k − Q∞ , Q∞ ) = − log 1 +
(c + k − Q∞ )α


32


∼−

(c + Q∞ )α
.
(c + k − Q∞ )α

If α ≤ 1 the series is divergent and thus limn→∞ `n (θ) = −∞. If α > 1 then the
series is convergent and thus, on Ω1 ,
lim `n (θ) =

n→∞

=
+

∞
X
k=0
T∞
X

Xk+1 log f (θ, Zk , k − Zk ) + (1 − Xk+1 ) log{1 − f (θ, Zk , k − Zk )}
Xk+1 log f (θ, Zk , k − Zk ) + (1 − Xk+1 ) log{1 − f (θ, Zk , k − Zk )}

k=0
∞
X

log f (θ, k − Q∞ , Q∞ ) .

T∞ +1

This implies that arg maxθ∈Θ `n (θ) = arg maxθ∈Θ,α>1 `n (θ) and that this argmax
is a random variable which is a function of the whole path, and does not depend
on the true value θ0 .
Acknowledgment We thanks Guy Theraulaz and Hugues Chaté for providing
their material framework and their field expertise to allow the first author to
collect the data of the Argentine ants experiments. These experiments are part
of the project TRACES supported by the CNRS. They were done during two
visits in April and July 2012 of the first author to the Centre de Recherches sur
la Cognition Animale (CRCA, Centre de Recherches sur la Cognition Animale,
UMR 5169, Paul Sabatier University, Toulouse), whose hospitality is gratefully
acknowledged.

References
Arganda, S., S. Nicolis, A. Perochain, C. Péchabadens, G. Latil, and A. Dussutour (2014): “Collective choice in ants: The role of protein and carbohydrates
ratios,” J. Insect Physiol., 69, 19–26.
Aron, S., J.-L. Deneubourg, S. Goss, and J. Pasteels (1990): Functional Selforganisation illustrated by Inter-nest Traffic in Ants : the Case of the Argentine
Ant, Springer-Verlag, chapter Biological Motion, 533–547.
Aron, S., J. Pasteels, and J.-L. Deneubourg (1989): “Trail-laying behaviour during
exploratory recruitment in the argentine ant, iridomyrmex humilis (mayr),”
Biol. of Behav., 14, 207–217.
Beckers, R., J.-L. Deneubourg, and S. Goss (1992): “Trails and u-turns in the
selection of a path by the ant lasius niger,” J. Theor. Biol., 159, 397–415.

33

Beckers, R., J.-L. Deneubourg, and S. Goss (1993): “Modulation of trail laying
in the ant lasius niger (hymenoptera: formicidae) and its role in the collective
selection of a food source,” J. of Insect Behav., 6, 751–759.
Dacunha-Castelle, D. and M. Duflo (1986): Probability and statistics. Vol. II, New
York: Springer-Verlag.
Davis, B. (1990): “Reinforced random walk,” Probab. Theory Related Fields, 84,
203–229.
Deneubourg, J.-L., S. Aron, S. Goss, and J. Pasteels (1990): “The self-organizing
exploratory pattern of the argentine ant,” J. of Insect Behav., 3, 159–168.
Dussutour, A., J.-L. Deneubourg, and V. Fourcassié (2005): “Amplification of
individual preferences in a social context: the case of wall-following in ants,”
Proc. R. Soc. London, Ser. B, 272, 705–714.
Freedman, D. (1965): “Bernard friedman’s urn,” Ann. Math. Stat., 36, 956–970.
Garnier, S., A. Guérécheau, M. Combe, V. Fourcassié, and G. Theraulaz (2009):
“Path selection and foraging efficiency in argentine ant transport networks,”
Behav. Ecol. Sociobiol., 63, 1167–1179.
Jeanson, R., F. Ratnieks, and J.-L. Deneubourg (2003): “Pheromone trail decay
rates on different substrates in the pharaohs ant, monomorium pharaonis,”
Physiol. Entomol., 28, 192–198.
Khanin, K. and R. Khanin (2001): “A probabilistic model for the establishment
of neuron polarity,” J. Math. Biol., 42, 26–40.
Le Goff, L. and O. Raimond (2015): “Vertex reinforced non-backtracking random
walks: an example of path formation,” URL http://arxiv.org/abs/1506.
01239, arXiv:1506.01239.
Nicolis, S. and J.-L. Deneubourg (1999): “Emerging patterns and food recruitment in ants: an analytical study,” J. Theor. Biol., 198, 575–592.
Nicolis, S. and A. Dussutour (2008): “Self-organization, collective decision making
and resource exploitation strategies in social insects,” Eur. Phys. J. B, 65, 379–
385.
Pasteels, J., J.-L. Deneubourg, and S. Goss (1987): “Transmission and amplification of information in a changing environment: The case of insect societies,”
Eds I. Prigogine & M. Sanglier. Gordes, Bruxelles.

34

Pemantle, R. (2007): “A survey of random processes with reinforcement,” Probability Surveys, 4, 1–79.
Pólya, G. (1931): “Sur quelques points de la théorie des probabilités,” Ann.
I.H.P., 1, 117–161.
Tarrès, P. (2011): “Localization of reinforced random walks,” URL http://
arxiv.org/abs/1103.5536, arXiv:1103.5536.
Thienen, W., D. Metzler, D.-H. Choe, and V. Witte (2014): “Pheromone communication in ants: a detailed analysis of concentration-dependent decisions in
three species,” Behav. Ecol. Sociobiol., 68, 1611–1627.
Van Vorhis Key, S. and T. Baker (1982): “Trail-following responses of the argentine ant, iridomyrmex humilis (mayr), to a synthetic trail pheromone component
and analogs,” J. Chem. Ecol., 8, 3–14.
Vittori, K., G. Talbot, J. Gautrais, V. Fourcassié, A. Araujo, and G. Theraulaz
(2006): “Path efficiency of ant foraging trails in an artificial network,” J. Theor.
Biol., 239, 507–515.
Wasserman, L. (2004): All of Statistics: A Concise Course in Statistical Inference,
Springer.

35

