1

A Bayesian framework for distributed estimation
of arrival rates in asynchronous networks

arXiv:1702.04939v1 [math.OC] 16 Feb 2017

Angelo Coluccia, member, IEEE and Giuseppe Notarstefano, member, IEEE

Abstract—In this paper we consider a network of agents monitoring a spatially distributed arrival process. Each node measures
the number of arrivals seen at its monitoring point in a given
time-interval with the objective of estimating the unknown local
arrival rate. We propose an asynchronous distributed approach
based on a Bayesian model with unknown hyperparameter,
where each node computes the minimum mean square error
(MMSE) estimator of its local arrival rate in a distributed
way. As a result, the estimation at each node “optimally” fuses
the information from the whole network through a distributed
optimization algorithm. Moreover, we propose an ad-hoc distributed estimator, based on a consensus algorithm for timevarying and directed graphs, which exhibits reduced complexity
and exponential convergence. We analyze the performance of
the proposed distributed estimators, showing that they: (i) are
reliable even in presence of limited local data, and (ii) improve
the estimation accuracy compared to the purely decentralized
setup. Finally, we provide a statistical characterization of the
proposed estimators. In particular, for the ad-hoc estimator, we
show that as the number of nodes goes to infinity its mean
square error converges to the optimal one. Numerical Monte
Carlo simulations confirm the theoretical characterization and
highlight the appealing performances of the estimators.
Index Terms—distributed estimation, Empirical Bayes, pushsum consensus, cyber-physical systems.

I. I NTRODUCTION
Arrival processes provide a useful description for events
occurring with some probability in a given time or space interval. Applications range from communications and transports
to medicine (e.g., diagnostic imaging) and astronomy (e.g.,
particle detection), [2]. From a statistical point of view the
prominent model for arrival processes is notoriously Poisson.
Indeed, even when independence between arrivals cannot
be assumed, the superposition of a large number of nonPoisson processes is approximately distributed as a Poisson
(Palm-Khintchine Theorem), [3]. Also, the limit distribution
of counting processes described by the Binomial distribution
is a Poisson, according to the law of rare events.
In modern network contexts (as, e.g., data, communication
and sensor networks or Intelligent Transportation Systems)
estimating the process intensity at different locations, i.e., the
Copyright (c) 2015 IEEE. Personal use of this material is permitted.
However, permission to use this material for any other purposes must be
obtained from the IEEE by sending a request to pubs-permissions@ieee.org.
Angelo Coluccia and Giuseppe Notarstefano are with the Department of
Engineering, Università del Salento, via Monteroni, 73100, Lecce, Italy,
{name.lastname}@unisalento.it.
A preliminary short version of this paper has appeared as [1]. This result is
part of a project that has received funding from the European Research Council
(ERC) under the European Unions Horizon 2020 research and innovation
programme (grant agreement No 638992 - OPT4SMART).

arrival rates at different nodes, is an important preliminary
problem to be solved in order to gain context awareness.
Clearly, each arrival rate can be estimated by performing a
decentralized Maximum Likelihood (ML) estimation based
only on the local arrivals at the given node. In this paper we
want to investigate how the estimation can be improved by
exploiting the cooperation among the nodes. Recently, a great
interest has been devoted to cooperation schemes in which
estimation is performed by a network of computing nodes in
a distributed way, rather than by collecting all the data in a
central unit, [4]–[6]. How to use the information from other
nodes in the network and how to design a distributed algorithm
merging such information will be the focus of the paper.
Distributed estimation has received a widespread attention
in the distributed computation literature, especially as a natural
application of linear (average) consensus algorithms, see, e.g.
[7], or the recent surveys [8], [9]. Nodes typically interact
iteratively with their neighbors by means of a “diffusionlike” process in which the estimation is improved by suitably
combining the estimations from neighboring nodes, [6], [10],
[11]. An incremental and a diffusive distributed algorithm
with finite time convergence are proposed in [12] for (static)
state estimation. Distributed optimization is also strictly related
to distributed estimation. In [5] a distributed Alternating Direction Method of Multipliers (ADMM) has been introduced
as a tool for distributed ML estimation of vector parameters
in a wireless sensor network. Notice that ADMM and other
distributed optimization algorithms, as, e.g., [13], can be used
as building blocks to solve parts of an estimation problem in
a distributed set-up. ML approaches for distributed estimation
of a commonly observed parameter have been proposed also
in [14]. In [15], [16] consensus-based algorithms have been
developed to simultaneously estimate a common parameter
measured by noisy sensors and classify sensor types. In [17]
consensus-based algorithms have been developed to estimate
global parameters in a linear Bayesian framework. In [18]
a distributed algorithm is proposed for adaptive Bayesian
estimation of a common parameter with known prior. Dynamic
methods have been also proposed in which the nodes keep
collecting new measurements while interacting with each
other. In [19] a diffusion-based Recursive Least Squares (RLS)
algorithm is proposed to estimate a constant parameter, but
with dynamically acquired measurements. Finally, in [20] and
[21] distributed ADMM-based algorithms are proposed for the
estimation of random signals and dynamical processes.
Differently from the above references, in our work we
consider a more general Bayesian framework that allows the
nodes to improve their local estimate, rather than reaching a

2

consensus on a common parameter. As detailed later, we will
consider a special model for continuous mixtures of Poisson
variables. Poisson mixtures have been widely used in the
arrival-rate estimation literature to model non-homogeneous
scenarios, see, e.g., [22], [23] as early references. The survey
[24] provides an extensive review of properties and applications for Poisson mixtures. In a Bayesian setting the classical
mixing (prior) distribution is the Gamma [24], [25], which
results in a closed-form posterior [26].
The main contribution of the paper is twofold. First, we
propose a distributed estimation scheme for arrival rates in
an asynchronous network, based on a hierarchical probabilistic framework. Specifically, we develop an Empirical Bayes
approach, in which the arrival rates are treated as random
variables, whose prior distribution is parametrized by an unknown hyperparameter to be determined via ML estimation. In
particular, we borrow from the centralized statistics literature
the classical Gamma-Poisson model (for Poisson mixtures),
assuming however the hyperparameter is unknown, and extend
it to a scenario with non-homogeneous sample sizes. We show
that the ML estimation of the hyperparameter is a separable
optimization problem, that can be solved in a distributed way
over the network by using a distributed optimization algorithm.
Thanks to this modeling idea, the local estimates are obtained
by taking advantage of the whole network data, thus improving
the accuracy especially when the amount of local data is
scarce. With this approach we are able to capture the fact that
arrival rates are not the outcomes of isolated phenomena, but
rather the expression of global properties of the process. For
this “optimal” estimator we characterize mean and variance at
steady-state (i.e., once agents have reached consensus on the
optimal solution) for networks with a large number of agents.
Second, we propose an alternative ad-hoc distributed estimator that, although suboptimal, performs comparably to the
optimal one. The main advantage of this ad-hoc estimator is
that the resulting algorithm has a simple update rule based on
linear consensus protocols, thus exhibiting the same appealing
exponential convergence. For the ad-hoc estimator we also
characterize transient (at any time-instant) mean and variance
for a given number of agents. Notably, we show that at steadystate and for large number of agents the ad-hoc estimator
attains the performance of the optimal one.
To strengthen these two contributions we perform a Monte
Carlo analysis and compare the theoretical expressions obtained for mean and variance with their sample counterparts.
The numerical computations confirm the theoretical analysis
and show that some key assumptions made for a rigorous, but
tractable analysis are not restrictive. Moreover, they highlight
other interesting features of the proposed distributed estimators. For example, the ad-hoc estimator achieves performances
close to the optimal one even for a limited number of nodes.
The paper is organized as follows. In Section II we introduce
the model of a monitoring network and set up the estimation
problem. In Section III we develop the proposed distributed
estimators based on the Empirical Bayes approach. The statistical performances of the two estimators are analyzed in
Section IV, while in Section V we perform a Monte Carlo
analysis to confirm the theoretical results.

II. M ONITORING NETWORK AND ARRIVAL RATE
ESTIMATION PROBLEM

