No. 34
Online ISSN: 1882-2274

COMPUTER SCIENCE MONOGRAPHS
A Publicationof
The Institute of Statistical Mathematics

SIMULATION AND ESTIMATION OF THE
NEYMAN-SCOTT TYPE SPATIAL CLUSTER MODELSbyU. Tanaka, Y. Ogata and K. Katsura

March 2008

Computer Science Monographs

Editors
Junji NANAKO(in chief )Satoshi ITO
Nobuhisa KASHIWAGI
Yumi TAKIZAWA

The series of Computer Science Monographs publishes the statisticalsoftware and the articles on computer applications originating at The
Institute of Statistical Mathematics.

All communications relating to thispublication should be addressed to the Editorial Office, The Institute of
Statistical Mathematics, 4-6-7 Minami-Azabu, Minato-ku, Tokyo 106-8569,
Japan.

The electronic version of this issue is downloadable at
http://www.ism.ac.jp/

Simulation and estimation of the
Neyman-Scott type spatial cluster models
Tanaka, U., Ogata, Y. and Katsura, K.
The Institute of Statistical Mathematics

Contents
1 Overview of the models

2

2 Simulation of the models
3
2.1 Thomas model . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.1.1 Simulation method . . . . . . . . . . . . . . . . . . . . 4
2.1.2 Inputs and outputs of the program [Thomas-Simulation] 4
2.2 Inverse-power type model . . . . . . . . . . . . . . . . . . . . 5
2.2.1 Simulation method . . . . . . . . . . . . . . . . . . . . 7
2.2.2 Inputs and outputs of the program [IP-Simulation] . . 7
2.3 Generalized Thomas model of type A . . . . . . . . . . . . . . 8
2.3.1 Simulation method . . . . . . . . . . . . . . . . . . . . 10
2.3.2 Inputs and outputs of the program [TypeA-Simulation] 10
2.4 Generalized Thomas model of type B . . . . . . . . . . . . . . 11
2.4.1 Simulation method . . . . . . . . . . . . . . . . . . . . 13
2.4.2 Inputs and outputs of the program [TypeB-Simulation] 13
2.5 Generalized Thomas model of type C . . . . . . . . . . . . . . 14
2.5.1 Simulation method . . . . . . . . . . . . . . . . . . . . 14
2.5.2 Inputs and outputs of the program [TypeC-Simulation] 16
3 Parameter estimation
3.1 Overview of the estimation procedure . . . . . . . . . . . . . .
3.1.1 Palm intensity functions . . . . . . . . . . . . . . . . .
3.1.2 Palm log-likelihood functions . . . . . . . . . . . . . .
3.2 Inputs and outputs of the program [Simplex-Thomas] . . . . .
1

17
17
19
20
22

3.3 Non-parametric estimate of the Palm intensity . . . . . . . . .
3.4 Inputs and outputs of the program [Palm-Thomas] . . . . . .
3.5 Inputs and outputs of the program [Simplex-IP] . . . . . . . .
3.6 Inputs and outputs of the program [Palm-IP] . . . . . . . . . .
3.7 Inputs and outputs of the program [Simplex-TypeA] . . . . . .
3.8 Inputs and outputs of the program [Palm-TypeA] . . . . . . .
3.9 Inputs and outputs of the program [Simplex-TypeB] . . . . . .
3.10 Inputs and outputs of the program [Palm-TypeB] . . . . . . .
3.11 Inputs and outputs of the program [Simplex-TypeC] . . . . . .
3.12 Inputs and outputs of the program [Palm-TypeC] . . . . . . .

1

24
25
26
29
30
33
35
37
39
41

Overview of the models

We provide FORTRAN programs for the maximum Palm likelihood estimation of parameters of the Neyman-Scott type spatial cluster model (Tanaka
et al. [7]). Neyman-Scott models are frequently used to apply them to spatial
cluster point patterns. The centers of clusters of point pattern, called parent
points, form a stationary Poisson point ﬁeld with intensity µ. Then, theNeyman-Scott model is deﬁned as the set of all oﬀspring points described asfollows (e.g., Stoyan and Stoyan [5], [6]):1. Each parent point produces a random number M of oﬀspring points
which are realized independently and identically according to the probabilities pm = P { M = m }, m = 0, 1, . . . .
2. Locations of oﬀspring points are distributed independently and identically according to a density function qτ (x, y), where (x, y) is coordinatesof oﬀspring points relative to their parent point, and the parameter τ
indicates a set of coeﬃcients of the density function.
Hereafter, we consider parameters θ = (µ, ν, τ ), where ν = E[M] =mpm ,as the parameters that characterize the various models of the Neyman-Scott
type. In particular, throughout this article, we assume that M has the
Poisson distribution with mean ν.
Furthermore, we consider multi-type cluster models, which are formed
as follows (Tanaka et al. [7]): the models allow two types of clusters withthe diﬀerent cluster sizes and distance scales. The intensities of centres are2µ1 and µ2 , respectively. The corresponding oﬀspring cluster size distributes
according to Poisson random variable with means ν1 and ν2 , respectively.
In order to estimate the parameters and select the competing models,
we have prepared some programs written in FORTRAN77. We also plotted
ﬁgures which are delineated in PDF format.
While such point patterns can be considered in spaces of arbitrary dimension, we here restrict ourself to planar point patterns in observed rectangularregion W = [0, Tx ] × [0, Ty ], assuming the periodic boundary condition. Furthermore, we restrict ourself to the standardized case where Tx = 1 andTy ≥ 1 without loss of generality except for the scaling. All the models of
point patterns are assumed to be stationary and isotropic. Therefore thedensity function qτ (x, y) is denoted by polar coordinate qτ (r) of distance r
from its parent point.
For the simulations, we use the pseudo-random number generation routinereferred to Wichmann and Hill [8].

2Simulation of the models

In order to simulate spatial clustered point pattern of the Neyman-Scottmodels and their extensions, this package includes the FORTRAN programs[Thomas-Simulate.f], [IP-Simulate.f], [TypeA-Simulate.f], [TypeBSimulate.f] and [TypeC-Simulate.f]. We consider the periodic boundary
condition (i.e., a two-dimensional torus) for each of simulations of models on
a rectangular with normalized scale as described below.

2.1

Thomas model

The Thomas model has the density function that is given by the bivariate
Gaussian distribution with zero mean and the covariance matrix σ 2 I where
I is 2 × 2 identity matrix as follows (e.g., Tanaka et al. [7]). Due to isotropy
this is described in polar coordinate as a functionrr2
qσ (r) = 2 exp − 2 , r ≥ 0,σ
2σ
of distance r from the corresponding parent’s location.

32.1.1

Simulation method

The simulation of this model is produced as follows: let random variable U
be independently and uniformly distributed in [0, 1]. We put
 rqσ (t) dt
U=
0


r2
= 1 − exp − 2 .
2σThen we haver=σ
−2 log(1 − U).

Let (xpi , yip), i = 1, 2, . . . , I be a coordinate of each parent point where
the integer I is generated from the Poisson random variable P oisson(µ)
with mean µ from now on. Then, for each i, the number of oﬀspring Ji is
generated by the random variable P oisson(ν) with mean ν. Then, using
series of diﬀerent uniform random numbers {U} for diﬀerent i and j, theoﬀspring coordinates (xij , yji ), j = 1, 2, . . . , Ji is given by
xij = xpi + r cos(2πU),yji = yip + r sin(2πU),owing to the isotropy condition of the distribution.
We adopt the following procedure of generation of the Poisson random
number M = P oisson(ν) (e.g., Devroye [2]): Given a positive number ν andlet a sequence of a random variable {Uk } be independently and uniformly
distributed in [0, 1], M is the smallest integer such thatM+1− log Uk > ν,

(2.1)

k=1

where log represents natural logarithm.
2.1.2

Inputs and outputs of the program [Thomas-Simulation]

