Table Of ContentMaximum pseudo-likelihood estimator for
nearest-neighbours Gibbs point processes
6
0
0
By Jean-Michel Billiot1, Jean-Franc¸ois Coeurjolly1,2
2
n and R´emy Drouilhet1
a
J Labsad, University of Grenoble 2, France
4
January 3, 2006
]
T
S
.
h
t
a
m
[
1
v
5
6
0
1
0
6
0
/
h
t
a
m
:
v
i
X
r
a
1Supportedby agrant from IMAG Project AlpB.
AMS 2000 subject classifications. Primary 60K35,60G55; secondary 62M30.
Key words and phrases : Gibbs point processes, nearest-neighbours models, Delaunay graph, pseudo-likelihood
method,minimum contrast estimators.
2Corresponding author
1
2 J.-M. Billiot, J.-F. Coeurjolly and R. Drouilhet
Abstract
This paper is devoted to the estimation of a vector θ parametrizing an energy
function associatedto some“Nearest-Neighbours”Gibbs pointprocess,via the
pseudo-likelihoodmethod.Wepresentsomeconvergenceresultsconcerningthis
estimator, that is strong consistency and asymptotic normality, when only a
singlerealizationisobserved.Sufficientconditionsareexpressedintermsofthe
local energy function and are verified on some examples.
1 Introduction
Gibbs point processes first appeared in the theory of statistical physics. Historical aspects
of the mathematical theory are covered briefly in Kallenberg (1983). The importance of the
Gibbs point process as a model building principle became widely recognized through these
works. Indeed, the class of Gibbs point processes is interesting because it allows to introduce
andstudyinteractions betweenpointsthroughthemodellingofanassociated potentialfunction.
This resultinggain explains their usein statistical physicsRuelle (1969), Feynman (1972)(when
takinginteractionsbetweenmolecules inmodelsofdilutegases intoaccount)orinecology (when
analysingcompetitionsbetweenplants).Withinthemechanicsstatisticsframework,Gibbsstates
are defined as solutions of the well known equilibrium equations refered to Dobrushin-Lanford-
Ruelle (D.L.R.) equations Dobrushin (1969), Lanford and Ruelle (1969). One way to introduce
Gibbs point processes consists in using a family of local specifications with respect to a weight
process. The Preston’s theorems (Preston (1976)) used precisely this approach in order to give
sufficient conditions on local specifications for the existence of Gibbs states.
Manyproposalstriedtoestimatethepotentialfunctionfromtheavailablepointpatterndata
generated by some Gibbs point processes. If the potential belongs to a parametric family model,
the most well-known methodology is the use of the likelihood function. The main drawback of
this approach is that the likelihood function contains an unknown scaling factor whose value
depends on the parameters and which is difficult to calculate. The first class of models on which
the estimation of the maximum likelihood has been undertaken is the class of pairwise interac-
tion pointprocesses.Ogata and Tanemura(1984)developedthemaximumlikelihood estimation
method based on numerical approximations of the likelihood. Penttinen (1984) used a similar
approach while applying a Monte Carlo method in a way to solve the likelihood equation by the
stochastic Newton-Raphson algorithm. Moyeed and Baddeley (1991) proposed another iterative
procedureforestimatingthemaximumlikelihoodestimator.Formaximumlikelihood byMarkov
chain Monte Carlo, see Geyer and Thompson (1992), Geyer (1999) and for U.L.A.N. conditions
formaximumlikelihood estimator,seeMase(1992).Analternativeapproachconsistsinavoiding
to optimize the likelihood function (because of the scaling factor problem) and introducing a
Pseudo-likelihood for some Gibbs point processes 3
pseudo-likelihood function instead. This idea originated from Besag (1974) in the study of lat-
tice processes. Besag et al. (1982) further considered this method for pairwise interaction point
process, while Jensen and Møller (1991) generalized it to the general class of Gibbs point pro-
cesses, see Mase (1995), Mase (1999), Jensen and Ku¨nsch (1994), Guyon (1991) for asymptotic
properties. A third way is the Takacs-Fiksel estimation method (Takacs (1986) Fiksel (1988)),
which relies on a characteristic property of Gibbs processes using Palm measure. Asymptotics
propertiesofTakacs-FikselestimatorarestudiedinHeinrich(1992),Billiot(1997).Acomparison
of these different procedures applied to the Strauss model is presented in Diggle et al. (1994).
The non parametric setting has been undertaken by Gl¨otzl and Rauschenschwandtner (1981)
and Diggle et al. (1987) (and the references therein). Heikkinen and Penttinen (1999) proposed
a semiparametric estimator based on Bayesian smoothing techniques. A general review of the
problem of statistical inference on spatial point processes can be found in the recent monograph
of Møller and Waggepetersen (2003).
The present study is devoted to “Nearest-Neighbour” Gibbs point models by combining
stochasticgeometryarguments(Stoyan et al.(1995))andcomputationalgeometryones(Preparata and Shamos
(1988),Edelsbrunner(1988),Boissonnat and Yvinec(1995)).SuchmodelsareintroducedbyBaddeley and Møller
(1989) where the neighbourhood relation depends on the realization of the process. Sufficient
conditions (expressed in terms of the energy function) for the existence of such processes are
proposed in Bertin et al. (1999b) and Bertin et al. (1999a), where some examples are also pro-
posed. The main one is a pairwise interaction point process where the neighbourhood relation
corresponds to the (slightly modified) Delaunay graph of the realization of the process.
In this paper, we study a pseudo-likelihood estimator for such processes. More precisely, our
framework is restricted to stationary Gibbs point processes based on energy function related
to some graph (for instance the Delaunay graph) such that the energy function is invariant by
translation and such that the local energy function is stable and quasi-local (or local). The main
results of this paper are convergence results (strong consistency and asymptotic normality) of
maximum pseudo-likelihood estimators in this framework. These results are obtained when only
a single realization is observed. Sufficient conditions are expressed in terms of the local energy
function (which makes the results quite general) for some large family of parametrized energy
functions. Among the different parametrizations, the exponential family is considered.
The paper is organized as follows. Section 2 is devoted to some background on Gibbs point
processes and to the description of our framework. The statistical model and the pseudo-
likelihood method are presented in Section 3. Consistency and asymptotic normality of the
maximum pseudo-likelihood estimator are respectively proved in Section 4 and Section 5. Fi-
nally, the different sufficient conditions ensuring convergence results are verified on some ex-
amples in Section 6. A short simulation is presented to check the effectiveness of maximum
pseudo-likelihood estimator.
4 J.-M. Billiot, J.-F. Coeurjolly and R. Drouilhet
2 Background on Gibbs point processes
2.1 Gibbs point processes
We define B, B to be respectively the Borel σ-field and the bounded Borel boolean ring.
b
Let Ω denotes the class of locally finite subsets of IRd. In particular, an element ϕ of Ω,
also called configuration (of points), could be represented as ϕ = δ which is a simple
i∈IN xi
counting Radon measure in IRd (i.e. all the points x of IRd are distinct) where for every Λ ∈
i P
B, δ (Λ) = 1 (x) is the Dirac measure and 1 (.) is the indicator function of a set A. This space
x Λ A
Ω is equipped with the vague topology, that is to say the weak topology for Radon measures
with respect to the set of continuous functions vanishing outside a compact set. We also define
the σ-field F spanned by the maps ϕ −→ ϕ(Λ),Λ ∈B , where ϕ(Λ) corresponds to the number
b
of points of ϕ in Λ due to the Radon measure representation of ϕ. The set of all configurations
in a measurable set Λ⊂ IRd will be denoted by Ω and the corresponding σ-field F is similarly
Λ Λ
defined. Furthermore, for any Λ∈ B ,
b
(Ω,F) = (Ω ,F )×(Ω ,F )
Λ Λ Λc Λc
where Λc = IRd\Λ denotes the complementary of Λ in IRd. Finally, Ω denotes the class of all
f
finite subsets of IRd.
A point process on IRd is a Ω-valued random variable, denoted by Φ, with probability distri-
bution P on (Ω,F). and the intensity measure Λ of P is defined as a measure on B such that
p
for any D ∈ B
Λ (D) = ϕ(D) P(dϕ).
p
ZΩ
In the stationary case, Λ (D) = λ ν(D) where the constant λ is called the intensity of P and
p p p
ν is the Lebesgue measure on IRd.
A Gibbs point process is usually defined using a family of local specifications with respect to
a weight process (often a stationary Poisson process with distribution Q and intensity λ = 1).
Q
Let Λ be a bounded region in IRd. For such a process, given some configuration ϕ on Λc, the
Λc
conditional probability on Λ is of the form, for any Y ∈ F :
1
Π (ϕ,Y) = exp(−V (ψ|ϕ ))1 (ψ∪ϕ )Q (dψ) 1 (ϕ),
Λ Z (ϕ) Λc Y Λc Λ RΛ
(cid:26) Λ ZΩΛ (cid:27)
where
Z (ϕ)= exp(−V (ψ|ϕ ))Q (dψ)
Λ Λc Λ
ZΩΛ
is called the partition function and R = {ϕ ∈ Ω : 0 < Z (ϕ)< ∞}.
Λ Λ
Whereas the finite energy function V (ϕ) measures the cost of any configuration, the local
energy V (ψ|ϕ) is defined as the energy required to add the points of ψ in ϕ :
V (ψ|ϕ)= V (ψ∪ϕ)−V (ϕ).
Pseudo-likelihood for some Gibbs point processes 5
Let us notice that when ψ reduces to one point x, we denote by a slight abuse V (x|ϕ) instead
of V ({x}|ϕ). It is well known that the collection of probability kernels (Π ) satisfies the set
Λ Λ∈Bb
of compatibility and measurability conditions which define a local specification in the Preston’s
sense (Preston (1976)). The main condition is the consistency :
′
ΠΛΠΛ′ = ΠΛ for Λ ⊂ Λ.
Notice that some conditions are needed to ensure the existence of a probability measure P
withrespect toany local energyV andanyweightprocess that satisfies theso-called Dobrushin-
Lanford-Ruelle (D.L.R.) equations :
P(Y|F )(ϕ) =Π (ϕ,Y) for P a.e. ϕ ∈ Ω for any Λ∈ B and Y ∈ F.
Λc Λ b
For the general theory of Gibbs point processes, the reader may refer to Kallenberg (1983);
Daley and Vere-Jones (1988); Stoyan et al. (1995) and the references therein.
2.2 Campbell and Palm measures and Gl¨otz Theorem
The reduced Campbell measure C! of P is a measure on B⊗F such that for any D ∈ B and
p
any Y ∈ F
C!(D×Y) = 1 (ϕ−δ )ϕ(dx) P(dϕ).
p Y x
ZΩZD
When some measurable function h from IRd×Ω on IRd is given, the following equation is often
called the refined Campbell theorem
h(x,ϕ−δ )P(dϕ) = h(x,ϕ)C!(d(x,ϕ)).
x p
ZΩx∈ϕ ZIRd×Ω
X
If the intensity measure Λ is σ−finite, then for Λ - a.a. x ∈ IRd, the distribution P! on (Ω,F)
p p x
exists. It is unique for Λ - a.a. x ∈IRd and such that
p
C!(D×Y) = P!(Y) Λ (dx) for any D ∈ B, Y ∈ F.
p x p
ZD
Then P! is called the reduced Palm distribution of the point process P with respect to point x.
x
Intuitively, the Palm distribution P is the conditional probability of configurations of the point
x
process given that the point x belongs to the realization ϕ. Therefore, we have
h(x,ϕ−δ )P(dϕ) = h(x,ϕ)P!(dϕ)Λ (dx).
x x p
ZΩx∈ϕ ZIRd×Ω
X
When the process is stationary, one may apply the previous equation by replacing the intensity
measure Λ by its expression in this case λ ν, and in this framework, Gl¨otzl (1980) proved that
p p
6 J.-M. Billiot, J.-F. Coeurjolly and R. Drouilhet
P ∈ G (V) if and only if thereduced CampbellmeasureC! is absolutely continuouswith respect
0 p
to ν ×P and :
dC! dP!
p (x,ϕ) = λ x (ϕ)= exp(−V (x|ϕ))
P
d(ν ×P) dP
where λ = exp(−V (x|ϕ))P (dϕ) is the intensity of the process P. In the particular case
P
Ω
when V (x|ϕ)R= 0, the point process corresponds to the stationary Poisson process Q. We know
from the Slivnyak’s theorem that Q! = Q which is one way of characterizing such process.
x
2.3 Description of some Gibbs models
This paper is mainly devoted to the statistical studyof some nearest-neighbours Gibbs point
processes first introduced in Baddeley and Møller (1989). More precisely, we are interested in
models based on energy function of the form
3
V (ϕ) = u(k)(ξ,ϕ), (1)
Xk=1ξ∈DXelk(ϕ)
where Del (ϕ) is the set of clique of order k of the Delaunay graph defined just below. For some
k
ϕ ∈ Ω in general position, one defines Del (ϕ) by the unique decomposition into triangles ψ
3
in which the convex hull of the circle C(ψ) does not contain any point of ϕ\ψ. The Delaunay
graph is then defined by the set of edges :
Del (ϕ) = ∪ P (ψ).
2 ψ∈Del3(ϕ) 2
In order to ensure the existence of such Gibbs state in IRd, Bertin et al. (1999b) prove that the
local stability and quasilocality properties (only expressed in terms of the energy function) are
sufficient conditions of Preston’s Theorem.
Without any additional modification, the previous model does not satisfy the previous as-
sumptions.We then introducesome subgraphs.First let us denote, for some triangle ψ,byD(ψ)
the diameter of the circle circumscribed of ψ and by β(ψ) the smallest angle of ψ.
Definition 1 Given any β ∈]0,π/3], we introduce the following particular subset of Del (ϕ) :
0 3
β
Del 0 (ϕ) = {ψ ∈ Del (ϕ) :β(ψ)> β }.
3,β 3 0
The β-Delaunay graph of order β of any configuration ϕ is the Delaunay subgraph defined by :
0
β
Del 0 (ϕ)= P (ψ).
2,β 2
ψ∈D[elβ0(ϕ)
3,β
The model obtained by replacing the original Delaunay graph by the β-Delaunay subgraph of
order β in (1) satisfies the previous sufficient conditions of Preston’s Theorem. From now on,
0
this model is called the β-Delaunay model. In this spirit,some other models may be defined (see
Pseudo-likelihood for some Gibbs point processes 7
e.g. Bertin et al. (1999b), Bertin et al. (1999a)) but we advice the reader to keep in mind the
β-Delaunay model as the main example in order to illustrate the statistical results developped
in this work.
TheframeworkofthispaperisrestrictedtostationaryGibbspointprocessesbasedonenergy
functionrelatedtosomegraph,denotedG (ϕ)forsomefiniteconfigurationϕ(G (ϕ)representing
2 k
the set of cliques of order k), of the form
Kmax
V (ϕ) = u(k)(ξ;ϕ) (2)
Xk=1 ξ∈XGk(ϕ)
Kmax
= θ(1)|ϕ|+ u(k)(ξ;ϕ) , when u(1) ≡ θ(1)
Xk=2 ξ∈XGk(ϕ)
and satisfying Assumptions E , Eloc or more generally Eqloc, E defined by :
1 2 2 3
E V (·) is invariant by translation.
1
Eloc Locality of the local energy : there exists some fixed range denoted by D such that for
2
any ϕ∈ Ω one has
V (0|ϕ)= V (0|ϕ∩B(0,D)).
qloc
E Quasi-locality of the local energy : there exists a nonnegative function ε vanishing
2
asymptotically such that for any ϕ ∈ Ω one has
|V (0|ϕ)−V (0|ϕ∩B(0,D))| < ε(D).
E Stability of the local energy : there exists K ≥ 0 such that for any ϕ ∈ Ω,
3
V (0|ϕ) ≥ −K.
This framework includes some classical point processes such as :
– models based on the usual complete graph G (ϕ) = P (ϕ)) with pairwise interaction
2 2
function satisfying a hard-core or inhibition condition and with finite range.
– k-nearest neighbours models with pairwise interaction function bounded and with finite
range (see Bertin et al. (1999c)).
– Widom-Rowlinson or area interaction model.
– ...
3 Statistical model and inference method
3.1 Statistical model
We consider Gibbs point processes with energy function V (·;θ) parametrized as follows
8 J.-M. Billiot, J.-F. Coeurjolly and R. Drouilhet
As a statistical model weconsider a parametrized versionof (2) wherethedifferentu(k)(ξ,ϕ)
depend on a vector parameters θ(k) and then denoted from now byu(k)(ξ;ϕ,θ(k)). It is assumed
that the vector of parameters θ = (θ ,...,θ ) = θ(1),θ(2),...,θ(Kmax) ∈ Θ where Θ is an
1 p+1
open bounded set of IRp+1. (cid:16) (cid:17)
Our data consist in the realization of a point process with energy function V (·;θ⋆) in a do-
main Λ ⊂ IRd satisfying Assumptions E to E . Thus,θ⋆ is the true parameter to be estimated.
1 3
The Gibbs measure will be denoted by Pθ⋆. From (2), we obtain easily the energy to insert a
point x in a configuration ϕ.
Kmax
V (x|ϕ;θ)= u(k)(ξ;θ(k),ϕ∪{x})− u(k)(ξ;θ(k),ϕ)
Xk=1 ξ∈Gk(ϕ∪X{x})\Gk(ϕ) ξ∈Gk(ϕ)X\Gk(ϕ∪{x})
(3)
Theoretical results presented in the next sections are valid for a general energy function
V (·;θ). But among this class of models, we will focus on energy functions described by (2) and
such that
u(k)(ξ;θ(k),ϕ)= θ(k)Tu(k)(ξ;ϕ).
The energy can be rewritten
Kmax Kmax
V (ϕ;θ) = θ(k)Tu(k)(ξ;ϕ) = θ(k)Tu(k)(ϕ) (4)
Xk=1 ξ∈XGk(ϕ) Xk=1
where for any finite configuration ϕ
u(k)(ϕ)= u(k)(ξ;ϕ) and u(ϕ) = (u (ϕ),...,u (ϕ)) = u(1)(ϕ),...,u(Kmax)(ϕ)
1 p+1
ξ∈XGk(ϕ) (cid:16) (cid:17)
For two finite configurations ϕ and ψ, by denoting
u(k)(ψ|ϕ) = u(k)(ψ∪ϕ)−u(k)(ϕ) and u(ψ|ϕ)= u(1)(ψ|ϕ),...,u(Kmax)(ψ|ϕ) (5)
(cid:16) (cid:17)
we have for any point x
Kmax Kmax
V (ϕ;θ)= θ(k)Tu(k)(ϕ)= θTu(ϕ) and V (x|ϕ;θ)= θ(k)Tu(k)(x|ϕ) = θTu(x|ϕ),
k=1 k=1
X X
(6)
whereu(x|ϕ)= (u (x|ϕ,...,u (x|ϕ) = (u(1)(x|ϕ),...,uKmax(x|ϕ)).Thelocal specification of
1 p+1
theGibbspointprocessassociated toanenergyfunctiondefinedby(6)belongstoanexponential
family.
Pseudo-likelihood for some Gibbs point processes 9
3.2 Pseudo-likelihood
Asprecisedintheintroduction,theideaofmaximumpseudo-likelihoodisduetoBesag(1975)
who first introduced the concept for Markov random fields in order to avoid the normalizing
constant. This work was then widely extended and Jensen and Møller (1991) (Theorem 2.2)
obtained a general expression for Gibbs point processes. With our notation and up to a scalar
factor the pseudo-likelihood defined for a configuration ϕ and a domain of observation Λ is
denoted by PL (ϕ;θ) and given by
Λ
PL (ϕ;θ) = exp − exp(−V (x|ϕ;θ))dx exp(−V (x|ϕ\x;θ)). (7)
Λ
(cid:18) ZΛ (cid:19)x∈ϕ
YΛ
It is more convenient to define (and work with) the log-pseudo-likelihood function, denoted by
LPL (ϕ;θ).
Λ
LPL (ϕ;θ) = − exp(−V (x|ϕ;θ))dx− V (x|ϕ\x;θ) (8)
Λ
ZΛ x∈ϕ
XΛ
3.3 Main statistical tools
Let us start by presenting a particular case of Campbell Theorem combined with Gl¨otz
Theorem that is widely usedin our futureproofs.For some finiteconfiguration ϕ (resp. for some
set G) and for all x, we denote by ϕ (resp. G ) the configuration ϕ (resp. the set G) translated
x x
of x.
Corollary 1 If the probability measure P is stationary and if the function h(·,·) (used in Camp-
bell Theorem) can be decomposed into h(x,ϕ) = 1(x ∈ Λ)g(x,ϕ) for Λ ⊂ IRd where g(·,·) is
such that g(x,ϕ ) = g(0,ϕ) for all x, then the refined Campbell theorem combined with Gl¨otz
x
Theorem allow us to obtain
E g(x,Φ\x) = |Λ| E g(0,Φ)exp(−V (0|Φ)) (9)
P P
(cid:16)x∈XΦΛ\x (cid:17) (cid:16) (cid:17)
Let us now present a version of an ergodic theorem obtained by Nguyen and Zessin (1979)
and widely used in this paper. Let D > 0 and denote by Λ the following fixed domain
0
e D D
Λ = z ∈ IR2,− ≤ |z| ≤ ,
0
2 2
( )
e e
where for all z ∈ IR2, |z| = max(z ,z ).
1 2
Theorem 2 (Nguyen and Zessin (1979)) Let {H ,G ∈ B } be a family of random vari-
G b
ables, which is covariant, that for all x∈ IRd,
H (ϕ ) = H (ϕ), a.s.
Gx x G
10 J.-M. Billiot, J.-F. Coeurjolly and R. Drouilhet
and additive, that is for every disjoint G ,G ∈ B ,
1 2 b
H = H +H , a.s.
G1∪G2 G1 G2
Let I be the sub−σ−algebra of F consisting of translation invariant (with probability 1) sets.
Assume there exists a nonnegative and integrable random variable Y such that |H | ≤ Y a.s.
G
for every convex G⊂ Λ . Then,
0
1 1
lim H = E(H |I), a.s.
n→+∞|Gn| Gn |Λ0| Λ0
for each regular sequence G → IRd.
n
4 Consistency of the maximum pseudo-likelihood estimator
Maximizing the pseudo-likelihood is equivalent to minimize U (θ) defined by
n
1
U (θ)= − LPL (ϕ;θ).
n |Λ | Λn
n
We denote byθ = θ (ϕ) the maximum pseudo-likelihood estimator based on the configuration
n n
ϕ, alternatively defined as
b b
θ (ϕ)= argmin U (θ)
n θ∈Θ n
In this section, the existence ofban ergodic measure is ensured, relatively to our framework,
qloc
by Assumptions E , E and E . The following Assumptions are needed to derive the almost
1 2 3
sure convergence of this estimator.
C (Λ ) is a regular sequence of domains such that Λ → IR2 as n → +∞.
1 n n≥1 n
C For all θ ∈ Θ,
2
V (0|·;θ) ∈L1(Pθ⋆).
C For all θ ∈Θ\θ⋆
3
Pθ⋆ {ϕ, V (0|ϕ;θ) 6= V (0|ϕ;θ⋆)} > 0
(cid:16) (cid:17)
C4 For all θ,θ′ ∈ Θ, there exists c > 0 such that Pθ⋆−almost surely, we have
|V (0|Φ;θ)−V 0|Φ;θ′ |≤ ||θ−θ′||cg(0,Φ) (10)
(cid:0) (cid:1)
where g(·,·) is a function such that for all x, g(0,Φ) = g(x,Φ ) and such that g(0,·) ∈
x
L1(Pθ⋆).
Remark 1 If one only assumes the existence of an ergodic measure, in particular without As-
sumption E (taken into account to express Assumptions C and C ) then
3 2 4