We consider a network of monitors having sensing, communication and computation capabilities. That is, each monitor
can measure the number of arrivals in a given measurement
time-scale (e.g., 1 second or 1 minute), share local data with
neighboring agents, and perform local computations on its own
and its neighbors’ data. The objective is to fuse the data in
order to improve the estimation of the arrival rates.
Formally, each node collects measurements asynchronously
in an observation window, over which the underlying process
can be assumed to be stationary, i.e., the set of rates can be
considered constant over the observation window. Clearly, the
arrival rates have stationary increments, so that the number
of arrivals in disjoint intervals are statistically independent.
A measurement consists of the number of arrivals detected
at a given location in the (common) time-scale interval.
Accordingly, for each monitor i ∈ {1, . . . , N } we introduce
the following variables:
• λi unknown arrival rate;
• yi,` the `-th collected measurement (number of arrivals
per time-scale interval);
• ni ∈ [1, nmax ] number of measurements yi,` collected
in the observation window (where nmax is the maximum
number of measurements that can be collected).
The conditional distribution of yi,` given λi is a Poisson
random variable with parameter λi , i.e., yi,` |λi ∼ Poisson(λi ).
All measurements are assumed to be independent.
def PN
We denote by n = i=1 ni the total number of measurements. If all the nodes have the same number of measurements,
i.e., ni = n/N for all i ∈ {1, . . . , N } we say that the network
is homogeneous.
In Figure 1 a scheme of the network measurement scenario
is depicted with the variables of interest.
observation window
y1,1

y1,2

i=1 (λ1)
n1=2

i=2 (λ2)
n2=4

i=N (λN)
nN=2

time-scale
interval
y2,1

yN,1

y2,2

y2,3

.
.
.

y2,4

yN,2

Figure 1. Scheme of the measurement scenario.

We assume that the network evolution is triggered by
a universal slotted time, t ∈ Z≥0 , not necessarily known
by the monitors. The monitors communicate according to a
time-dependent directed communication graph t 7→ G(t) =
({1, . . . , N }, E(t)), where {1, . . . , N } are the monitor identifiers and the edge set E(t) describes the communication
among monitors: (i, k) ∈ E(t) if monitor i communicates
to k at time t ∈ Z≥0 . For each node i, the nodes sending
information to i at time t, i.e., the set of k ∈ {1, . . . , N } such

3

that (k, i) ∈ E(t) is the set of in-neighbors of i at time t,
and is denoted by NiI (t). We make the following minimal
assumption on the graph connectivity. First, we recall that
a fixed directed graph is strongly connected if for any pair
of nodes, i and j, there exists a directed path (i.e., a set of
consecutive edges) from i to j.

λ1
λ2

λN

λi ~ Gamma(a,b)

Poiss(λ2)
statistical relationships

Poiss(λ1)

Assumption II.1 (Uniform joint strong connectivity).
an integer Q ≥ 1 such that the graph
 There exists
S(t+1)Q−1
{1, . . . , N }, τ =tQ
E(τ ) is strongly connected ∀ t ≥ 0.
It is worth remarking that this network setup is very general,
since it naturally embeds asynchronous scenarios as well as
missing measurements due to, e.g., sensor failures.

communication
links

Poiss(λN)

Figure 2. Scheme of the Bayesian model

A. Empirical Bayes approach in monitoring networks
III. D ISTRIBUTED ARRIVAL RATE ESTIMATION
VIA E MPIRICAL BAYES
In a decentralized set-up, in which nodes do not communicate, each node could estimate λi based on the sample
{yi,` }`=1,...,ni by simply computing the empirical mean of the
available measurements. That is, the decentralized estimator is
λ̂dec
i =

ni
σi
1 X
yi,` =
ni
ni
`=1

def Pni
where σi = `=1
yi,` .
Notice that, the decentralized estimator turns out to be the
ML estimator of λi when node i can use only its own data.
However, decentralized estimation yields reliable estimates
only when the number of samples ni is large enough.
In our heterogenous set-up, it may happen that some nodes
satisfy such a condition, while other ones do not have enough
data, thus resulting in a poor estimation. In this paper we
propose a distributed estimation scheme in which every node,
especially the ones with fewer measurements, take advantage
from cooperating with neighboring nodes.
How the measurements at other nodes can help the local
estimation at a given node is a nontrivial issue and needs to
be investigated by means of a suitable probabilistic framework.
Specifically, we adopt a Bayesian model in which all the
unknown arrival rates λi s are i.i.d. random variables ruled
by a common probability distribution that captures the spatial
variability of the process intensity.
This model belongs to the family of Poisson mixtures and
is broadly accepted in the (centralized) statistical literature for
modeling non-homogeneous scenarios, see, e.g., [22], [23],
and [24, and references therein] for a through and extensive
survey of properties and applications. As customary, see,
e.g., [26], we adopt as common distribution the conjugate
prior of the Poisson, which is the Gamma distribution. This
choice allows one to obtain a close-form expression for the
posterior (predictive) distribution. However, notice that, our
model extends the classical Gamma-Poisson mixture since the
sample sizes ni , i ∈ {1, . . . , N }, can be different and thus
exhibits even more flexibility.
A scheme of the proposed Bayesian framework for our
network scenario is depicted in Figure 2.

In applying a Bayesian estimation approach to a network
context, the assumption that the prior distribution is fully
known to all monitors is rather strong and may be a severe
limitation in realistic scenarios. To overcome this limitation
we adopt the Empirical Bayes approach in which only the
class of the prior is known, i.e., λi ∼ Gamma(a, b), where
the shape parameter a is known, but the scale parameter b is
unknown. The assumption that a is known, while only b is
unknown, says that only the shape of the Gamma distribution
(determined by the parameter a) is known, while the scaling is
not. This assumption is reasonable in many applications, since
it is a way to embed a rough information on the phenomenon,
and is customary for the sake of mathematical tractability, [26].
The hyperparameter b can be estimated via a ML procedure.
To this aim, we need the joint distribution of all measurements
{yi,` }`=1,...,ni for each agent i. The likelihood function is the
product of the marginal distributions of all agents
L(y1 , . . . , yN |b) =

N
Y

p(yi |b)

(1)

i=1

where yi = [yi,1 · · · yi,ni ]T . The marginal distribution of
agent i is derived from the joint distribution of yi and λi ,
!
Z ∞ Y
ni
p(yi |b) =
f (yi,` |λi , b) p(λi |b)dλi
0