In order to execute the FORTRAN program [Thomas-Simulate.f] for thissimulation, we have to prepare the following control ﬁle that speciﬁes the4necessary parameter values and others, as input variables. Then we have
output ﬁles as speciﬁed below.
INPUTS:
The control ﬁle [Thomas.param] includes following variables in the fourlines.
• 1st line. The ﬁrst three positive integer variables ix, iy and iz in thisﬁle are initial seeds for a sequence of uniform random numbers. We
can get diﬀerent independent point patterns by changing these seeds.
• 2nd line. The variable Ty (≥ 1) is provided here for the standardized
coordinates of points in the rectangular region [0, 1] × [0, Ty ] without
loss of generality except for the scaling.
• 3rd line. The parameter values µ, ν and σ are provided in the orderfrom the left to the right.
• 4th line is ignored in this case.
OUTPUTS:• Calculated records of the simulation program are printed on line. For
convenience of users to check, we left the print in the ﬁle [Thomasparents-oﬀspring.print], which shows numbers of oﬀspring for eachparent.
• The ﬁles named by [Thomas-parents.xy] and [Thomas-oﬀspring.xy]
provide the coordinates of simulated parent points and their oﬀspring
points, respectively, in the rectangular region [0, 1] × [0, Ty ].
Here we simulate parent points with µ = 50.0 (cf., Figure 1) and of theThomas oﬀspring with (µ, ν, σ) = (50.0, 30.0, 0.03) (cf., Figure 2).

2.2

Inverse-power type model

The Inverse-power type density function of the distance is given as follows(Tanaka et al. [7]):qp,c (r) =cp−1 (p − 1),
(r + c)pr ≥ 0,
5p > 1,

c > 0,

0.0

0.2

0.4

0.6

0.8

1.0

Parent points

0.0

0.2

0.4

0.6

0.8

1.0

Figure 1: Simulated parent points of the Thomas model with µ = 50.0.

0.0

0.2

0.4

0.6

0.8

1.0

Thomas model

0.0

0.2

0.4

0.6

0.8

1.0

Figure 2: Simulated oﬀspring points of the Thomas model with (µ, ν, σ) =(50.0, 30.0, 0.03).

6

where the parameters p and c represent a decay order and scaling factor withrespect to the distance, respectively.
2.2.1

Simulation method

The simulation of the inverse-power type model is carried out as follows:
Let random variable U be independently and uniformly distributed in [0, 1].
For any r ≥ 0,
 r
Qp,c (r) :=qp,c (t) dt
0(r + c)1−p − c1−p
1−p
p−1
1−p
= 1 − c (r + c) .
= cp−1 (p − 1)

Here, we put Qp,c (r) = U. From this, we have
r = c (1 − U)1/(1−p) − 1 .
Similarly, coordinate of the oﬀspring points (xij , yji ), j = 1, 2, . . . , P oisson(ν)
with the inverse-power type is given for each i = 1, 2, . . . , P oisson(µ),
xij = xpi + r cos(2πU),yji = yip + r sin(2πU),using series of diﬀerent uniform random numbers {U} for diﬀerent i and j.
2.2.2

Inputs and outputs of the program [IP-Simulation]

In order to execute the program [IP-Simulate.f] for this simulation, we haveto prepare the following control ﬁle that speciﬁes the necessary parametervalues and others, as inputs. Then we have outputs of ﬁles as speciﬁed below.
INPUTS:
The control ﬁle [IP.param] includes following variables in the ﬁve lines.
• 1st line. The ﬁrst three positive integer variables ix, iy and iz in thisﬁle are initial seeds for a sequence of uniform random numbers. We
can get diﬀerent independent point patterns by changing these seeds.
7• 2nd line. The variable Ty (≥ 1) is provided here for the standardized
coordinates of points in the rectangular region [0, 1] × [0, Ty ] without
loss of generality except for the scaling.
• 3rd line. The parameter values µ, ν, p and c are provided in the orderfrom the left to the right.
• 4th line is ignored in this case.
• 5th line is ignored in this case.
OUTPUTS:• Calculated records of the simulation program are printed on line. For
convenience of users to check, we left the print in the ﬁle [IP-parentsoﬀspring.print], which shows numbers of oﬀspring for each parent.
• The ﬁles named by [IP-parents.xy] and [IP-oﬀspring.xy] provide
the coordinates of simulated parent points and their oﬀspring points,
respectively, in the rectangular region [0, 1] × [0, Ty ].
Here we simulate parent points with µ = 50.0 (cf., Figure 3) and of
the inverse-power type oﬀspring with (µ, ν, p, c) = (50.0, 30.0, 1.5, 0.005) (cf.,Figure 4).

2.3

Generalized Thomas model of type A

The generalized Thomas model of type A, or the Type A model for short,
is given by the mixed Gaussian density function with σ1 and σ2 as follows(Tanaka el al. [7]):




ar
r2
r2
(1 − a)r
exp −
exp −
qτ (r) =+, τ = (a, σ1 , σ2 ),
σ1 2
2σ1 2
σ2 2
2σ2 2
where the parameter a indicates the mixture ratio of the density function.

8

0.0

0.2

0.4

0.6

0.8

1.0

Parent points

0.0

0.2

0.4

0.6

0.8

1.0

Figure 3: Simulated parent points of the inverse-power type with µ = 50.0.

0.0

0.2

0.4

0.6

0.8

1.0

Inverse−power type

0.0

0.2

0.4

0.6

0.8

1.0

Figure 4: Simulated oﬀspring points of the inverse-power type with(µ, ν, p, c) = (50.0, 30.0, 1.5, 0.005).

9

2.3.1

Simulation method

Parents’ conﬁguration and numbers of the oﬀspring cluster sizes are generated
by the same way as the Thomas model (cf., (2.1)).
Let random variable Uk , k = 1, 2 be independently and uniformly distributed in [0, 1]. Then r satisﬁes as follows:r=
σ1 −2 log(1 − U1 ), U2 ≤ a,
σ2 −2 log(1 − U1 ), otherwise.

Then, by the isotropy condition, for i = 1, 2, . . . , P oisson(µ), coordinateof the oﬀspring points (xij , yji ), j = 1, 2, . . . , P oisson(ν) is given for eachi = 1, . . . , P oisson(µ),
xij = xpi + r cos(2πU),yji = yip + r sin(2πU),using series of diﬀerent uniform random numbers {U1 , U2 , U} for diﬀerent iand j.
2.3.2

Inputs and outputs of the program [TypeA-Simulation]

In order to execute the program [TypeA-Simulate.f] for this simulation,
we have to prepare the following control ﬁle that speciﬁes the necessary
parameter values and others, as inputs. Then we have outputs of ﬁles as
speciﬁed below.
INPUTS:
The control ﬁle [TypeA.param] includes following variables in the ﬁvelines.
• 1st line. The ﬁrst three positive integer variables ix, iy and iz in thisﬁle are initial seeds for a sequence of uniform random numbers. We
can get diﬀerent independent point patterns by changing these seeds.
• 2nd line. The variable Ty (≥ 1) is provided here for the standardized
coordinates of points in the rectangular region [0, 1] × [0, Ty ] without
loss of generality except for the scaling.

10• 3rd line. The parameter values µ, ν, a, σ1 and σ2 are provided in theorder from the left to the right.
• 4th line is ignored in this case.
• 5th line is ignored in this case.
OUTPUTS:• Calculated records of the simulation program are printed on line. For
convenience of users to check, we left the print in the ﬁle [TypeAparents-oﬀspring.print], which shows numbers of oﬀspring for eachparent.
• The ﬁles named by [TypeA-parents.xy] and [TypeA-oﬀspring.xy]
provide the coordinates of simulated parent points and their oﬀspring
points, respectively, in the rectangular region [0, 1] × [0, Ty ].
Here we simulate parent points with µ = 50.0 (cf., Figure 5) and of
the mixed Thomas oﬀspring with (µ, ν, a, σ1 , σ2 ) = (50.0, 30.0, 0.3, 0.005, 0.1)
(cf., Figure 6).

2.4

Generalized Thomas model of type B

Furthermore, we provide the extension of the Neyman-Scott model allowing
two types of clusters with the diﬀerent cluster sizes and the distance scales.
In particular, here, we assume two types of clusters with the diﬀerent
sizes and distributions of the Thomas type. The intensities of the respective
parents (centres) are µ1 and µ2 , respectively. The respective clusters havethe Poisson means of the oﬀspring points ν1 and ν2 . The density functionsof the oﬀspring relative to respective parent are given by




r
r2
r
r2
, qσ2 (r) = 2 exp −,(2.2)
qσ1 (r) = 2 exp −
σ1
2σ1 2
σ2
2σ2 2
where σ1 and σ2 indicate their scales. If the mean sizes of the two typesof clusters are the same, then it is hereafter referred to as the generalizedThomas model of type B, or the Type B model for short (Tanaka et al. [7]).

11

0.0

0.2

0.4

0.6

0.8

1.0

Parent points

0.0

0.2

0.4

0.6

0.8

1.0

Figure 5: Simulated parent points of the Type A model with µ = 50.0.

0.0

0.2

0.4

0.6

0.8

1.0

Type A model

0.0

0.2

0.4

0.6

0.8

1.0

Figure 6: Simulated oﬀspring points of the Type A model with(µ, ν, a, σ1 , σ2 ) = (50.0, 30.0, 0.3, 0.005, 0.1).

12

2.4.1

Simulation method

Consider the two types of the Thomas model with parameters (µ1 , ν, σ1 ) and(µ2 , ν, σ2 ). Parents’ conﬁguration and numbers of the oﬀspring cluster sizes
are generated by the two types of uniformly distributed parents (xki , xki ) withi = 1, 2, . . . , P oisson(µk ) for k = 1, 2, respectively.
Then, using series of diﬀerent uniform random numbers {U} for diﬀerentk,ii and j, the oﬀspring coordinates (xk,ij , yj ) of the parents (k, i) with k = 1, 2
and j = 1, 2, . . . , P oisson(ν) is given bykxk,ij = xi + rk cos(2πU),

yjk,i = yik + rk sin(2πU),
where
rk = σk−2 log(1 − Uk ),k = 1, 2,with diﬀerent random numbers {Uk , U} for diﬀerent k, i and j.
2.4.2

Inputs and outputs of the program [TypeB-Simulation]

In order to execute the program [TypeB-Simulate.f] for this simulation,
we have to prepare the following control ﬁle that speciﬁes the necessary
parameter values and others, as inputs. Then we have outputs of ﬁles as
speciﬁed below.
INPUTS:
The control ﬁle [TypeB.param] includes following variables in the fourlines.
• 1st line. The ﬁrst three positive integer variables ix, iy and iz in thisﬁle are initial seeds for a sequence of uniform random numbers. We
can get diﬀerent independent point patterns by changing these seeds.
• 2nd line. The variable Ty (≥ 1) is provided here for the standardized
coordinates of points in the rectangular region [0, 1] × [0, Ty ] without
loss of generality except for the scaling.
• 3rd line. The parameter values µ1 , µ2 , ν, σ1 and σ2 are provided in
the order from the left to the right.
13• 4th line is ignored in this case.
OUTPUTS:• Calculated records of the simulation program are printed on line. For
convenience of users to check, we left the print in the ﬁle [TypeBparents-oﬀspring.print], which shows numbers of oﬀspring for eachparent.
• The ﬁles named by [TypeB-parents.xy] and [TypeB-oﬀspring.xy]
provide the coordinates of simulated parent points and their oﬀspring
points, respectively, in the rectangular region [0, 1] × [0, Ty ].
Here we simulate parent points with µ1 = 10.0 and µ2 = 40.0 (cf., Figure7) and of the two types of the Thomas oﬀspring with (µ1 , µ2 , ν, σ1 , σ2 ) =(10.0, 40.0, 30.0, 0.01, 0, 03) (cf., Figure 8).
It is possible to change this simulation program for the general case where
there is no restrictions among the parameters µ1 , µ2 , ν1 , ν2 , σ1 and σ2 , but
the estimation by means of the Palm log-likelihood does not identify theparameters (cf., Tanaka et al. [7]).

2.5Generalized Thomas model of type C

If ν1 = ν2 and each density function qσi , i = 1, 2 have the same forms (2.2),
then, in particular, for σ1 and σ2 we may constrain ν1 and ν2 to be σ2,ν2 = ν1
σ1and it is hereafter referred to as the generalized Thomas model of type C or
the Type C model, for short (Tanaka et al. [7]).
2.5.1Simulation method

The simulation method of the Type C model is similar to that of the TypeB model because the Type C model consists of two types of the Thomasmodel with the diﬀerent cluster sizes and the distance scales. Consider thetwo types of the Thomas model with parameters (µ1 , ν1 , σ1 ) and (µ2 , ν2 , σ2 ).
14

0.0

0.2

0.4

0.6

0.8

1.0

Parent points

0.0

0.2

0.4

0.6

0.8

1.0

Figure 7: Simulated two types of parent points of the Type B model with
µ1 = 10.0 and µ2 = 40.0.

0.0

0.2

0.4

0.6

0.8

1.0

Type B model

0.0

0.2

0.4

0.6

0.8

1.0

Figure 8: Simulated two types of oﬀspring points of the Type B model with(µ1 , µ2 , ν, σ1 , σ2 ) = (10.0, 40.0, 30.0, 0.01, 0.03).

15

Parents’ conﬁguration and numbers of the oﬀspring cluster sizes are generated by the two types of uniformly distributed parents (xki , xki ) with i =1, 2, . . . , P oisson(µk ) for k = 1, 2, respectively.
Then, using series of diﬀerent uniform random numbers {U} for diﬀerentk,ii and j, the oﬀspring coordinates (xk,ij , yj ), j = 1, 2, . . . , P oisson(νk ) of theparents (k, i) with k = 1, 2 is given bykxk,ij = xi + rk cos(2πU),

yjk,i = yik + rk sin(2πU),
where
rk = σk−2 log(1 − Uk ),k = 1, 2,with diﬀerent random numbers {Uk , U} for diﬀerent k, i and j.
2.5.2

Inputs and outputs of the program [TypeC-Simulation]

In order to execute the program [TypeC-Simulate.f] for this simulation,
we have to prepare the following control ﬁle that speciﬁes the necessary
parameter values and others, as inputs. Then we have outputs of ﬁles as
speciﬁed below.
INPUTS:
The control ﬁle [TypeC.param] includes following variables in the fourlines.
• 1st line. The ﬁrst three positive integer variables ix, iy and iz in thisﬁle are initial seeds for a sequence of uniform random numbers. We
can get diﬀerent independent point patterns by changing these seeds.
• 2nd line. The variable Ty (≥ 1) is provided here for the standardized
coordinates of points in the rectangular region [0, 1] × [0, Ty ] without
loss of generality except for the scaling.
• 3rd line. The parameter values µ1 , µ2 , ν1 , ν2 , σ1 and σ2 are provided
in the order from the left to the right.
• 4th line is ignored in this case.
16

OUTPUTS:• Calculated records of the simulation program are printed on line. For
convenience of users to check, we left the print in the ﬁle [TypeCparents-oﬀspring.print], which shows numbers of oﬀspring for eachparent.
• The ﬁles named by [TypeC-parents.xy] and [TypeC-oﬀspring.xy]
provide the coordinates of simulated parent points and their oﬀspring
points, respectively, in the rectangular region [0, 1] × [0, Ty ].
Here we simulate parent points with µ1 = 5.0 and µ2 = 9.0 (cf., Figure 9)
and then respective oﬀsprings with (µ1 , µ2 , ν1 , ν2 , σ1 , σ2 ) = (5.0, 9.0, 30.0, 150.0,
0.01, 0.05) (cf., Figure 10).
It is possible to change this simulation program for the general case where
there is no restrictions among the parameters µ1 , µ2 , ν1 , ν2 , σ1 and σ2 , but
the estimation by means of the Palm log-likelihood does not identify theparameters (Tanaka et al. [7]).