`=1
∞

λσi e−ni λi λa−1
e−λi /b
i
Qi ni
=
dλi
Γ(a)ba
0
`=1 yi,` !

σi +a
Γ(σi + a)
b
Q
=
ni
Γ(a)ba `=1
yi,` ! ni b + 1
Z

(2)

By using eq. (2) into eq. (1) the likelihood is rewritten as:
σi +a
N 
1 Y
b
L(y1 , . . . , yN |b) ∝ aN
.
(3)
b i=1 ni b + 1
Thus, the ML estimator b̂ML of b can be found by solving the
following optimization problem
(
)

N
X
b
ML
b̂ = arg min aN log b −
(σi + a) log
ni b + 1
b∈R+
i=1
= arg min
b∈R+

N
X
i=1

f (b; ni , σi ).

(4)

4

The problem can be solved in closed-form only for the
homogeneous case where all ni s are equal, i.e., for ni = n/N ,
recalling that n is the total number of measurements. In
this case the ML estimator of b based on the entire set of
measurements is given by
N

b̂

hom

σ
1 X
σi =
=
an i=1
an

def PN
def PN
where n = i=1 ni , and σ = i=1 σi .
After obtaining an estimate for b, the Empirical Bayes
estimator of the arrival rate λi that minimizes the Mean Square
Error (MMSE) can be obtained by computing the conditional
mean of the posterior distribution p(λi |yi , b). The latter is
given by the ratio between the joint pdf p(yi , λi |b) and the
marginal pdf p(yi |b) as from eq. (2), i.e.,
ni b+1 
−σi −a
b
λσi i +a−1 e−λi b
(5)
p(λi |yi , b) =
Γ(σi + a)
ni b + 1

b
Eq. (5) is a Gamma pdf with parameters (σi +a, ni b+1
), hence
the Empirical Bayes MMSE estimator of each λi is
ML
λ̂EB
]=
i = E[λi |yi , b̂

b̂ML
b̂ML ni + 1

(a + σi ).

(6)

Remark III.1. It is worth highlighting that the Bayesian
estimate is especially useful when ni is small. In fact, nodes
improve the quality of their local estimate by combining the
frequentist estimation (based only on local observations) with
a correction term (based on a prior global knowledge), which
is estimated in a cooperative way. Indeed, we can rewrite
b̂ML ni
λ̂EB
= ρ nσii + (1 − ρ)ab̂ML , with ρ = b̂ML
. When the
i
ni +1
local information is abundant (ni → ∞ and thus ρ → 1), the
MMSE estimator (6) tends towards λ̂dec
i , meaning that when
ni is large no further information can be inferred from the
network. Conversely, when local information is scarce, i.e.,
the sample size ni is small, even one, the MMSE estimator
(6) approaches the estimate of the global mean E[λi ] = ab.
In the following we will also consider an ad-hoc estimator
obtained by using b̂hom instead of the optimal b̂ML , i.e.,
def

λ̂ad-hoc
=
i

b̂hom
b̂hom

ni + 1

(a + σi ) =

σ
(a + σi ).
an + σni

Clearly, the estimator λ̂EB
follows the “rationale” of the
i
Empirical Bayes approach and, therefore, has performance
guarantees inherited from the ML procedure. Conversely, the
ad-hoc estimator is an alternative that at the moment has the
only advantage of having a closed-form expression, whose
performances need to be understood. In the rest of the paper
we will show that, not only this estimator leads to a simpler
and faster distributed algorithm, but also that for large number
of agents performs as λ̂EB
i .

separable cost (i.e., the sum of N local costs), hence it can be
solved by using available distributed optimization algorithms
for asynchronous networks [13], [27]. We propose a distributed
estimator in which each node implements the local update rule
of the chosen distributed optimization algorithm.
The ML estimation problem (4) is not guaranteed to be
convex in general. This is quite common in the estimation
literature. However, since the function is coercive, there exists
(at least) a minimizer and, thus, it is reasonable to apply
descent algorithms, as the one in [13], which, under suitable
conditions, guarantee convergence to a local minimizer.
The Empirical Bayes distributed estimator is as follows.
At each t ∈ Z≥0 , each agent i stores a local state ξi (t), an
EB
ML
and an estimate λ̂EB
estimate b̂ML
i (t) of λ̂i . The
i (t) of b̂
node initializes its local state ξi to an initial value ξi0 chosen
according to the distributed optimization algorithm in use, and
sets b̂ML
i (0) = σi /(ani ) (which would be the solution of (4)
if i were the only agent). Then it updates its estimate of b̂ML
by using the local update rule of the chosen distributed optimization algorithm, and updates the current estimate λ̂EB
i (t) by
using (6). The algorithm is defined formally in the following
table. For each t ∈ Z≥0 , let {ξˆk (t)}k∈NiI (t) be the collection
of states of the in-neighbors of node i, opt_local the local
update of the chosen distributed optimization algorithm, and
γ(t) an algorithm parameter as, e.g., a time-varying step-size.
Algorithm 1 Empirical Bayes distributed estimator
Initialization: ξi (0) = ξi0 , b̂ML
i (0) =
Iterate:

ˆ
(b̂ML
i (t + 1), ξi (t + 1)) =


ˆi (t), {ξˆk (t)}
I
opt_local b̂ML
(t),
ξ
;
γ(t)
,
k∈Ni (t)
i
λ̂EB
i (t + 1) =

From eq. (6) it is clear that each agent can compute the
Empirical Bayes MMSE estimator provided it knows b̂ML .
Optimization problem (4), giving the ML estimator of b, has a

b̂ML
i (t + 1)
b̂ML
i (t + 1)ni + 1

(a + σi ).

As an example, we show the opt_local function for the
distributed subgradient-push method proposed in [13]. To be
consistent with the notation in [13] we let ξi = (vi , yi , xi ),
which is initialized to (vi (0), yi (0), xi (0)) = (1, 1, xi0 ), with
xi0 an arbitrary initial value. Also, dk (t) denotes the number
of out-neighbors of node k at time t.
function ξi (t + 1) =

ˆ
ˆ
opt_local b̂ML
i (t), ξi (t), {ξk (t)}k∈N I (t) ; γ(t)
i

vi (t + 1) =

X
k∈NiI (t)∪{i}

yi (t + 1) =

X
k∈NiI (t)∪{i}

b̂ML
i (t + 1) =

B. Distributed estimators

σi
ani .

xk (t)
dk (t)
yk (t)
dk (t)

vi (t + 1)
yi (t + 1)

xi (t + 1) = vi (t + 1) − γ(t + 1)∇f (b̂ML
i (t + 1); ni , σi )

with ∇f (b; ni , σi ) = b(niab+1) (ani b − σi ).
If problem (4) has a unique minimizer, the distributed opti-

5

mization algorithm guarantees that all nodes reach consensus
on the global minimizer. That is,
ML
lim b̂ML
,
i (t) = b̂

for all i ∈ {1, . . . , N }.

t→∞

From the convergence properties of the chosen distributed optimization algorithm it follows immediately that the proposed
distributed estimator asymptotically computes at each node i
the Empirical Bayes MMSE estimator of λi .
However, most of the available distributed optimization
algorithms, as the ones in [13], [27], need the tuning of a
global parameter (we denoted it γ), and typically exhibit a
sub-exponential convergence even in static graphs. To overcome these drawbacks, we propose an alternative distributed
estimator with reduced complexity that, although suboptimal,
will be shown to perform comparably to the optimal one.
The ad-hoc distributed estimator is defined as follows.
For each t ∈ Z≥0 , each node i ∈ {1, . . . , N } stores in memory
hom
two local states si (t) and ηi (t), an estimate b̂hom
,
i (t) of b̂
ad-hoc
ad-hoc
and an estimate λ̂i
(t) of λ̂i
. Let wik (t) ∈ R≥0 be a
set of weights such that wik (t) > 0 if (i, k) ∈ E(t) or k = i,
and wik (t) = 0 otherwise. The ad-hoc distributed estimator is
given in the following table.

For the algorithm to converge we need the following assumption together with the Assumption II.1 (uniform joint
connectivity of the communication digraph).
Assumption III.2 (Properties of W (t)). For
Pneach t ∈ Z≥0 ,
the matrix W (t) is column stochastic, i.e., i=1 wik (t) = 1,
and there exists a positive constant α > 0 such that wii (t) > α
and wik (t) ∈ {0} ∪ [α, 1].
Remark III.3. It is worth noting that the column stochasticity
assumption above is not the usual assumption used in linear
consensus algorithms in which row stochasticity is assumed.
Lemma III.4. Let t 7→ G(t) be a uniformly jointly strongly
connected graph (Assumption II.1) and {W (t)}t≥0 a sequence of matrices satisfying Assumption III.2. Then denoting
Φ(t) = W (t)W (t − 1) · · · W (0) with (i, k) element φik (t),
the following holds true.
(i) The matrix sequence {W (t)}t≥0 is weakly ergodic, i.e.,
lim |φik (t) − φih (t)| = 0

t→∞

for all i, k, h ∈ {1, . . . , N }.
(ii) There exists some µ > 0 such that for all t ≥ 0,
1T Φ(t) ≥ µ > 0, i.e., for any i ∈ {1, . . . , N },
N
X

Algorithm 2 ad-hoc distributed estimator
Initialization: si (0) = σi , ηi (0) = ni , b̂hom
i (0) =
b̂hom
i (0)
σi /(ani ), λ̂ad-hoc
(0)
=
(a
+
σ
i ).
i
b̂hom
i (0)ni +1
Iterate:
X
si (t + 1) =
wik (t)sk (t)

The result is well-known and can be found, e.g., in [13].
Further references on this result under the same or different
connectivity assumptions are [28]–[31].

k∈NiI (t)∪{i}

X

ηi (t + 1) =

Proposition III.5. Let Assumption II.1 hold. Then the ad-hoc
distributed estimator (Algorithm 2) satisfies

wik (t)ηk (t)

lim b̂hom
i (t)
t→∞
lim λ̂ad-hoc
(t)
i
t→∞

k∈NiI (t)∪{i}

b̂hom
i (t + 1) =
λ̂ad-hoc
(t
i

+ 1) =

1 si (t + 1)
a ηi (t + 1)
b̂hom
i (t + 1)
b̂hom
i (t + 1)ni + 1

(a + σi ).

s(t + 1) = W (t)s(t)

def

T

def

with s(0) = [σ1 . . . σN ] = σ, η(0) = [n1 . . . nN ] = n and
W (t) the matrix with elements wij (t). Let us denote
Φ(t) = W (t − 1) · · · W (0)

(8)

the state transition matrix associated to each one of the linear
systems (7), so that
s(t) = Φ(t)s(0)
η(t) = Φ(t)η(0).

si (t) → φ̄i (t)

= λ̂ad-hoc
.
i

(9)

N
X

sk (0) = φ̄i (t)

k=1

and
ηi (t) → φ̄i (t)

(7)

η(t + 1) = W (t)η(t)

= b̂hom

Proof. From Lemma III.4, {W (t)}t≥0 is weakly ergodic so
that there exists a sequence of stochastic vectors φ̄(t) =
[φ̄1 (t) . . . φ̄N (t)] such that for each i ∈ {1, . . . , N }, φik (t) →
φ̄i (t) for all k ∈ {1, . . . , N }. From eq. (9) it follows that

We can rewrite the update of ηi (t) and si (t) by using an
aggregate dynamics. That is, let η(t) = [η1 (t) . . . ηN (t)]T
and s(t) = [s1 (t) . . . sN (t)]T be the aggregate states, their
dynamics is given by

T

φik (t) > µ.

k=1

N
X
k=1

Thus,

si (t)
ηi (t)

→

PN
k=1 σk (0)
PN
,
k=1 nk (0)

N
X

σk (0)

k=1

ηk (0) = φ̄i (t)

N
X

nk (0).

k=1

so that the proof follows.

Remark III.6 (Connectivity assumption). The result of Proposition III.5 can be proven also if Assumption II.1 is replaced
by Assumption 1 in [31], that is if {W (t)}t≥0 > 0 is a
stationary and ergodic sequence of stochastic matrices with
positive diagonals, and E[W ] is irreducible. The result can be
proven by following the same line of proof developed therein.
Thus, the ad-hoc estimator could be implemented also in a
monitoring network with stochastic gossip communication.

6

Remark III.7. The proposed algorithms are based on the assumption that the arrival-rates are constant in the observation
window. Therefore, one can apply the algorithms iteratively by
recomputing the estimates on different windows of data, thus
getting time-varying arrival-rates.
To conclude this section, we point out that the update of the
subgradient push optimization algorithm includes a push-sum
consensus step, i.e., a diffusive update based on a column
stochastic matrix (with coefficients wik (t) = 1/dk (t) for
k ∈ Ni (t)), which is the same used in the ad-hoc distributed
estimator (Algorithm 2). However, in the subgradient push
this update is part of a gradient descent step. In fact, the role
and the evolution of the involved variables, i.e., (si , ηi ) and
(vi , yi ) respectively, are different as well as the convergence
rates. Indeed, the ad-hoc estimator exhibits the exponential
convergence of linear√consensus protocols as opposed to the
much slower O(ln t/ t) rate of the subgradient-push [13].
These considerations, together with the lower computational
burden, make the ad-hoc distributed estimator appealing even
though not optimal. In the following section we will show that
it actually performs very closely to the MMSE estimator.
IV. E STIMATOR PERFORMANCE ANALYSIS
In this section we analyze the performance of the proposed
distributed estimators. In particular, for the Empirical Bayes
distributed estimator, lacking a closed form for the update rule
of b̂ML and, in turn, of λ̂EB , we are able to derive only steadystate (t → ∞) and asymptotic (N → ∞) bounds. Conversely,
for the ad-hoc distributed estimator a transient analysis (at any
t ∈ Z≥0 ) can be derived for any N ∈ N.
To develop the performance analysis, we consider a “special” agent, we label it as j, that does not participate to the
computation of b̂ML (respectively b̂hom ). Under this assumption,
it turns out that σj is independent of b̂ML (respectively b̂hom ).
Clearly, σj is also independent of any local estimate b̂ML
i (t)
hom
(respectively b̂i (t)), i ∈ {1, . . . , N }, at any t ∈ Z≥0 .
Consistently, here we can simply assume that for the comad-hoc
putation of λ̂EB
(t)), agent j uses the
j (t) (respectively λ̂j
local estimate of b computed by one of its in-neighbors, that
hom
I
is, e.g., b̂hom
j (t) = b̂i (t) for some i ∈ Nj (t) ⊂ {1, . . . , N }.
Notice that, in practice, the analysis developed for such
an agent j holds approximately for any node in the network
participating to the distributed computation. Indeed, due to
the large number of agents N , for any agent i (running
the distributed algorithm), b̂ML (respectively b̂hom ) and σi are
very weakly correlated. The validity of this statement will be
corroborated by the Monte Carlo analysis in the next section.
We start by deriving the Cramer-Rao lower bound (CRB)
for any unbiased estimator b̂ of b, i.e.,
VAR[b̂] ≥ CRB(b).
Lemma IV.1. The CRB for the estimation of the hyperparameter b is given by
b/a
.
(10)
CRB(b) = PN
n
i

i=1 ni b+1

The proof is reported in Appendix A.

A. Analysis of the Empirical Bayes distributed estimator
For the Empirical Bayes distributed estimator we can analyze the performance only once consensus on the optimal
value b̂ML , and thus on λ̂EB
j , has been reached. Specifically, let
us recall that
λ̂EB
j =

b̂ML
b̂ML nj + 1

(a + σj ).

(11)

Despite no analytical expression is available for the ML
estimator b̂ML (nor for its moments) the asymptotic analysis
(N → ∞) of the corresponding MMSE estimator λ̂EB
can
j
be obtained by using the properties of ML estimation. In
particular, any ML estimator is asymptotically unbiased and
efficient, hence E[b̂ML ] → b and VAR[b̂ML ] → CRB(b) as
N → ∞. From equation (10),
CRB(b) = PN

b/a

ni
i=1 ni b+1

≤ PN

b/a

1
i=1 nmax b+1

=

N →∞

b/a
N

1
nmax b+1

−−−−−→ 0 ,

(12)

so that VAR[b̂ML ] → 0 as N → ∞.
These results on the asymptotic properties of b̂ML can be
used to prove the following theorem characterizing the asymptotic behavior of the Empirical Bayes distributed estimator .
Due to the nonlinear dependency of λ̂EB from b̂ML , we will
perform an approximate analysis by considering the Taylor
expansion of λ̂EB around E[b̂ML ]. For tractability we will
consider respectively the second-order and the first-order approximations for mean and variance. The numerical analysis in
Section V will confirm the validity of such an approximation.
Theorem IV.2. Consider a network of monitors as in Section II running the Empirical Bayes distributed estimator.
Then, as N → ∞, it holds true
b
E[λ̂EB
(a + nj λj ).
(13)
j |λj ] −→
1 + nj b
to second order and
VAR[λ̂EB
j |λj ]


−→

b
1 + nj b

2
nj λ j ,

(14)

to first order.
B. Analysis of the ad-hoc distributed estimator
The closed-form update of the ad-hoc distributed estimator
allows us to perform a more detailed analysis. In particular,
we are able to characterize mean and variance of the local
estimator λ̂ad-hoc
(t) during the algorithm evolution (transient
j
analysis) and for any fixed value of the number of nodes N .
As we have done for the Empirical Bayes distributed
estimator, also for the ad-hoc one we first characterize the
estimator of the hyperparameter b.
Before addressing the transient analysis, we compute mean
and variance of the consensus value b̂hom .
Proposition IV.3. The homogeneous estimator of b, b̂hom , is
unbiased, i.e.,
E[b̂hom ] = b,

7

assumption that the measurements of agent j do not contribute
2
to the computation of b̂hom
j (t), we clearly have

and has variance
VAR[b̂hom ] =

N
2 X

b
b
+ 2
an an

n2i .

(15)

hom
E[b̂hom
j (t)|λj ] = E[b̂j (t)] = b.

i=1

Moreover, the estimator b̂hom is consistent, i.e., it converges in
probability to the true value b as N → ∞.

(18)

and
hom
VAR[b̂hom
j (t)|λj ] = VAR[b̂j (t)].

The proof is given in Appendix C.
Remark IV.4. The estimator b̂hom is the ML estimator in the
homogenous case, i.e., when ni = n/N for all i ∈ {1, . . . , N },
and it attains the Cramer-Rao Bound (CRB) not only asymptotically, but for any N . It follows by substituting ni = n/N
in (10) and (15). Interestingly, we will show in the numerical
analysis that even in the non-homogeneous scenario the ML
estimator b̂ML approaches the CRB at any fixed N .
The following lemma gives a characterization of the transient local estimates.
Lemma IV.5. Let Assumption II.1 and Assumption III.2 hold.
For all i ∈ {1, . . . , N }, the local estimator b̂hom
i (t) used in the
ad-hoc distributed estimator (Algorithm 2) is unbiased at any
t ∈ Z≥0 , i.e.,
E[b̂hom
i (t)] = b,
and has variance
VAR[b̂hom
i (t)]

Remark IV.7. It is worth noticing once more that although
for a tractable, rigorous analysis we need the assumption that
agent j does not contribute to the distributed computation, in
practice the results hold with good approximation also in the
scenario in which agent j contributes to the distributed computation. This is due to the weak impact of single measurements
onto the aggregate quantities. In fact, for example, following
analogous calculations as in the proof of Lemma IV.5, the
conditional mean in this latter case turns out to be


λj
φjj (t)nj
hom
b−
,
(19)
E[b̂j (t)|λj ] = b −
ηj (t)
a
where φjj (t) is, again, the element (j, j) of Φ(t) and we have
also used that, conditioned to the λj of agent j,

λj i = j
E[yi,` |λj ] =
.
ab i 6= j