3Parameter estimation

3.1

Overview of the estimation procedure

It is very desirable to estimate model parameters as acculate as possible.
However, it is diﬃcult to obtain maximum likelihood estimates due to thestructure and property of Neyman-Scott cluster models (e.g., Tanaka etal. [7]). In order to obtain the better estimates than the contract estimatesusing the diﬀerence between K-functions of the model and of the empirical,
we shall consider a likelihood function based on the Palm intensity function(e.g., Tanaka et al. [7]) below.
Here the parameter estimation due to the Palm likelihood function in (3.6)
and (3.7) is provided for various models. For the maximization procedure(equivalent to minimization without loss of generality), we adopt the Simplex
estimation method (e.g., Kowalik and Osborne [3]) to maximize the Palm
likelihood function (or minimize the negative Palm likelihood function).
The simplex means a usual geometric ﬁgure, to simplify, which is a setof n + 1 points in the n-dimensional Euclidean space. The principal idea of
this method is that we can easily form a new simplex from the current one
17

0.0

0.2

0.4

0.6

0.8

1.0

Parent points

0.0

0.2

0.4

0.6

0.8

1.0

Figure 9: Simulated two types of parent points of the Type C model with
µ1 = 5.0 and µ2 = 9.0.

0.0

0.2

0.4

0.6

0.8

1.0

Type C model

0.0

0.2

0.4

0.6

0.8

1.0

Figure 10: Simulated two types of oﬀspring points of the Type C model with(µ1 , µ2 , ν1 , ν2 , σ1 , σ2 ) = (5.0, 9.0, 30.0, 150.0, 0.01, 0.05).

18

by reﬂecting one point in the hyperplane spanned by the remaining points.
If we choose for this purpose the vertex of the simplex at which the function
is greatest, we can expect that at the reﬂected vertex the function value will
be lower. If this is the case then we can continue the process and move oursimplex closer to the minimum.
In order to calculate the function, we need to give the derived formulaof the Palm intensity function λo due to the Ripley’s K-function (Tanaka etal. [7]) as follows:λo (r) =λ dK(r), r > 0,2πr dr(3.1)

where λ = µν.
For any Neyman-Scott cluster models, the K-function satisﬁes
m[2]Fτ (r), r ≥ 0,
(3.2)
λK(r) = λπr2 +νwhere m[2] = E[M(M − 1)] and Fτ (r) are the second-order factorial moment
and the distribution function of distance between the oﬀspring points inthe same cluster, respectively (cf., Stoyan and Stoyan [5], [6]). If M =P oisson(ν), then m[2] = ν 2 .
The distribution Fτ (r) of the distance between the oﬀspring points of thesame cluster is given by 2 r/2  r+r1
r1 + r2 2 − r 2
1arccosFτ (r) = 2qτ (r2 ) dr2 qτ (r1 ) dr1
2r1 r20r−r1 π 2 ∞  r1 +r
r1 + r2 2 − r 2
1arccos+qτ (r2 ) dr2 qτ (r1 ) dr1
π2r1 r2
r/2
r1
 r/2  r−r1+qτ (r2 ) dr2 qτ (r1 ) dr1 ,
(3.3)
0

r1

for any r ≥ 0 and the parameter τ (e.g., Tanaka et al. [7]).
3.1.1Palm intensity functions

Although analytical form of the Palm intensity function is not available forthe inverse-power type model and the Type A model due to the complexityof the fτ (r) that is the derivative of the Fτ (r), it is available for the Thomasmodel, Type B model and Type C model. According to the Tanaka et al. [7],the Palm intensity functions of the following models are calculated as follows:19

Example 3.1 (Thomas model). (e.g., Daley and Vere-Jones [1]) For any
r ≥ 0,νr2
λo (r) = µν +exp − 2 .
4πσ 2
4σ
Example 3.2 (Type B model). For any r ≥ 0,




a
ν
r2
r2
(1 − a)exp −exp −λo (r) = λ ++,4π σ1 2
4σ1 2
σ2 2
4σ2 2

(3.4)

where λ = ν(µ1 + µ2 ) is the total population size and a = µ1 /(µ1 + µ2 ) is theratio of the parent points of the smaller sized cluster to the total ones.
Example 3.3 (Type C model). For any r ≥ 0,




1 aν1
r2
r2
(1 − a)ν2
exp −
exp −
λo (r) = λ ++,4π σ1 2
4σ1 2
σ2 2
4σ2 2

(3.5)

where λ = µ1 ν1 + µ2ν2 is the total population size and a = µ1 ν1 /λ is the ratioof the all oﬀspring points of smaller sized cluster to the total population size.
3.1.2

Palm log-likelihood functions

Procedures of parameter estimation due to the Palm intensity functions depend on models.
Thus, the Palm log-likelihood function with the analytical Palm intenstiyfunction is given as follows (e.g., Tanaka et al. [7]): Rlog L(µ, ν, τ ) =log(N (W )λo (rij )) − N(W )2πλo (r)r dr,0{i,j;i=j,rij <R}(3.6)
where N(W ) is the number of points in W , and the sum is taken for any
pair of points xi and xj , i = j such that the distance rij between xi and xj issmaller than R, where R is suﬃciently greater than the minimum value of rsuch that λo (r) is equal to the intensity λ = µν but assume R ≤ 1/2 under
the periodic boundary condition.
Here we note that R is shown from a non-parametric Palm intensity
function as below according to the Tanaka et al. [7], where R is taken as 1/2
20

at the largest. Also, we note that the constant log N(W ) of the ﬁrst term
can be neglected in maximization of (3.6).
The followings analytic form of the log Palm-likelihood are provided for
the Thomas, Type B and Type C models (cf., Tanaka et al. [7]).
Example 3.4 (Thomas model). The Palm log-likelihood function of theThomas model is analytically calculated because the Thomas model has theanalytical form of the Palm intensity function. By (3.1), (3.2) and (3.6) or(3.7), we have



1
rij 2
log L(µ, ν, σ) =log ν µ +exp − 24πσ 24σ{i,j;i=j,rij <R}R22− N(W )ν πµR + 1 − exp − 2,
4σwith R = 1/2 which means the half of the tx (TX) in side length of thenormalized rectangular.
Example 3.5 (Type B model). The Palm log-likelihood function of theType B model is analytically calculated because the Type B model has theanalytical form of the Palm intensity function. By (3.4) and (3.6), we havelog L(λ, α, β, σ1 , σ2 )
1log λ +=
4π
{i,j;i=j,rij <R}





rij 2
rij 2
β
α
exp −
+ 2 exp −
σ1 2
4σ1 2
σ2
4σ2 2





R2
R2
+ β 1 − exp −− N(W ) πR λ + α 1 − exp −
4σ1 2
4σ2 2
2

,

with R = 1/2, where α = aν and β = (1 − a)ν.
Example 3.6 (Type C model). The Palm log-likelihood function of theType C model is analytically calculated because the Type C model has the21

analytical form of the Palm intensity function. By (3.5) and (3.6), we havelog L(λ, α, β, σ1 , σ2 )
1=log λ +
4π
{i,j;i=j,rij <R}





α
rij 2
rij 2
β
exp −
+ 2 exp −
σ1 2
4σ1 2
σ2
4σ2 2





R2
R2
− N(W ) πR λ + α 1 − exp −+ β 1 − exp −
4σ1 2
4σ2 2
2

,

with R = 1/2, where α = aν1 and β = (1 − a)ν2 .
On the other hand, general model including the inverse-power type andType A model such that Palm intensity function can not be calculated analytically due to the complexity of the diﬀerential of Fτ (r) require the numericalcalculation of the Palm likelihood function. For such models, we need totake the alternative form without explicit representation of the Palm intensity function (Tanaka et al. [7]) as follows: By (3.1), the representation of
the function using the K-function is written bylog L(µ, ν, τ ) =log{N (W )λo (rij )} − N(W )µνK(R),(3.7){i,j;i=j,rij <R}

where we calculate (3.1), (3.2) and (3.3) numerically.
To summarize, the programs [Simplex-Thomas.f], [Simplex-TypeB.f]
and [Simplex-TypeC.f] of the Thomas related models use analytical formsof the functions, while the programs [Simplex-IP.f] and [Simplex-TypeA.f]of other type models have to use numerical integration and diﬀerence tocompute the functions (see below), which need very long c.p.u. time in theminimization procedure.

3.2

Inputs and outputs of the program [Simplex-Thomas]

In order to execute the program [Simplex-Thomas.f] for a dataset, we haveto prepare the following control ﬁle that speciﬁes the necessary parametervalues, name of data and other variables, as inputs. Then we have outputs
of ﬁles as speciﬁed below.
INPUTS:
22

The following variables are given in the four lines of the control ﬁle[Thomas.param] that is also used for the simulation program of the Thomasmodel.
• 1st line. All the variables in this line are ignored in the present case.
• 2nd line. The variable Ty (≥ 1) is provided here for the standardized
coordinates of points in the rectangular region [0, 1] × [0, Ty ] without
loss of generality except for the scaling. The coordinates of points in
the rectangular region are standardized. Note here that the coordinatesof given data are considered the periodic boundary condition.
• 3rd line. The initial guesses of the parameters µ, ν and σ are provided
in order to get the MPLE values.
• 4th line describes the name of dataset that are standardized by therectangular length of the smaller side that should be set as the xcoordinate.
Here we set the true values (µ0 , ν0 , σ0 ) = (50.0, 30.0, 0.03) to apply thesimulated data [Thomas-oﬀspring.xy] (see Figure2).
OUTPUTS:The output ﬁle [Thomas.simplex.print] provides the process of minimizing the negative Palm log-likelihood function until the values converge to
the MPLE values for the given data. The outputs are arrayed as follows:• In the 1st column, characters update or testfn is written, where theupdate indicates that − log L value in the 2nd column has attained
the minimum so far, otherwise not.
• The 3rd column and the rests show the values of the corresponding
parameters as indicated in the 1st row.
• From the 2nd row to the last row show the optimization procedure
starting from the initial values in the 2nd row. The last row shows
the minimized − log L and the corresponding MPLE values.
The other output ﬁle [Thomas.simplex.x.print] shows the same optimization procedure by the simplex method as indicated in above output ﬁle23

but with the normalized parameters depending on the initial scaling guess ofthe parameters. Namely, the actual estimates are obtained by the indicatedx-values times the corresponding initial values.
• The 1st column shows the smallest − log L value at the stage so far.
• The 2nd column shows the standard deviations among function valuesof the vertices x(i) of the simplex.
• The 3rd column and the rests show the values of the corresponding
normalized parameters x(i) as indicated in the 1st row.
• From the 2nd row to the last row show the optimization procedure
starting from the initial values in the 2nd row, all of which shouldbe 1.0. The last row shows the minimized − log L and the normalized variables corresponding to the MPLE values relative to the initial
estimates. The 2nd column is also expected to be small enough.

3.3

Non-parametric estimate of the Palm intensity

The programs [Palm-Thomas.f], [Palm-IP.f], [Palm-TypeA.f], [PalmTypeB.f] and [Palm-TypeC.f] below include the subroutine that calculates the non-parametric Palm intensity function estimated directly from a
set of point pattern data.
The non-parametric Palm intensity function is realized as follows: let
ψ1 , . . . , ψn be the Cartesian vector coordinates of given points. Then thediﬀerences of the vector coordinates
Ψi,j ≡ ψj − ψi ,i = 1, . . . , n,

j = 1, . . . , n,i = j,is regarded as a realized conﬁguration under the Palm intensity function, theprocedure of which is the so-called diﬀerence process.
As a non-parametric estimation of the Palm intensity function λ0 (r) in
the polar coordinate, we count the number N(r, r + ∆) of such vectors ineach annular set whose radii ends are r and r +∆. Namely, roughly speaking,the Palm intensity function notices not conﬁgurations but diﬀerences of thevector coordinates of the dataset.
Then the ratio of this number to the area of the annulus, that is, N(r, r +∆)/N(W )/(2πr∆), provides an estimate of the Palm intensity function of r
in the polar coordinates.
24

3.4

Inputs and outputs of the program [Palm-Thomas]

INPUTS:
The control ﬁle [Thomasparam.palm] speciﬁes the followings in orderto delineate non-parametric and parametric Palm intensity functions thatare corresponding to the speciﬁed values.
• 1st line describes the dataset name of point pattern to delineate thenon-parametric Palm intensity function.
• 2nd line. A value of ∆ shows a width for the non-parametric Palm
intensity function. The value depends on the number of given point
pattern data. Here we set ∆ = 0.001. The variable Ty (≥ 1) is provided
here for the standardized coordinates of points in the rectangular region[0, 1] × [0, Ty ] without loss of generality except for the scaling.
• From the 3rd line to the bottom line (not more than 10th lines) describe respective value of the given parameters (µi , νi , σi ), i = 1, 2, . . . , I≤ 8, for respective estimates of the Palm intensity functions of theThomas model. If there are nothing in these lines (3rd-10th lines),
the following outputs provide only the non-parametric Palm intensity
function, namely, only the 1st and 2nd columns.

OUTPUTS:
Then we have the outputs-ﬁle [Palm-Thomas.txt] in the following format.
• The 1st column provides the distance r = j∆, where j = 1, 2, . . . , [R/∆],
where [.] is the Gauss’ symbol and R = 1/2 is given in the program forthe normalized rectangular region for the point pattern.
• The 2nd column provides the corresponding values of the non-parametric
Palm intensity function of r, which is normalized by the total intensity estimate (the mean number of points in a unit area) of the point
pattern data.
• The 3rd-last columns show the normalized Palm intensity functions
λo (r)/λ̂ calculated from the given sets of parameter values (µi, νi , σi ).
25

Here, as an example, we provide the non-parametric Palm intensity function with ∆ = 0.001 and two Palm intensity functions corresponding to thetrue and MPLE parameters, respectively (Figure 11).

3.5

Inputs and outputs of the program [Simplex-IP]

This program computes MPLEs via numerical calculation of the Ripley’s
K-function for the Palm-likelihood function (3.7). For the detail, Tanaka etal. [7].
INPUTS:
The following variables are given in the ﬁve lines of the control ﬁle[IP.param], which is also used for the simulation program of the inversepower type.
• 1st line. The ﬁrst three variables in this line are ignored in the present
case. The last column variable indicates iskip for the fast likelihood
but rough approximation of the initial estimates, where iskip = 1 isset in the ﬁnal procedure. Because the calculation takes a long time
due to large number of repetitions of the numerical integration (seeTanaka et al. [7]), a suitable initial estimate is eﬀective in the minimizing procedure. For example, iskip = 1000 is set to obtain a goodinitial estimate. The iskip calculate the Palm intensity function in the
log-likelihood function for every iskip-th ri,j in the ordered distancesof the pairs i and j, and the smaller iskip value provides the moreaccurate estimate but takes longer time for the computation.
• 2nd line. The variable Ty (≥ 1) is provided here for the standardized
coordinates of points in the rectangular region [0, 1] × [0, Ty ] without
loss of generality except for the scaling. The coordinates of points in
the rectangular region are standardized. Note here that the coordinatesof given data are considered the periodic boundary condition.
• 3rd line. The initial guesses of the parameters µ, ν, p and c areprovided in order to get the MPLEs.
• 4th line describes the name of dataset that are standardized by therectangular length of the smaller side that should be set as the xcoordinate.
26

1.5
1.0
0.5

Palm intensity

2.0

3.0

Thomas model

0.001

0.005

0.020

0.100

0.500

r

Figure 11: Palm intensity functions of the Thomas model due to the parameters in Table 1 in doubly logarithmic scales, dots are the non-parametricestimate, the black and gray curves are of the true and MPLE parameters,
respectively, and the horizontal line shows the normalized intensity of thesimulated point pattern. All the plots are normalized with respect to themean intensity (the mean number of events per unit area).