φjj (t)nj
n
PN
PN
2
2 2 Since as t → ∞,
→ nj  1, then the difference
ηj (t)
b2
b
k=1 φik (t) nk
k=1 φik (t) nk
= P
2 + P
2 between equations (19) and (18) is practically negligible.
a
a
N
N
φ
(t)n
k
k=1 ik
k=1 φik (t)nk
In the next theorem we provide explicit transient expressions
(16) for the conditional mean and variance of λ̂ad-hoc (t). Again, we
j

where φik (t) is the element (i, k) of the state transition matrix
Φ(t) defined in (8).
Moreover, as t → ∞, the variance of b̂hom
i (t) converges
exponentially to the variance of b̂hom , (15), and satisfies
hom
|VAR[b̂hom
]| ≤
i (t)] − VAR[b̂

b(1 + 2bnmax )
δ(Φ(t)), (17)
µa

where µ > 0 is given in Lemma III.4 and δ(Φ(t)) is a (proper)
coefficient of ergodicity1 (exponentially decaying with time)
def
defined as δ(A) = max max |ai1 ,k − ai2 ,k |.
k

i1 ,i1

The proof is reported in Appendix D.
Remark IV.6. It is worth noting that µ and δ(Φ(t)) in
equation (17) typically depend on N (and thus also on n).
The interesting aspect of the given result is that the bound
provided in the previous lemma relates the convergence rate
of VAR[b̂hom ] to parameters of the communication graph
and the diffusion protocol as µ and δ(Φ(t)). For example,
(t+1)Q−1
if ∪τ =tQ
G(τ ) is balanced, the matrix Φ(t) is doubly
stochastic, hence µ = 1.
With this characterization of the hyperparameter estimator,
we are ready to analyze the estimator λ̂ad-hoc
(t). Recalling the
j
1A

coefficient of ergodicity is a function τ (·) continuous on the set of
raw (respectively column) stochastic matrices satisfying 0 ≤ τ (A) ≤ 1. It is
proper if τ (A) = 0 if and only if A = 1v T , with 1 = [1 . . . 1]T and v a
stochastic vector, [29], [32].

use second-order and first-order approximations respectively.
Theorem IV.8. Consider a network of monitors as in Section II running the ad-hoc distributed estimator (Algorithm 2).
Then the transient conditional mean and variance of λ̂ad-hoc
(t)
j
are given by
E[λ̂ad-hoc
(t)|λj ] ≈
j


b
nj
hom
(a + nj λj )
−
VAR[
b̂
(t)]
(20)
j
1 + nj b (1 + nj b)3
and
VAR[λ̂ad-hoc
(t)|λj ] ≈
j

2
b
nj
hom
−
VAR[
b̂
(t)]
nj λ j
j
1 + nj b (1 + nj b)3

2


2nj
2
+
VAR[b̂hom
j (t)] (a + nj λj ) + nj λj
3
(1 + nj b)
(21)
where ≈ indicates respectively the second-order and the
first-order approximations, and VAR[b̂hom
j (t)] is given in (16).
Moreover, the asymptotic value λ̂ad-hoc
satisfies, as N → ∞,
j
E[λ̂ad-hoc
|λj ] −→
j

b
(a + nj λj )
1 + nj b

(22)

2 As stated before, agent j uses as b̂hom (t) the estimate b̂hom (t) of some
j
i
neighbor, i ∈ NjI ⊂ {1, . . . , N } participating to the distributed computation.

8

to second order and
VAR[λ̂ad-hoc
|λj ]
j


−→

b
1 + nj b

2

0.3

nj λ j

(23)
0.104

Remark IV.9. By comparing equations (22)-(23) with (13)(14), one can notice that at steady-state (t → ∞) and for large
number of agents (N → ∞) the mean and variance of the adhoc distributed estimator approach the ones of the optimal
Empirical Bayes. Moreover,

 the right-hand side of (23) can
be rewritten as

λj
nj

nj b
1+nj b

2

R MSE [ b̂ hom
(t)]
i

0.25

to first order.

0.102
0.2

0.1

0.098
95

100

105

0.15

0.1

, which is clearly smaller than the

variance of the decentralized estimator VAR[λ̂dec
j |λj ] =

λj
nj .

V. N UMERICAL PERFORMANCE ANALYSIS
In this section we analyze the performance of the proposed
estimators. Starting from the theoretical characterization developed in the previous section, we perform a Monte Carlo
analysis confirming the theoretical bounds and adding other
insights on the performance of the estimators.
As performance metric we adopt the Root Mean Square
Error (RMSE), thus taking into account both bias and variance
of the estimators. We recall that for an estimator ω̂ of a
parameter ω, the RMSE is defined as
p
RMSE[ω̂] = E[(ω̂ − ω)2 ]
p
= VAR[ω̂] + E2 [ω̂] − 2ωE[ω̂] + ω 2 .
(24)
Clearly, if the estimator is unbiased the RMSE
p coincides with
the standard deviation, i.e., RMSE[ω̂] = VAR[ω̂].
The statistical RMSE, RMSE[ω], will be compared with
the sample value obtained through the Monte Carlo trials,
RMSEM C , computed as
v
u
M
u 1 X
(ω̂[m] − ω)2 ,
RMSEM C = t
M m=1
where M is the number of trials and ω̂[m] is the mth estimate
of ω.
In the following we set M = 5 × 104 and, to generate the
random values, we use a Gamma distribution with parameters
a = 10 and b = 1, which gives values of λi in the range
[1, 25] with 99.98% probability.
In order to challenge the ad-hoc distributed estimator we
focus on a strongly inhomogeneous network scenario. That
is, we consider a network in which half of the nodes have
the maximum number of measurements in the observation
window, ni = nmax (we set nmax = 50), and the remaining
ones only one measurement, ni = 1.
We start by analyzing the transient performance of the adhoc distributed estimator. Consistently with the theoretical
analysis, we first focus on the time evolution of the RMSE
hom
of b̂hom
i (t). Notice that, since
qb̂i (t) is an unbiased estimator
it holds RMSE[b̂hom
VAR[b̂hom
i (t)] =
i (t)]. We compare the
evolution of the sample RMSE with the theoretical expression
obtained from (16). For this analysis we set N = 20 and
consider two possible communication models. In the first one
we consider a fixed, directed and very sparse communication

20

40

60

80

100
t im e

120

140

160

180

Figure 3. Time evolution of RMSE[b̂hom
i (t)] under fixed communication
graph. The solid lines indicate the theoretical expressions, while the dash-dot
lines are the ones obtained via Monte Carlo simulations. The dashed horizontal
line is the theoretical consensus value.

graph defined as follows: starting from a directed cycle, we
add the edges (3, 1), (3, 2), (4, 1) and (4, 2) to unbalance the
graph. In Figure 3 we plot the evolution of RMSE[b̂hom
i (t)] for
4 nodes, two of them with nmax measurements and two with
one measurement. The theoretical curve predicts very accurately the sample RMSE obtained by the Monte Carlo trials,
as highlighted in the inset. As expected, the RMSE of the
different nodes converges to the consensus value RMSE[b̂hom ]
obtained from (15).
As a second scenario we challenge the algorithm on a timevarying topology. Namely, we consider a graph obtained by
extracting at each time-instant an Erdős-Rényi graph with
parameter 0.01. We choose a small value, so that at a given
instant the graph is disconnected with high probability, but it
turns out to be uniformly jointly connected with Q = 36.
In Figure 4 we again compare the theoretical evolutions
of RMSE[b̂hom
i (t)] with their sample counterparts. We can
highlight two main differences with respect to the previous
scenario. The curves have some constant portions showing
that nodes can be isolated for some time-intervals. However,
the convergence is faster compared to the fixed scenario. This
can be explained by the higher density of the union graph
in the time-varying scenario as opposed to the sparsity of the
fixed graph. In fact, we noticed that increasing the Erdős-Rényi
graph parameter increases the convergence speed.
Now we focus on the transient behavior of RMSE[λ̂ad-hoc
(t)]
i
for a generic node i. Notice that we use the index i, rather
than j, meaning that in the Monte Carlo trials the node under
investigation participates to the computation of b̂hom
i (t). This
allows us to show that the uncorrelation assumption in Section IV is in fact reasonable. Moreover, the computations will
also confirm the validity of the low-order approximation made
to derive the theoretical expressions of mean and variance.
In Figure 5 we compare the theoretical and sample evolution
of RMSE[λ̂ad-hoc
(t)]. The theoretical curve is obtained by
i
plugging equations (20) and (21) into equation (24). We
compute the curves for two different values of N , namely
N = 20 and N = 50. The difference between the theoretical
and sample curves is already minimal for N = 20 (showing a

9

0.45

0.4

R MSE [ b̂ hom
(t)]
i

0.35

0.3

0.25

0.2

0.15

0.1
10

20

30

40
t im e

50

60

70

80

Figure 4. Time evolution of RMSE[b̂hom
i (t)] under time-varying communication graph. The solid lines indicate the theoretical expressions, while the
dash-dot lines are the ones obtained via Monte Carlo simulations. The dashed
horizontal line is the theoretical consensus value.

2.3
2.2

N=20 Monte Carlo
N=20 theoretical
N=20 steady−state
N=50 Monte Carlo
N=50 theoretical
N=50 steady−state

1.72
1.7

2.1 1.68

RMSE

1.66
2 1.64
1.62
1.9

1.6
40

60

80

1.8
1.7
1.6
0

50

100

150

200
time

250

300

350

400

Figure 5. Time evolution of RMSE[λ̂ad-hoc
(t)] under fixed communication
i
graph for N = 20 (red lines) and N = 50 (black lines). The solid lines
indicate the theoretical expressions, while the dash-dot lines are the ones
obtained via Monte Carlo simulations. The dashed horizontal lines are the
theoretical steady-state values.

very weak correlation between b̂hom
i (t) and σi ) and completely
disappears for N = 50 (showing that the correlation has
practically no more influence). We want to stress that running
the same computation for a node not participating to the
distributed computation (hence matching the uncorrelation
assumption) the theoretical and sample curves are indistinguishable even for N = 20. This suggests that the low-order
approximation does not affect the goodness of the prediction.
It is worth noting that by increasing N the steady-state
value, RMSE[λ̂ad-hoc
], decreases, since the hyperparameter b
i
is estimated by means of a larger sample. This aspect will
be better highlighted in the following asymptotic analysis in
which we focus on how RMSE[λ̂ad-hoc
] varies with N .
i
In the asymptotic analysis we consider both the Empirical
Bayes distributed estimator and the ad-hoc distributed estimator by comparing again the predicted theoretical values with
the sample counterparts. As in the transient analysis we first
focus on the estimation of b.
In Figure 6 we plot the sample RMSE of the two estimators
(ML and homogeneous) and compare them with the theoretical
value of the homogeneous estimator and with the Cramer-Rao

Bound (CRB). As expected, for each fixed N the RMSE of
the ML estimator is lower than the homogeneous one due
to the (strong) inhomogeneity of the network. Once again,
we recall that the homogeneous estimator coincides with the
ML estimator only when the network is homogeneous (i.e.,
ni = n/N for any i ∈ {1, . . . , N }). As already experienced
in the transient analysis, the theoretical values of RMSE[b̂hom ]
practically coincide with the sample ones. Although for the
ML estimator we have no theoretical expression for fixed N ,
the picture shows a very interesting property. That is, the ML
estimator achieves the CRB not only asymptotically (N →
∞) as predicted by the theory, but also for each fixed N .
Interestingly, in accordance to the theoretical results in the
previous section, also the homogeneous estimator achieves the
CRB as N goes to infinity.
Finally, we analyze the RMSE of the estimators of λi . We
consider an agent with one measurement (ni = 1) and use the
most frequent value for the arrival rate, i.e., the mode of the
Gamma distribution, λi = (a − 1)b.
In Figure 7 we plot the sample RMSE of the Empirical
Bayes and ad-hoc estimators. We also plot the (sample)
values of the decentralized estimator. We decided to normalize
all the curves to the theoretical value of the decentralized
estimator in order to highlight the improvements of the proposed distributed estimators. Clearly, the sample values of the
decentralized estimator are approximately equal to one, with
minor fluctuations only due to the finite number of samples.
We compare the sample curves with the theoretical curve of
the homogeneous estimator and with the theoretical asymptotic
value as N → ∞. The theoretical curve is obtained as follows:
E[λ̂ad-hoc
] and VAR[λ̂ad-hoc
] can be computed by plugging
i
i
hom
VAR[b̂ ] from (15) in (20) and (21); then, RMSE[λ̂ad-hoc
] is
i
obtained by plugging E[λ̂ad-hoc
] and VAR[λ̂ad-hoc
] into (24).
i
i
Again, although computed under the uncorrelation assumption and neglecting higher order terms, the theoretical expression RMSE[λ̂ad-hoc
] predicts very accurately the sample values
i
(cross markers and solid curve).
The plot confirms how the distributed estimators take advantage from the network growth although the local sample
remains constant (even ni = 1). Indeed, the RMSE decreases
as N grows. The Empirical Bayes distributed estimator, being
the optimal estimator, always outperforms the ad-hoc distributed estimator. However, as predicted by the theoretical
analysis, the two estimators achieve the same asymptotic limit
as N → ∞. Moreover, it is interesting to notice that the RMSE
of the two estimators practically coincide already for N = 16,
thus strengthening the already appealing features of the ad-hoc
distributed estimator found from the theoretical analysis (i.e.,
easier computation and asymptotic optimality).
VI. C ONCLUSIONS
In this paper we have proposed a novel distributed scheme,
based on a hierarchical framework, for the Bayesian estimation of arrival rates in asynchronous monitoring networks.
The proposed distributed approach allows each node to gain
information from the network and thus outperforms the decentralized estimator, especially when the node local information

10

is obtained from the likelihood (3) as:

 2
∂ log L(y1 , . . . , yN |b)
,
I(b) = −E
∂b2
and, computing the derivatives,
"
#
N
aN X 2ni b + 1
I(b) = −E
−
(σi + a)
b2
b2 (ni b + 1)2
i=1

0.35
b̂
0.3

Mont e C ar lo

b̂ h o m Mont e C ar lo
b̂ h o m t he or e t ic al
CRB

0.25

R MSE

ML

0.2

0.15

N

0.1

=−
0.05

0

1

2

10

10
N

Figure 6. Sample RMSE for the homogeneous (triangles markers) and ML
(circles markers) estimators of b as a function of the number of agents N .
The sample values are compared with the Cramer-Rao Bound (dashed line)
and the theoretical RMSE of the homogeneous estimator (solid line).

Now, the mean of σi turns out to be E[σi ] = E[E[σi |λi ]] =
E[ni λi ] = ni ab, so that, after some manipulation

N 
a X
2ni b + 1
(ni b + 1)
I(b) = 2
−1 +
b i=1
(ni b + 1)2
=

1

so that the proof follows.
B
λ̂ E
Mont e C ar lo
i

0.9
N or m aliz e d R MSE

N
a X ni b
,
b2 i=1 ni b + 1

λ̂ id e c Mont e C ar lo

0.95

λ̂ ai d − h o c Mont e C ar lo
0.85

λ̂ ia d − h o c t he or e t ic al
as y mpt ot ic limit

0.8
0.75

B. Proof of Theorem IV.2
The MMSE estimator of λj in (11) can be written as λ̂EB
j =
def
b̂ML
ML
ML def
Z(b̂ )Y (σj ), where Z(b̂ ) = 1+b̂ML n and Y (σj ) = a+σj .
j

Due to the independence between b̂ML and σj , we have

0.7

E[λ̂EB
j |λj ] = E[ZY |λj ] = E[Z|λj ]E[Y |λj ]

0.65
0.6

(25)

and

0.55
0.5

aN X 2ni b + 1
+
(E [σi ] + a).
b2
b2 (ni b + 1)2
i=1

1

2

10

10

2
VAR[λ̂EB
j |λj ] = E [Z|λj ]VAR[Y |λj ]

N

Figure 7. Sample RMSE for the ad-hoc (triangle markers) and Empirical
Bayes (circle markers) estimators normalized with respect to the decentralized
estimator of λi (dot markers) as a function of the number of agents N . The
theoretical RMSE of the ad-hoc estimator (solid line) and the theoretical limit
(dashed line) are shown for comparison.

is scarce. In particular, the distributed estimator consists of
the convex combination of a global information, computed
through a distributed optimization algorithm, and a function
of the local data. Then we have proposed an ad-hoc distributed
estimator that performs closely to the optimal Empirical Bayes
estimator, but is much simpler to implement and exhibits
faster (exponential) convergence. We have analyzed the two
estimators and provided expressions for mean and variance
as the network size goes to infinity, showing that in this
asymptotic situation the ad-hoc estimator achieves the same
RMSE of the optimal one. Moreover, for the ad-hoc estimator
we have provided transient expressions for mean and variance.
A numerical Monte Carlo analysis has been performed to
corroborate the theoretical results and highlight the interesting
features of the two distributed estimators.
A PPENDIX
A. Proof of Lemma IV.1
1
,
The Cramer-Rao Bound is defined as CRB(b) = I(b)
where I(b) is the Fisher Information. The Fisher Information

+VAR[Z|λj ](E2 [Y |λj ] + VAR[Y |λj ]). (26)
The conditional moments of Y (σj ) are easily obtained
E[Y (σj )|λj ] = a + nj λj

(27)

VAR[Y (σj )|λj ] = nj λj .

(28)

and
ML

For the nonlinear function Z(b̂ ) we resort to an approximate
analysis. That is, we consider the Taylor expansion for the
ML
def
moments of the function Z(b̂ML ) = 1+b̂b̂ML n around the mean
j

value E[b̂ML |λj ] = E[b̂ML ].
E[Z(b̂ML )|λj ] = E[Z(b̂ML )]
h
≈ E Z(E[b̂ML ]) + Z 0 (E[b̂ML ])(b̂ML − E[b̂ML ])
i
1
(29)
+ Z 00 (E[b̂ML ])(b̂ML − E[b̂ML ])2 ,
2
where we have neglected terms of order higher than two in
the expansion of Z, hence the ≈ symbol. Then
1
E[Z(b̂ML )|λj ] ≈ Z(E[b̂ML ]) + Z 00 (E[b̂ML ])VAR[b̂ML ]. (30)
2
2n
Observing that Z 0 (s) = (1+n1j s)2 and Z 00 (s) = − (1+njjs)3 ,
and plugging (27) and (30) in (25), we get
E[λ̂EB
j |λj ] ≈
"
(a + nj λj )

E[b̂ML ]
E[b̂ML ]nj + 1

−

nj
(E[b̂ML ]nj + 1)3

#
ML

VAR[b̂

] .

11

From the asymptotic properties of the ML estimators, we know
that for N → ∞, E[b̂ML ] → b and VAR[b̂ML ] → 0, so that
equation (13) follows.
Similarly, it turns out that

2
VAR[Z(b̂ML )|λj ] ≈ Z 00 (E[b̂ML ]) VAR[b̂ML ].
(31)
Using (30)-(31) into eq. (26) we obtain
!2
E[b̂ML ]
EB
nj λj
VAR[λ̂j |λj ] ≈
E[b̂ML ]nj + 1
"
#2


2nj
VAR[b̂ML ] (a + nj λj )2 + nj λj .
+
ML
3
(E[b̂ ]nj + 1)
Using again the asymptotic properties of b̂ML , equation (14)
follows, thus concluding the proof.
C. Proof of Proposition IV.3
First, notice that E[σi |λi ] = ni λi , thus
"N
#
"N
#
X
X
E[σ] = E[E[σ|λ]] = E
E[σi |λi ] = E
ni λi
i=1

=

N
X

i=1

ni ab = nab,

i=1

so that E[b̂hom ] = b.
To prove the second part,
i ] = ni λi
PNrecall that VAR[σi |λ
P
N
and, thus, VAR[σ|λ] =
VAR[σ
|λ
]
=
i
i
i=1
i=1 ni λi .
Using the law of total variance, we have
VAR[σ] = E[VAR[σ|λ]] + VAR[E[σ|λ]]
"N
#
"N
#
X
X
=E
VAR[σi |λi ] + VAR
E[σi |λi ]
=E

i=1
"N
X

i=1

#

"

ni λi + VAR

i=1

N
X

#
ni λi

D. Proof of Lemma IV.5
To prove that b̂hom
i (t) is unbiased at any t ∈ Z≥0 , first let us
recall that the aggregate states η(t) = [η1 (t) . . . ηN (t)]T and
s(t) = [s1 (t) . . . sN (t)]T evolve according to the dynamics (7)
def
def
with s(0) = [σ1 . . . σN ]T = σ, η(0) = [n1 . . . nN ]T = n and
W (t) the column stochastic matrix with elements wij (t). The
evolution of s(t) and η(t) is given by (9) where Φ(t) is the
deterministic state transition matrix defined in (8).
Next, observe that the update ηi (t) depends only on the
initial value ni and thus is deterministic. Therefore, using the
update in Algorithm 2, we have
1
E[si (t)]
E[b̂hom
i (t)] =
aηi (t)
Using the evolution of the aggregate state and denoting ei the
ith canonical vector (e.g., e1 = [1 0 . . . 0]T ), we have
1
E[eTi (Φ(t)σ)]
E[b̂hom
i (t)] =
aηi (t)
1
=
eT Φ(t)E[σ].
aηi (t) i
Noting that E[σi ] = abni it follows
1
eT Φ(t)ab n
E[b̂hom
i (t)] =
aηi (t) i
b T
=
e η(t) = b
ηi (t) i
where the last two steps follow respectively from (9) (Φ(t)n =
η(t)) and from eTi η(t) = ηi (t).
Next, we show that the transient variance is given by
equation (16).
Using again the update in Algorithm 2, it holds
1
VAR[b̂hom
VAR[eTi (Φ(t)σ)].
i (t)] = 2
a ηi (t)2
Noting that eTi Φ(t) is the ith row of Φ(t), it follows
"N
#
X
1
hom
φik (t)σk ,
VAR
VAR[b̂i (t)] = 2
a ηi (t)2
k=1

i=1

= abn + ab2

N
X

n2i .

i=1

The variance of the
estimator is given by

 σ homogeneous
, so that equation (15) follows.
VAR[b̂hom ] = VAR an
By Markov’s inequality and Lemma IV.3
n σ
o 1
h
i
Pr
− b ≥  ≤ VAR b̂hom ,
an

and the variance can be bounded as follows
N
h
i
b
b2 X 2
VAR b̂hom =
+ 2
n
an an i=1 i

≤

N
b
b2 X 2
+
n
aN
aN 2 i=1 max

=

b
nmax 2 b2 N →∞
+
−−−−−→ 0
aN
aN

so that b̂hom → b in probability, thus concluding the proof.

where φik (t) is the element (i, k) of the matrix Φ(t). By the
law of total variance we have that
( "
"N
##
X
1
hom
VAR[b̂i (t)] = 2
E VAR
φik (t)σk λ
a ηi (t)2
k=1
" "N
##)
X
+VAR E
φik (t)σk λ
1
= 2
a ηi (t)2

k=1
N
X

( "
E

φik (t) nk λk

k=1

"
+VAR

#
2

N
X

#)
φik (t)nk λk

,

k=1

where we have used E[σk ] = nk λk and VAR[σk ] = nk λk .
Finally, recalling that E[λk ] = ab, VAR[λk ] = ab2 and the
λk s are independent, it turns out
!
N
N
X
X
b
hom
2
2 2
φik (t) nk + b
φik (t) nk ,
VAR[b̂i (t)] =
aηi (t)2
k=1

k=1

12

so that equation (16) follows.
To prove the asymptotic result, we work
hom
|VAR[b̂hom
]| by using (16) and (15)
i (t)] − VAR[b̂

from the definition of the coefficient of ergodicity, it follows
out

hom
|VAR[b̂hom
]| =
i (t)] − VAR[b̂


PN
2
φik (t) nk
1
b
=
 Pk=1
2 − 
a
n
N
k=1 φik (t)nk


PN
PN
2 2
2
2
φik (t) nk
b 
k=1 nk 
+  Pk=1

2 −
a
n2
N
φ
(t)n
ik
k
k=1

!2 
N
N
X
X
1
b 
=
n
φik (t)2 nk −
φik (t)nk 
ηi (t)2 an
k=1
k=1


!2 N
N
N
X
X
b2  2 X
+ 2 n
φik (t)2 n2k −
φik (t)nk
n2k 
an
k=1

k=1

k=1

PN
By using the definition of n =
n and writing
P
2 P
h=1
P h

N
N
N
=
k=1 φik (t)nk
k=1 φik (t)nk
h=1 φih (t)nh
hom
|VAR[b̂hom
]| =
i (t)] − VAR[b̂

=

b
1
2
ηi (t) an
−

N
X

N
X

nh

h=1

φik (t)nk

k=1

b2
+ 2
an
−

N
X

N
X

N
X

φih (t)nh

h=1

`=1
N
X
`=1

n`

N
X

+

b
an2

k=1

φi` (t)n`

N
X

!
n2k

k=1

nh n` n2k φik (t)2 − φih (t)φi` (t)

h=1 `=1 k=1

"

N N
1
b XX
nh nk φik (t) |φik (t) − φih (t)|
ηi (t)2 an
2

+

b
an2

h=1 k=1
N
N X
N
XX


nh n` n2k φik (t) |φik (t) − φih (t)|

h=1 `=1 k=1

+ |φik − φi` (t)|



The local update of the ad-hoc distributed estimator (Algorithm 2) can be written as λ̂ad-hoc
(t) = X(b̂hom
j
j (t))Y (σj ),
def

where X(b̂hom
j (t)) =

b̂hom
j (t)
1+b̂hom
j (t)nj

def

and Y (σj ) = a + σj . Using

again the independence between b̂hom
j (t) and σj , we can write
ad-hoc
E[λ̂ad-hoc
(t)|λ
]
and
VAR[
λ̂
(t)|λ
j
j ] as
j
j
(32)

+ VAR[X|λj ](E2 [Y |λj ]+VAR[Y |λj ]). (33)

+ φik (t)φi` (t) − φik (t)φi` (t))
≤

E. Proof of Theorem IV.8

φik (t)2 n2k

h=1 k=1

N X
N X
N
X

thus concluding the proof.

VAR[λ̂ad-hoc
(t)|λj ] = E2 [X|λj ]VAR[Y |λj ]
j

N N
1
b XX
=
nh nk φik (t) (φik (t) − φih (t))
ηi (t)2 an
2

δ(Φ(t))b(1 + 2bnmax )
µa

and

φih (t)nh

nh

h=1

hom
|VAR[b̂hom
]| ≤
i (t)] − VAR[b̂

E[λ̂ad-hoc
(t)|λj ] = E[XY |λj ] = E[X|λj ]E[Y |λj ]
j

!

h=1
N
X

where we have simplified the commonPfactor ηi (t) =
P
N
N
k=1 φik (t)nk . Finally, by Lemma III.4,
k=1 φik (t) ≥ µ,
with µ > 0 so that

φik (t)2 nk

k=1
N
X

hom
|VAR[b̂hom
]|
i (t)] − VAR[b̂


1
b
b2 2
≤ PN
nδ(Φ(t)) + 2 n nmax 2δ(Φ(t))
an
an
k=1 φik (t)nk

#

where the last inequality follows by using the triangular
inequality. Now, weak ergodicity of W (t) implies that for any
i, k, h ∈ {1, . . . , N } it holds |φik (t) − φih (t)| ≤ M λt for
some M > 0 and λ > 0, see, e.g., [13, Corollary 8], so that
the exponential convergence of VAR[b̂hom
i (t)] follows. Then,

Considering, as in the previous theorem, the Taylor expansion
b̂hom
def
j (t)
for the moments of the function X(b̂hom
j (t)) = 1+b̂hom (t)n
j

j

around the mean E[b̂hom
j (t)|λj ] = b, we obtain
hom
E[X(b̂hom
j (t))|λj ] = E[X(b̂j (t))]
1 00
hom
hom
≈ X(E[b̂hom
j (t)]) + X (E[b̂j (t)])VAR[b̂j (t)]
2
b
nj
=
−
VAR[b̂hom
(34)
j (t)]
1 + nj b (1 + nj b)3

where VAR[b̂hom
j (t)] is the one given in (16) and, again, the
≈ symbol indicates that we have neglected higher-order terms
in the expansion of X. Hence, plugging (27) and (34) in (32),
equation (20) follows. Similarly,

2
0
hom
VAR[X(b̂hom
VAR[b̂hom
j (t))|λj ] ≈ X (E[b̂j (t)])
j (t)]
!

2
N
1
b
b2 X 2
=
+
n .
(35)
(1 + nj b)2
an an2 i=1 i
Using (34)-(35) into (33), equation (21) follows, thus concluding the first part of the proof.
To prove the asymptotic results, we just need to recall from
Proposition IV.3 that as N → ∞ VAR[b̂hom ] → 0, so that
equations (22) and (23) follow, thus concluding the proof.

13

R EFERENCES
[1] A. Coluccia and G. Notarstefano, “Distributed bayesian estimation of arrival rates in asynchronous monitoring networks,” in IEEE International
Conference on Acoustics, Speech and Signal Processing (ICASSP), 2014,
pp. 5050–5054.
[2] D. L. Snyder and M. I. Miller, Random Point Processes in Time and
Space. Springer, 1991.
[3] D. Heyman and M. Sobel, Stochastic Models in Operations Research:
Stochastic optimization, ser. Dover Books on Computer Science Series.
Dover Publications, 2003.
[4] S. Barbarossa, S. Sardellitti, and P. Di Lorenzo, Distributed detection
and estimation in wireless sensor networks, ser. Communications and
Radar Signal Processing. Academic Press Library in Signal Processing,
October 2013, vol. 2, ch. 7, pp. 329–408.
[5] I. D. Schizas, A. Ribeiro, and G. B. Giannakis, “Consensus in ad hoc
WSNs with noisy links Part I: Distributed estimation of deterministic
signals,” IEEE Transactions on Signal Processing, vol. 56, no. 1, pp.
350–364, 2008.
[6] F. S. Cattivelli and A. H. Sayed, “Diffusion lms strategies for distributed
estimation,” IEEE Transactions on Signal Processing, vol. 58, no. 3, pp.
1035–1048, 2010.
[7] S. Barbarossa and G. Scutari, “Bio-inspired sensor network design,”
IEEE Signal Processing Magazine, vol. 24, no. 3, pp. 26–35, 2007.
[8] F. Garin and L. Schenato, “A survey on distributed estimation and control
applications using linear consensus algorithms,” Networked Control
Systems, pp. 75–107, 2011.
[9] P. Frasca, H. Ishii, C. Ravazzi, and R. Tempo, “Distributed randomized
algorithms for opinion formation, centrality computation and power
systems estimation: A tutorial overview,” European Journal of Control,
2015.
[10] S. Sardellitti, M. Giona, and S. Barbarossa, “Fast distributed average
consensus algorithms based on advection-diffusion processes,” IEEE
Transactions on Signal Processing, vol. 58, no. 2, 2010.
[11] D. E. Marelli and M. Fu, “Distributed weighted least-squares estimation
with fast convergence for large-scale systems,” Automatica, vol. 51, pp.
27–39, 2015.
[12] F. Pasqualetti, R. Carli, and F. Bullo, “Distributed estimation via iterative
projections with application to power network monitoring,” Automatica,
vol. 48, no. 5, pp. 747–758, 2012.
[13] A. Nedic and A. Olshevsky, “Distributed optimization over time-varying
directed graphs,” arXiv preprint arXiv:1303.2289, 2013.
[14] S. Barbarossa and G. Scutari, “Decentralized maximum-likelihood estimation for sensor networks composed of nonlinearly coupled dynamical
systems,” IEEE Transactions on Signal Processing, vol. 55, no. 7, 2007.
[15] A. Chiuso, F. Fagnani, L. Schenato, and S. Zampieri, “Gossip algorithms
for simultaneous distributed estimation and classification in sensor
networks,” IEEE Journal of Selected Topics in Signal Processing, vol. 5,
no. 4, pp. 691–706, 2011.
[16] F. Fagnani, S. M. Fosson, and C. Ravazzi, “Input driven consensus algorithm for distributed estimation and classification in sensor networks,” in
50th IEEE Conference on Decision and Control and European Control
Conference (CDC-ECC). IEEE, 2011, pp. 6654–6659.
[17] D. Varagnolo, G. Pillonetto, and L. Schenato, “Distributed consensusbased bayesian estimation: sufficient conditions for performance characterization,” in American Control Conference. IEEE, 2010, pp. 3986–
3991.
[18] P. Di Lorenzo and S. Barbarossa, “Distributed least mean squares
strategies for sparsity-aware estimation over gaussian markov random
fields,” in IEEE International Conference on Acoustics, Speech and
Signal Processing (ICASSP), 2014, pp. 5472–5476.
[19] F. Cattivelli, C. G. Lopes, and A. H. Sayed, “Diffusion recursive
least-squares for distributed estimation over adaptive networks,” IEEE
Transactions on Signal Processing, vol. 56, no. 5, 2008.
[20] G. Mateos, I. D. Schizas, and G. B. Giannakis, “Distributed recursive
least-squares for consensus-based in-network adaptive estimation,” IEEE
Transactions on Signal Processing, vol. 57, no. 11, pp. 4583–4588, 2009.
[21] I. D. Schizas, G. B. Giannakis, S. I. Roumeliotis, and A. Ribeiro,
“Consensus in ad hoc WSNs with noisy links Part II: Distributed
estimation and smoothing of random signals,” IEEE Transactions on
Signal Processing, vol. 56, no. 4, pp. 1650–1666, 2008.
[22] D. A. Freedman, “Poisson processes with random arrival rate,” The
Annals of Mathematical Statistics, pp. 924–929, 1962.
[23] W. A. Massey, G. A. Parker, and W. Whitt, “Estimating the parameters of
a nonhomogeneous poisson process with linear rate,” Telecommunication
Systems, vol. 5, no. 2, pp. 361–388, 1996.

[24] D. Karlis and E. Xekalaki, “Mixed poisson distributions,” International
Statistical Review, vol. 73, pp. 35–58, 2005.
[25] C. Withers and S. Nadarajah, “On the compound poisson-gamma
distribution,” Kybernetika, vol. 47, no. 1, pp. 15–37, 2011.
[26] E. Lehmann and G. Casella, Theory of Point Estimation. Springer,
1998.
[27] F. Zanella, D. Varagnolo, A. Cenedese, G. Pillonetto, and L. Schenato,
“Asynchronous newton-raphson consensus for distributed convex optimization,” in 3rd IFAC Workshop on Distributed Estimation and Control
in Networked Systems (NecSys’12), 2012.
[28] J. N. Tsitsiklis, D. P. Bertsekas, M. Athans et al., “Distributed asynchronous deterministic and stochastic gradient optimization algorithms,”
IEEE Transactions on Automatic Control, vol. 31, no. 9, pp. 803–812,
1986.
[29] E. Seneta, Non-negative matrices and Markov chains. Springer, 2006.
[30] A. Jadbabaie, J. Lin, and A. S. Morse, “Coordination of groups of mobile
autonomous agents using nearest neighbor rules,” IEEE Transactions on
Automatic Control, vol. 48, no. 6, pp. 988–1001, 2003.
[31] F. Bénézit, V. Blondel, P. Thiran, J. Tsitsiklis, and M. Vetterli, “Weighted
gossip: Distributed averaging using non-doubly stochastic matrices,”
in IEEE International Symposium on Information Theory Proceedings
(ISIT). IEEE, 2010, pp. 1753–1757.
[32] N. H. Vaidya, C. N. Hadjicostis, and A. D. Dominguez-Garcia, “Distributed algorithms for consensus and coordination in the presence of
packet-dropping communication links-part ii: Coefficients of ergodicity
analysis approach,” arXiv preprint arXiv:1109.6392, 2011.

Angelo Coluccia (M’13) received the Eng. degree in
Telecommunication Engineering (summa cum laude)
in 2007 and the PhD degree in Information Engineering in 2011, both from the University of Salento,
Lecce, Italy. Former researcher at Forschungszentrum Telekommunikation Wien, Vienna, since 2008
he has been engaged in research projects on traffic
analysis, security and anomaly detection in operational cellular networks. He is currently Assistant Professor at the Dipartimento di Ingegneria
dell’Innovazione, University of Salento, where he
teaches the course of Telecommunication Systems. His research interests
are signal processing, communications and wireless networks, in particular
cooperative sensing/estimation approaches for localization and other (possibly
distributed) applications.

Giuseppe Notarstefano has been an Assistant Professor (Ricercatore) at the Università del Salento
(Lecce, Italy) since February 2007. He received the
Laurea degree “summa cum laude” in Electronics
Engineering from the Università di Pisa in 2003
and the Ph.D. degree in Automation and Operation
Research from the Università di Padova in April
2007. He has been visiting scholar at the Universities
of Stuttgart, California Santa Barbara and Colorado
Boulder. His research interests include distributed
optimization, cooperative control in multi-agent networks, applied nonlinear optimal control, and trajectory optimization and
maneuvering of aerial and car vehicles. He serves as an Associate Editor
in the Conference Editorial Board of the IEEE Control Systems Society
and for the European Control Conference, IFAC World Congress and IEEE
Multi-Conference on Systems and Control. He coordinated the VI-RTUS team
winning the International Student Competition Virtual Formula 2012. He is
recipient of an ERC Starting Grant 2014.