Table 1: MPLE values of the Thomas model from the simulated data.
Parameter
True value
MPLE

µ
ν
σ
50.0 30.0 0.03
51.3 26.1 0.0291

27

• 5th line indicates upper limit value X2 in place of ∞ in the integralin (3.3), which is provided in such a way that X2 > R/2 where R isfrom (3.6). Here we set X2 = 0.3 (> 1/4).
Here we set the true values (µ0 , ν0 , p0 , c0 ) = (50.0, 30.0, 1.5, 0.005) to apply
the simulated data [IP-oﬀspring.xy] (see Figure4).
OUTPUTS:The output ﬁle [IP.simplex.print] provides the process of minimizing
the negative Palm log-likelihood function until the values converge to the
MPLEs for the given data. The outputs are arrayed as follows:• In the 1st column, characters update or testfn is written, where theupdate indicates that − log L value in the 2nd column has attained
the minimum so far, otherwise not.
• The 3rd column and the rests show the values of the corresponding
parameters as indicated in the 1st row.
• From the 2nd row to the last row show the optimization procedure
starting from the initial values in the 2nd row. The last row shows
the minimized − log L and the corresponding MPLEs.
The other output ﬁle [IP.simplex.x.print] shows the same optimization procedure by the simplex method as indicated in above output ﬁle butwith the normalized parameters depending on the initial scaling guess ofthe parameters. Namely, the actual estimates are obtained by the indicatedx-values times the corresponding initial values.
• The 1st column shows the smallest − log L value at the stage so far.
• The 2nd column shows the standard deviations among function valuesof the vertices x(i) of the simplex.
• The 3rd column and the rests show the values of the corresponding
normalized parameters x(i) as indicated in the 1st row.
• From the 2nd row to the last row show the optimization procedure
starting from the initial values in the 2nd row, all of which should be
1.0. The last row shows the minimized − log L and the normalized
variables corresponding to the MPLEs relative to the initial estimates.
The 2nd column is also expected to be small enough.
28

In the directory, in addition, the output ﬁles [IP1000.simplex.print]
and [IP1000.simplex.x.print] for the case of iskip = 1000 are also provided for the sake of users to check within a shorter c.p.u. time.

3.6

Inputs and outputs of the program [Palm-IP]

The program [Palm-IP.f] calculates the Palm intensity functions due to thenumerical diﬀerential calculation of the Fτ (r) (cf., Tanaka et al. [7]) for given
parameter sets of the inverse-power type associated with the non-parametric
Palm intensity function estimated directly from any point pattern data.
INPUTS:
The following variables are given in the control ﬁle [IPparam.palm].
• 1st line describes the dataset name of point pattern to delineate thenon-parametric Palm intensity function.
• 2nd line. A common value of ∆ shows a width for the non-parametric
Palm intensity function and the increment of numerical diﬀerential calculation of the Fτ (r). The value depends on the number of given point
pattern data. Here we set ∆ = 0.001. The variable Ty (≥ 1) is provided here for the standardized coordinates of points in the rectangularregion [0, 1] × [0, Ty ] without loss of generality except for the scaling.
• 3rd line indicates upper limit value X2 in place of ∞ in the integralin (3.3), which is provided in such a way that X2 > R/2 where R isfrom (3.6). Here we set X2 = 0.3 (> 1/4).
• From the 4th line to the bottom line (not more than 10th lines)
describe respective value of the given parameters (µi , νi , pi , ci ), i =1, 2, . . . , I ≤ 7, for respective estimates of the Palm intensity functionsof the inverse-power type. If there are nothing in these lines (4th10th lines), the following outputs provide only the non-parametric Palm
intensity function, namely, only the 1st and 2nd columns.

OUTPUTS:
Then we have the outputs-ﬁle [Palm-IP.txt] in the following format.

29• The 1st column provides the distance r = j∆, where j = 1, 2, . . . , [R/∆],
where [.] is the Gauss’ symbol and R = 1/2 is given in the program forthe normalized rectangular region for the point pattern.
• The 2nd column provides the corresponding values of the non-parametric
Palm intensity function of r, which is normalized by the total intensity estimate (the mean number of points in a unit area) of the point
pattern data.
• The 3rd-last columns show the normalized Palm intensity functions
λo (r)/λ̂ calculated from the given sets of parameter values (µi , νi , pi , ci ).
Here, as an example, we provide the non-parametric Palm intensity function with ∆ = 0.001 and two Palm intensity functions corresponding to thetrue and MPLE parameters, respectively (Figure 12).

3.7

Inputs and outputs of the program [Simplex-TypeA]

The same as the inverse-power type, this program also computes MPLEs via
the Ripley’s K-function, numerically.
INPUTS:
The following variables are given in the ﬁve lines of the control ﬁle [TypeA.param], which is also used for the simulation program of the Type A
model.
• 1st line. The ﬁrst three variables in this line are ignored in the present
case. The last column variable indicates iskip for the fast likelihood
but rough approximation of the initial estimates, where iskip = 1 isset in the ﬁnal procedure. Because the calculation takes a long time due
to large number of repetitions of the numerical integration (see Tanaka
et al., [7]), a suitable initial estimate is eﬀective in the minimizing procedure. Therefore iskip = 1000, for instance, is set to obtain a goodinitial estimate. The iskip calculate the Palm intensity function in the
log-likelihood function for every iskip-th ri,j in the ordered distancesof the pairs i and j, and the smaller iskip value provides the moreaccurate estimate but takes longer time for the computation.

30

10.0
5.0
2.0
0.5

1.0

Palm intensity

20.0

Inverse−power type

0.001

0.005

0.020

0.100

0.500

r

Figure 12: Palm intensity functions of the inverse-power type due to the parameters in Table 2 in doubly logarithmic scales, dots are the non-parametricestimate. The curves of the true and MPLE parameters are overlapped toeach other, and the horizontal line shows the normalized intensity of the simulated point pattern. All the plots are normalized with respect to the mean
intensity (the mean number of events per unit area).

Table 2: MPLE values of the inverse-power type from the simulated data.
Parameter
True value
MPLE

µ
ν
p
c
50.0 30.0 1.5
0.005
56.5 23.9 1.60 0.00557

31• 2nd line. The variable Ty (≥ 1) is provided here for the standardized
coordinates of points in the rectangular region [0, 1] × [0, Ty ] without
loss of generality except for the scaling. The coordinates of points in
the rectangular region are standardized. Note here that the coordinatesof given data are considered the periodic boundary condition.
• 3rd line. The initial guesses of the parameters µ, ν, a, σ1 and σ2 areprovided in order to get the MPLEs.
• 4th line describes the name of dataset that are standardized by therectangular length of the smaller side that should be set as the xcoordinate.
• 5th line indicates upper limit value X2 in place of ∞ in the integralin (3.3), which is provided in such a way that X2 > R/2 where R isfrom (3.6). Here we set X2 = 0.3 (> 1/4).
Here we set the true values (µ, ν, a, σ1 , σ2 ) = (50.0, 30.0, 0.3, 0.005, 0.1) toapply the simulated data [TypeA-oﬀspring.xy] (see Figure6).
OUTPUTS:The output ﬁle [TypeA.simplex.print] provides the process of minimizing the negative Palm log-likelihood function until the values converge to
the MPLEs for the given data. The outputs are arrayed as follows:• In the 1st column, characters update or testfn is written, where theupdate indicates that − log L value in the 2nd column has attained
the minimum so far, otherwise not.
• The 3rd column and the rests show the values of the corresponding
parameters as indicated in the 1st row.
• From the 2nd row to the last row show the optimization procedure
starting from the initial values in the 2nd row. The last row shows
the minimized − log L and the corresponding MPLEs.
The other output ﬁle [TypeA.simplex.x.print] shows the same optimization procedure by the simplex method as indicated in above output ﬁlebut with the normalized parameters depending on the initial scaling guess ofthe parameters. Namely, the actual estimates are obtained by the indicatedx-values times the corresponding initial values.
32• The 1st column shows the smallest − log L value at the stage so far.
• The 2nd column shows the standard deviations among function valuesof the vertices x(i) of the simplex.
• The 3rd column and the rests show the values of the corresponding
normalized parameters x(i) as indicated in the 1st row.
• From the 2nd row to the last row show the optimization procedure
starting from the initial values in the 2nd row, all of which should be
1.0. The last row shows the minimized − log L and the normalized
variables corresponding to the MPLEs relative to the initial estimates.
The 2nd column is also expected to be small enough.
In the directory, in addition, the output ﬁles [TypeA1000.simplex.print]
and [TypeA1000.simplex.x.print] for the case of iskip = 1000 are also
provided for the sake of users to check within a shorter c.p.u. time.

3.8

Inputs and outputs of the program [Palm-TypeA]

The program [Palm-TypeA.f] calculates the Palm intensity functions dueto the numerical diﬀerential calculation of the Fτ (r) for given parameter
sets of the Type A model associated with the non-parametric Palm intensity
function estimated directly from any point pattern data.
INPUTS:
The following are given in the control ﬁle [TypeAparam.palm].
• 1st line describes the dataset name of point pattern to delineate thenon-parametric Palm intensity function.
• 2nd line. A common value of ∆ shows a width for the non-parametric
Palm intensity function and the increment of numerical diﬀerential calculation of the Fτ (r). The value depends on the number of given point
pattern data. Here we set ∆ = 0.001. The variable Ty (≥ 1) is provided here for the standardized coordinates of points in the rectangularregion [0, 1] × [0, Ty ] without loss of generality except for the scaling.
• 3rd line indicates upper limit value X2 in place of ∞ in the integralin (3.3), which is provided in such a way that X2 > R/2 where R isfrom (3.6). Here we set X2 = 0.3 (> 1/4).
33

5.0
2.0
0.5

1.0

Palm intensity

10.0

20.0

Type A model

0.001

0.005

0.020

0.100

0.500

r

Figure 13: Palm intensity functions of the Type A model due to the parameters in Table 3 in doubly logarithmic scales, dots are the non-parametricestimate, the black and gray curves are of the true and MPLE parameters,
respectively, and the horizontal line shows the normalized intensity of thesimulated point pattern. All the plots are normalized with respect to themean intensity (the mean number of events per unit area).

Table 3: MPLE values of the Type A model from the simulated data.
ParameterµνaTrue value 50.0 30.0 0.3
MPLE56.4 23.6 0.355

34

σ1
0.005
0.00520

σ2
0.1
0.107

• From the 4th line to the bottom line (not more than 10th lines)
describe respective value of the given parameters (µi , νi , ai , σ1,i , σ2,i ),
i = 1, 2, . . . , I ≤ 7, for respective estimates of the Palm intensity functions of the Type A model. If there are nothing in these lines (4th-10thlines), the following outputs provide only the non-parametric Palm intensity function, namely, only the 1st and 2nd columns.

OUTPUTS:
Then we have the outputs-ﬁle [Palm-TypeA.txt] in the following format.
• The 1st column provides the distance r = j∆, where j = 1, 2, . . . , [R/∆],
where [.] is the Gauss’ symbol and R = 1/2 is given in the program forthe normalized rectangular region for the point pattern.
• The 2nd column provides the corresponding values of the non-parametric
Palm intensity function of r, which is normalized by the total intensity estimate (the mean number of points in a unit area) of the point
pattern data.
• The 3rd-last columns show the normalized Palm intensity functions
λo (r)/λ̂ calculated from the given sets of parameter values (µi , νi , ai , σ1,i ,
σ2,i ).
Here, as an example, we provide the non-parametric Palm intensity function with ∆ = 0.001 and two Palm intensity functions corresponding to thetrue and MPLE parameters, respectively (Figure 13).

3.9

Inputs and outputs of the program [Simplex-TypeB]

INPUTS:
The following variables are given in the four lines of the control ﬁle[TypeB.param], which is also used for the simulation program of the TypeB model.
• 1st line. All the variables in this line are ignored in the present case.

35• 2nd line. The variable Ty (≥ 1) is provided here for the standardized
coordinates of points in the rectangular region [0, 1] × [0, Ty ] without
loss of generality except for the scaling. The coordinates of points in
the rectangular region are standardized. Note here that the coordinatesof given data are considered the periodic boundary condition.
• 3rd line. The initial guesses of the parameters µ1 , µ2 , ν, σ1 and σ2
are provided in order to get the MPLEs.
• 4th line describes the name of dataset that are standardized by therectangular length of the smaller side that should be set as the xcoordinate.
Here we set the true values (µ1 , µ2 , ν, σ1 , σ2 ) = (10.0, 40.0, 30.0, 0.01, 0.03)
to apply the simulated data [TypeB-oﬀspring.xy] (see Figure 8).
OUTPUTS:The output ﬁle [TypeB.simplex.print] provides the process of minimizing the negative Palm log-likelihood function until the values converge to
the MPLEs for the given data. The outputs are arrayed as follows:• In the 1st column, characters update or testfn is written, where theupdate indicates that − log L value in the 2nd column has attained
the minimum so far, otherwise not.
• The 3rd column and the rests show the values of the corresponding
parameters as indicated in the 1st row.
• From the 2nd row to the last row show the optimization procedure
starting from the initial values in the 2nd row. The last row shows
the minimized − log L and the corresponding MPLEs
The other output ﬁle [TypeB.simplex.x.print] shows the same optimization procedure by the simplex method as indicated in above output ﬁlebut with the normalized parameters depending on the initial scaling guess ofthe parameters. Namely, the actual estimates are obtained by the indicatedx-values times the corresponding initial values.
• The 1st column shows the smallest − log L value at the stage so far.
36• The 2nd column shows the standard deviations among function valuesof the vertices x(i) of the simplex.
• The 3rd column and the rests show the values of the corresponding
normalized parameters x(i) as indicated in the 1st row.
• From the 2nd row to the last row show the optimization procedure
starting from the initial values in the 2nd row, all of which should be
1.0. The last row shows the minimized − log L and the normalized
variables corresponding to the MPLEs relative to the initial estimates.
The 2nd column is also expected to be small enough.

3.10

Inputs and outputs of the program [Palm-TypeB]

The program [Palm-TypeB.f] calculates the Palm intensity functions forgiven parameter sets of the Type B model associated with the non-parametric
Palm intensity function estimated directly from any point pattern data.
INPUTS:
The following variables are given in the control ﬁle [TypeBparam.palm].
• 1st line describes the dataset name of point pattern to delineate thenon-parametric Palm intensity function.
• 2nd line. A value of ∆ shows a width for the non-parametric Palm
intensity function. The value depends on the number of given point
pattern data. Here we set ∆ = 0.001. The variable Ty (≥ 1) is provided
here for the standardized coordinates of points in the rectangular region[0, 1] × [0, Ty ] without loss of generality except for the scaling.
• From the 3rd line to the bottom line (not more than 10th lines)
describe respective value of the given parameters (µi , νi , ai , σ1,i , σ2,i ),where µi = µ1,i + µ2,i , i = 1, 2, . . . , I ≤ 8, for respective estimates of
the Palm intensity functions of the Type B model. If there are nothing
in these lines (3rd-10th lines), the following outputs provide only the
non-parametric Palm intensity function, namely, only the 1st and 2ndcolumns.

37

2.0
0.5

1.0

Palm intensity

5.0

10.0

Type B model

0.001

0.005

0.020

0.100

0.500

r

Figure 14: Palm intensity functions of the Type B model due to the parameters in Table 4 in doubly logarithmic scales, dots are the non-parametricestimate, the black and gray curves are of the true and MPLE parameters,
respectively, and the horizontal line shows the normalized intensity of thesimulated point pattern. All the plots are normalized with respect to themean intensity (the mean number of events per unit area).

Table 4: MPLE values of the Type B model from the simulated data.
Parameter
True value
MPLE

µ
ν
a
50.0 30.0 0.20
44.1 30.7 0.151

38

σ1
σ2
0.01
0.03
0.00945 0.0268

OUTPUTS:
Then we have the outputs-ﬁle [Palm-TypeB.txt] in the following format.
• The 1st column provides the distance r = j∆, where j = 1, 2, . . . , [R/∆],
where [.] is the Gauss’ symbol and R = 1/2 is given in the program forthe normalized rectangular region for the point pattern.
• The 2nd column provides the corresponding values of the non-parametric
Palm intensity function of r, which is normalized by the total intensity estimate (the mean number of points in a unit area) of the point
pattern data.
• The 3rd-last columns show the normalized Palm intensity functions
λo (r)/λ̂ calculated from the given sets of parameter values (µi , νi , ai , σ1,i ,
σ2,i ).
Here, as an example, we provide the non-parametric Palm intensity function with ∆ = 0.001 and two Palm intensity functions corresponding to thetrue and MPLE parameters, respectively (Figure 14).

3.11

Inputs and outputs of the program [Simplex-TypeC]

INPUTS:
The following variables are given in the four lines of the control ﬁle[TypeC.param], which is also used for the simulation program of the TypeC model.
• 1st line. The variables in this line are ignored in the present case.
• 2nd line. The variable Ty (≥ 1) is provided here for the standardized
coordinates of points in the rectangular region [0, 1] × [0, Ty ] without
loss of generality except for the scaling. The coordinates of points in
the rectangular region are standardized. Note here that the coordinatesof given data are considered the periodic boundary condition.
• 3rd line. The initial guesses of the parameters µ1 , µ2 , ν1 , ν2 , σ1 andσ2 are provided in order to get the MPLEs.
39

• 4th line describes the name of dataset that are standardized by therectangular length of the smaller side that should be set as the xcoordinate.
Here, for the initial guess, we set the true values (µ1 , µ2 , ν1 , ν2 , σ1 , σ2 ) =(5.0, 9.0, 30.0, 150.0, 0.01, 0.05) to apply the simulated data [TypeC-oﬀspring.xy] (see Figure10).
OUTPUTS:The output ﬁle [TypeC.simplex.print] provides the process of minimizing the negative Palm log-likelihood function until the values converge to
the MPLEs for the given data. The outputs are arrayed as follows:• In the 1st column, characters update or testfn is written, where theupdate indicates that − log L value in the 2nd column has attained
the minimum so far, otherwise not.
• The 3rd column and the rests show the values of the corresponding
parameters as indicated in the 1st row.
• From the 2nd row to the last row show the optimization procedure
starting from the initial values in the 2nd row. The last row shows
the minimized − log L and the corresponding MPLEs.
The other output ﬁle [TypeC.simplex.x.print] shows the same optimization procedure by the simplex method as indicated in above output ﬁlebut with the normalized parameters depending on the initial scaling guess ofthe parameters. Namely, the actual estimates are obtained by the indicatedx-values times the corresponding initial values.
• The 1st column shows the smallest − log L value at the stage so far.
• The 2nd column shows the standard deviations among function valuesof the vertices x(i) of the simplex.
• The 3rd column and the rests show the values of the corresponding
normalized parameters x(i) as indicated in the 1st row.

40

• From the 2nd row to the last row show the optimization procedure
starting from the initial values in the 2nd row, all of which should be
1.0. The last row shows the minimized − log L and the normalized
variables corresponding to the MPLEs relative to the initial estimates.
The 2nd column is also expected to be small enough.

3.12

Inputs and outputs of the program [Palm-TypeC]

The program [Palm-TypeC.f] calculates the Palm intensity functions forgiven parameter sets of the Type C model associated with the non-parametric
Palm intensity function estimated directly from any point pattern data.
INPUTS:
The following variables are given in the controle ﬁle [TypeCparam.palm].
• 1st line describes the dataset name of point pattern to delineate thenon-parametric Palm intensity function.
• 2nd line. A value of ∆ shows a width for the non-parametric Palm
intensity function. The value depends on the number of given point
pattern data. Here we set ∆ = 0.001. The variable Ty (≥ 1) is provided
here for the standardized coordinates of points in the rectangular region[0, 1] × [0, Ty ] without loss of generality except for the scaling.
• From the 3rd line to the bottom line (not more than 10th lines)
describe respective value of the given parameters (λi , ν1,i , ai , σ1,i , σ2,i ),where λi = µ1,i ν1,i + µ2,i ν2,i , i = 1, 2, . . . , I ≤ 8, for respective estimatesof the Palm intensity functions of the Type C model. If there arenothing in these lines (3rd-10th lines), the following outputs provide
only the non-parametric Palm intensity function, namely, only the 1st
and 2nd columns.
OUTPUTS:
Then we have the outputs-ﬁle [Palm-TypeC.txt] in the following format.
• The 1st column provides the distance r = j∆, where j = 1, 2, . . . , [R/∆],
where [.] is the Gauss’ symbol and R = 1/2 is given in the program forthe normalized rectangular region for the point pattern.
41• The 2nd column provides the corresponding values of the non-parametric
Palm intensity function of r, which is normalized by the total intensity estimate (the mean number of points in a unit area) of the point
pattern data.
• The 3rd-last columns show the normalized Palm intensity functions
λo (r)/λ̂ calculated from the given sets of parameter values (λi , ν1,i , ai , σ1,i ,
σ2,i ).
Here, as an example, we provide the non-parametric Palm intensity function with ∆ = 0.001 and two Palm intensity functions corresponding to thetrue and MPLE parameters, respectively (Figure 15).

42

2.0
0.5

1.0

Palm intensity

5.0

10.0

Type C model

0.001

0.005

0.020

0.100

0.500

r

Figure 15: Palm intensity functions of the Type C model due to the parameters in Table 5 in doubly logarithmic scales, dots are the non-parametricestimate, the black and gray curves are of the true and MPLE parameters,
respectively, and the horizontal line shows the normalized intensity of thesimulated point pattern. All the plots are normalized with respect to themean intensity (the mean number of events per unit area).

Table 5: MPLE values of the Type C model from the simulated data.
Parameter
True valueMPLE

λ1500.01281.0ν1aσ1
σ230.0 0.10 0.01
0.05
43.1 0.30 0.0149 0.0517

43

Acknowledgement
We are grateful to Jiancang Zhuang for his advices on R commands in drawing the ﬁgures.

References[1] Daley, D.J. and Vere-Jones, D. (2008). An Introduction to the Theory
of Point Processes, Vol. II, General Theory and Structure, 2nd Ed.,
Springer-Verlag, New York.
[2] Devroye, L. (1986). Non-Uniform Random Variate Generation,
Springer-Verlag, New York.
[3] Kowalik, J. and Osborne, M. R. (1968). Methods for unconstrained optimization problems, American elsevier publishing company, Inc. New
York.
[4] Press, W. H., Teukolsky, S. A., Vettering, W. T. and Flannery, B. P.
(1992). Numerical Recipes in Fortran 77, 2nd Ed., Cambridge University
Press, Cambridge.
[5] Stoyan, D., Kendall, W. S. and Mecke, J. (1995). Stochastic Geometry
and its Applications, 2nd Ed., John Wiley & Sons, Chichester.
[6] Stoyan, D. and Stoyan, H. (1994). Fractals, Random Shapes and Point
Fields, John Wiley & Sons, Chichester.
[7] Tanaka, U., Ogata, Y. and Stoyan, D. (2008). Parameter estimation andmodel selection for Neyman-Scott point processes, Biometrical Journal
50, 43–57.
[8] Wichmann, B. A. and Hill, I. D. (1982). An Eﬃcient and Portable
Pseudo-random Number Generator, J. R. Statist. Soc. C, 31, 188–190.

44

The Institute of Statistical Mathematics
4-6-7 Minami-Azabu, Minato-ku
Tokyo 106-8569, Japan