Int. J. Anal. Appl. (2023), 21:51 

 

 

Received: Apr. 24, 2023. 

2020 Mathematics Subject Classification. 05B30, 60G55, 62K99. 

Key words and phrases. design of experiments; computer experiment designs; point process; continuum random cluster 

model; connected component process; Markov point processes; nearest-neighbour Markov point process; Strauss 

process; Markov Chain Monte Carlo (MCMC). 

 

https://doi.org/10.28924/2291-8639-21-2023-51 © 2023 the author(s) 

ISSN: 2291-8639  

1 

 

New Computer Experiment Designs Using Continuum Random Cluster Point Process 

Hichem Elmossaoui*, Nadia Oukid 

Department of Mathematics, Faculty of Science, University Saad Daheb, Blida, Algeria 

*Correspondence: elmossaoui.hichem@yahoo.com 

ABSTRACT. In this paper, we propose a new approach for building computer experiment designs using the continuum random 

cluster point process, also referred to as the connected component Markov point process. Our method involves generating 

designs through the Markov Chain Monte Carlo method (MCMC) and the Random Walk Metropolis Hastings algorithm 

(RWMH algorithm), which can be easily scaled to meet various objectives. We have conducted a comprehensive study on 

the convergence of the Markov chain and compared our approach with existing computer experiment designs. Overall, our 

approach offers a novel and flexible solution for constructing computer experiment designs. 

 

1. INTRODUCTION 

The development of modelling techniques, boosted by increased computing power, has led to 

the development of simulators with challenging complexity. The main difficulties stem from the high 

computing cost of the simulator, and from the size of the problem to deal with. It may become 

impossible to consider the direct use of the simulator with some applications, even after significantly 

reducing the size of the system. Hence, the alternative option would be to use one or more functions 

instead of the simulator. These functions are generally relatively simple, and can be obtained by 

means of approximation or interpolation from computer experiment designs. 

In order to achieve a more thorough exploration of the parameter space and obtain 

information throughout the entire experimental area, we propose a method for constructing computer 

https://doi.org/10.28924/2291-8639-21-2023-51
mailto:elmossaoui.hichem@yahoo.com


2 Int. J. Anal. Appl. (2023), 21:51 

 

experiment designs with points uniformly distributed in the unit hypercube. To accomplish this, we 

employ the Continuum Random Cluster Process (CRCP) ([1], [2], [3], [4]) to simulate the 𝑛 

computer experiments, which make up the computer experiment designs. Although the process is not 

Markovian under Ripley-Kelly [5], it is Markovian with respect to a new neighborhood relationship 

that depends on the configuration 𝑥 (as defined in Definition 2.1) and the relation ∼ 𝑥 (as defined 

in Definition 2.3) for nearest neighbours. 

Franco (2008) [6] introduced computer experiment designs based on Strauss point processes 

that incorporate the concept of interaction between pairs of points. Elmossaoui et al (2020) ([7], 

[8]) proposed computer experiment designs using Marked Strauss point processes that can achieve 

multiple objectives simultaneously. The first objective is related to the distribution of points, while 

the second objective concerns the specification of the marks of those points. In contrast to the 

aforementioned approaches, CRCPs are an alternative that allows for the creation of models with 

more or less regular spatial distributions without constraints on parameters. To generate the designs 

proposed in this work, we will use simulation techniques via MCMC methods and Metropolis-Hastings 

([9], [10]) algorithm. There are several sub-types of Metropolis-Hastings algorithms, depending on 

the chosen transition density 𝑞. In this regard, we will develop the algorithm to generate a Markov 

chain of the random walk type, where the density 𝑞 is centred on the current value of the chain and 

is symmetrical. The proposed value 𝑥𝑡+1  takes the form 𝑥𝑡+1 + , where ε is a random perturbation 

that follows the 𝑞 density distribution and is independent of 𝑥𝑡. 

This paper is organised as follows: section 2 is about some general definitions and notations. 

Section 3 concerns the construction of new computer experiment designs based on the use of CRCP 

by means of MCMC method and Metropolis-Hastings algorithm. Section 4 is about the study of the 

convergence of the algorithm. Finally, in section 5 we have compared our results with the existing 

ones.  

2. PRELIMINARIES AND GENERAL DEFINITIONS 

Let (Ω, ℬ, P) be a probability space that models the random aspects of the experiments. Let 

𝜒 be a nonempty set equipped with a Euclidean distance 𝑑, making it a complete and separable 

metric space. In most cases, 𝜒 will be equal to [0,1]𝑝 (a subset of ℝ𝑝), where 𝑝 is the number of 

continuous factors of interest (𝑝 ≥ 1). We will use μ to denote the Lebesgue measure associated 

with this space, considered with its Borel σ-algebra ℬ. 



3 Int. J. Anal. Appl. (2023), 21:51 

 

Definition 2.1. A completion 𝑥 of a point process 𝑋 on 𝜒 is defined as any locally finite collection of 

points from 𝜒, 𝑥 = ( 𝑥1, 𝑥2, … , 𝑥𝑛 ), with 𝑥𝑖 ∈ (𝜒, 𝑑) and 𝑖 ∈ ℕ. In other words, it is a part  𝑥 ⊂ 𝜒 , so 

that 𝑥 ∩ B is finite for every part B of ℬ that is Borelean and limited.  

Let 𝑁𝑙𝑓 denote the set of the locally finite point configurations, x, y … are configurations 

on 𝜒 , and 𝑥𝑖 , 𝑦𝑖 …  points from these configurations. 

Definition 2.2. A point process on 𝜒  is an application  𝑋 of a probability space (Ω, ℬ, P) within the 

set of the locally finite point configurations 𝑁𝑙𝑓, so that for every limited Borelean B, the number of 

points 𝑁(B) = 𝑁𝑋 (B) of points of  𝑋 falling on B is a finite random discrete variable. 

In this definition, 𝜒 can be replaced with a general complete metric space. However, it is 

important to note that the implementation of a point process is, at most, countable and without 

accumulation points. If 𝜒 is bounded, then  𝑁𝑋 (𝜒)  is almost surely finite and the point process is 

said to be finite. We shall consider here only simple point processes that do not allow the repetition 

of points, in which case the realization 𝑥 of the point process coincides with a subset of 𝜒. 

Definition 2.3. Let ∼ be a binary relation that is symmetrical and reflexive on χ. Two points 𝑥𝑖 and 

𝑥𝑖′ are said to be neighbours if  𝑥𝑖 ∼ 𝑥𝑖′. The neighbourhood of  𝑦 ⊂ 𝑥 is given by: 

∂(𝑦|𝑥) = {𝑥𝑖 ∈  𝑥 𝑎𝑛𝑑 𝑥𝑖 ∉  𝑦 ∶ ∃ 𝑥𝑖′  ∈  𝑦 𝑠𝑜 𝑡ℎ𝑎𝑡  𝑥𝑖′  ∼ 𝑥𝑖 } 

For example, 𝑥𝑖 ∼ 𝑥𝑖′ if  𝑑(𝑥𝑖 , 𝑥𝑖′ ) ≤ 𝑟 is the r-neighbourhood relation on (𝜒, 𝑑),                  and 

𝑟 > 0 being a fixed radius. 

Definition 2.4. Let 𝐵(𝑥) = {⋃ 𝐵(𝑥𝑖 ,
𝑟

2
)𝑥𝑖∈𝑥 } ∩ 𝜒 be the reunion of balls centered at points 𝑥𝑖 of  𝑥, of 

radius  
𝑟

2
, and limited at 𝜒. Two points 𝑥𝑖 and 𝑥𝑖′ of  𝑥 are said to be connected for 𝑥 if  𝑥𝑖  and 𝑥𝑖′  

are in the same connected component of 𝐵(𝑥). Such relation shall be noted  𝑥𝑖  ~𝑥  𝑥𝑖′ . 

Definition 2.5. Let 𝑋 be a point process of density 𝜋 in relation to a Poisson point process of law 

𝜋𝜈 (∙) and intensity 𝜈(∙). The process 𝑋 is the nearest-neighbour Markov point processes with regard 

to the relation ~𝑥 if, for every configuration  𝑥 ∈ 𝑁
𝑙𝑓, where 𝜋(𝑥) > 0, then we have what follows:   

(i) 𝜋(𝑦) > 0 for every  𝑦 ⊂ 𝑥 . 

(ii) ∀ 𝑥𝑖 ∈ 𝜒, then : 𝜆(𝑥𝑖 , 𝑥) =
𝑓((𝑥∪{𝑥𝑖 })

𝑓(𝑥)
 depends only on 𝑥𝑖,  ∂({𝑥𝑖 }|𝑥⋃{𝑥𝑖 }) ∩ 𝑥, and the 

two relations ~𝑥 and ~𝑥⋃{𝑥𝑖} . 

The quotient 𝜆(𝑥𝑖 , 𝑥) being the Papangelou conditional density. 



4 Int. J. Anal. Appl. (2023), 21:51 

 

3. COMPUTER EXPERIMENT DESIGN USING NEAREST-NEIGHBOUR MARKOV POINT 

PROCESS 

Identically to the way we proceeded in ([7], [8]), we consider each experiment to be a point 

or particle defined on the interval [0,1]𝑝 , and each configuration 𝑥 as an experiment design. We 

thus simulate 𝑛 experiments to implement a continuum random cluster process. It is worth noting 

that the continuum random cluster process has interaction potentials. These interactions are defined 

by the neighbourhood properties as given by a Markov chain on the nearest neighbours [3]. The 

interaction potential used is the connectedness interaction. These object processes are important in 

modelling repulsive phenomena. The probability density of the process is given as: 

𝜋(𝑥) = 𝑘 𝛽𝑛(𝑥)𝛾ℎ(𝑥) 

where k is a positive normalization constant which makes 𝜋 a density, 𝛽 a positive scaling parameter, 

𝑛(𝑥) denotes the number of points of the configuration 𝑥, 𝛾 is a repulsion parameter such as 0 <

𝛾 < 1, and ℎ(𝑥) = −𝑎(𝑥), where 𝑎(𝑥) refers to the area of the reunion of balls 𝐵(𝑥), or  ℎ(𝑥) =

−𝑐(𝑥), with 𝑐(𝑥) referring to the number of connected components of 𝐵(𝑥). 

Through this article, we aim to study the connected component process with a probability 

density: 

𝜋(𝑥) = 𝑘 𝛽𝑛(𝑥)𝛾−𝑐(𝑥) 

We note that the connected component process is not a Markov process under Ripley and 

Kelly [5]. Actually, two points of 𝜒 can be neighbours with the relation ∼ 𝑥 while being arbitrarily far 

from one another in 𝜒 with regard to the Euclidean distance. Connectedness through connected 

components shall link two points in case there exists a string of balls of radius 𝑟, centred at points 

𝑥𝑖  of 𝑥, and joining with one another. The Papangelou conditional intensity function of the process 

with connectedness interaction is given as follows [11]: 

𝜆(𝑥𝑖 , 𝑥) = 𝛽𝛾
𝑐(𝑥)−𝑐(𝑥∪{𝑥𝑖}) 

Nevertheless, 𝜆(𝑥𝑖 , 𝑥) can depend on a point 𝑥𝑖′  ∈  𝑥  arbitrarily far from 𝑥𝑖 in 𝜒, and we can 

select a configuration 𝑦 so that, for every 𝑥 =  𝑦 ∪ {𝑥𝑖′ }, 𝑐(𝑦) = 2 and 𝑐(𝑥) =  𝑐(𝑥 ∪ {𝑥𝑖 }) = 1. 

Hence, for every r > 0, the process is not Markovian under Ripley-Kelly for the usual r-neighbourhood 

relation. However, if 𝑥𝑖′  ∈  𝑥  is not connected to 𝑥𝑖 in 𝐵(𝑥 ∪ {𝑥𝑖  }),  then 𝑥𝑖′  will not contribute to 

(3.1) 

(3.2) 



5 Int. J. Anal. Appl. (2023), 21:51 

 

the difference 𝑐(𝑥) − 𝑐(𝑥 ∪ {𝑥𝑖 }), and 𝜆(𝑥𝑖 , 𝑥) will not depend on 𝑥𝑖′ , in this case the process is 

nearest-neighbour Markov point process for the nearest neighbour ~𝑥 relation. 

3.1. Generating computer experiment designs via the connected component processes using 

MCMC method and RWMH algorithm 

MCMC method allows the simulation of a 𝜋 density using an ergodic Markov 

chain {𝑋0, 𝑋1, … , 𝑋𝑁𝑀𝐶𝑀𝐶 }; the former being a stationary distribution. To establish such a chain, we 

use RWMH algorithm. The basic idea of the algorithm is to build a transition 𝑃𝑀𝐻  , which is 𝜋-

reversible.  𝑃𝑀𝐻   shall, then, be 𝜋-invariant. This shall be done in two steps:  

• Change proposition transition: we start by proposing a change x → y according to a density 

𝑞(𝑥,∙). 

• Change acceptance probability: change is, then, accepted with the probability 𝑎(𝑥, 𝑦),  a : 

Ω × Ω → [0,1]. 

RWMH algorithm includes an instrumental Markov chain the transition density of which 

depends on the current state; more precisely it satisfies the equation: 𝑞(𝑥, 𝑦) = 𝑞(𝑦 − 𝑥). In addition, 

we can choose an instrumental density 𝑁(𝑥, 𝜎2). So, given a current state 𝑥, the algorithm generates 

a potential state 𝑦 stemming from a normal distribution centered at x and having 𝜎2 as a variance. 

If the new state is accepted the next potential state shall be generated according to a distribution 

𝑁(𝑦, 𝜎2). Else, 𝑥 shall be maintained as current state and another possible state y shall be proposed 

according to a distribution 𝑁(𝑥, 𝜎2). The algorithm, as defined, generates a Markov chain the 

transition probabilities of which are given as follows [12]: 

𝑃𝑀𝐻 (𝑥, 𝑦) = 𝑎(𝑥, 𝑦)𝑞(𝑥, 𝑦) + 1(𝑥 = 𝑦) [1 − ∫ 𝑎(𝑥, 𝑧)𝑞(𝑥, 𝑧)𝑑𝑧
𝛺

] 

𝑃𝑀𝐻 is  𝜋-reversibility if and only if  the balance equation below is satisfied: 

∀𝑥, 𝑦 ∈  𝛺 ∶  𝜋(𝑥) × 𝑞(𝑥, 𝑦) × 𝑎(𝑥, 𝑦)  =  𝜋(𝑦) × 𝑞(𝑦, 𝑥) × 𝑎(𝑦, 𝑥) 

For a couple (𝑥, 𝑦), let’s define Metropolis-Hastings acceptance probability as :  

𝑎(𝑥, 𝑦)  =
𝜋(𝑦) × 𝑞(𝑦, 𝑥)

𝜋(𝑥) × 𝑞(𝑥, 𝑦)
 

The RWMH algorithm with the distribution 𝑁(𝑥, 𝜎2) has an additional property conferred by 

the symmetry of the instrumental density, mainly due to the fact that 𝑞(𝑥, 𝑦) =  𝑞(𝑦, 𝑥). This 

characteristic simplifies the calculation of a: 



6 Int. J. Anal. Appl. (2023), 21:51 

 

𝑎(𝑥, 𝑦) =  
𝜋(𝑦)

𝜋(𝑥)
=

𝑘𝛽𝑛(𝑦)𝛾−𝑐(𝑦)

𝑘𝛽𝑛(𝑥)𝛾−𝑐(𝑥)
= 𝛾𝑐(𝑥)−𝑐(𝑦) 

3.2. Algorithm for the construction of proposed experiment designs  

The computer experiment designs proposed in this work, referred to as connected component 

designs (CCD), are generated through the algorithm presented below, which is a variation of the 

random walk Metropolis Hastings algorithm. 

Algorithm 

• Initialisation  

– Choose an initial experiment design 𝑥 = ( 𝑥1, 𝑥2, … , 𝑥𝑛 ) according to a given 

probability distribution. For instance the normal probability distribution. 

– Take 𝑋0 = 𝑥.  

• For 𝑁 = 1, 2, . . . , 𝑁𝑀𝐶𝑀𝐶 

      For 𝑘 =, …,  𝑛  

– Randomly choose a spin s uniformly on {1, … , 𝑛} 

– Simulate an experiment 𝑦𝑗 according to a proposed distribution q~𝑁(𝑥𝑠, 𝜎
2). Thus, 

we take as a new configuration: 

𝑦 = ( 𝑥1, 𝑥2, … 𝑥𝑠−1, 𝑦𝑗 , 𝑥𝑠+1 … , 𝑥𝑛 ) 

– Calculate the acceptance probability: 

𝑎(𝑥, 𝑦) =  𝑀𝑖𝑛 (1, 𝛾𝑐(𝑥)−𝑐(𝑦)) 

 

– Take  𝑥 = {
𝑦 with a probability      𝑎

    𝑥 with a probability  1 −  𝑎
  

        End for k 

 

                 Take 𝑋𝑁 = 𝑥.  

            End for 𝑁 

For N=1000, figure 1 shows the convergence towards a configuration that satisfies the 

connected component property starting from an initial configuration of 30 points, drawn from a 

normal distribution with mean 0 and variance 1. 



7 Int. J. Anal. Appl. (2023), 21:51 

 

 

Figure 1. Left, an initial configuration of 30 points with c(x)=6. Right, a final configuration 

with c(x)=30 ( γ =0,01, r=0,19 and 𝜎 = 0,05 ). 

 

3.3. Influence of parameters  

The figure below shows the influence of the parameter 𝑟 on the final distribution. The choice 

of the radius proves to be important. A radius that is too small generates a distribution without 

interaction, but with numerous deficiencies. However, a radius that is too large leads to a distribution 

with clusters. 

Figure 2. Left, a configuration of 30 points γ =0.01 and r =0.1. Right, a configuration for γ=0.01 

and r =0.3 



8 Int. J. Anal. Appl. (2023), 21:51 

 

The same as with the interaction radius, it is important to adequately set the attraction 

parameter γ. The figure below shows that it is easier to generate a distribution that meets the 

criteria of filling the space with a strong attraction parameter.  

Figure 3. Left, a configuration of 30 points γ =0.01 and r =0.19. Right, a configuration for γ =0.1 

and r =0.19 

 

4. STUDY OF CONVERGENCE 

In this section, we shall prove the convergence of the sequence of computer experiment 

designs {𝑋0, 𝑋1, … , 𝑋𝑁𝑀𝐶𝑀𝐶 }  generated with the construction algorithm previously introduced towards 

the invariant distribution π of a connected component process. This sequence is the realization of a 

Markov chain with transition kernel: 

𝑃(𝑥, 𝑦) = 𝑃𝑀𝐻
𝑛 (𝑥, 𝑦) 

Moreover, it is important to know whether the distribution of the last generated design 

𝑋𝑁𝑀𝐶𝑀𝐶  is close to the distribution 𝜋. To this end, let us state the main results that are of interest 

to us here:  

Proposition 4.1. On a finite space, the transition kernel 𝑃 = 𝑃𝑀𝐻
𝑛  of the Markov chain (𝑋𝑁 )𝑁 ≥ 0  

obtained from the construction algorithm is irreducible and positive recurrent. The distribution 𝜋 is 

the unique stationary distribution of 𝑃; 𝑃 being aperiodic and a primitive kernel. 

Proof of proposition 4.1. Let us prove that the transition mechanism 𝑃𝑀𝐻 meets the three following 

conditions related to the distribution 𝜋 of the connected component process, defined in (3.2):  𝜋 − 

reversibility, 𝜋 −stationary and 𝜋 −irreductibility. 



9 Int. J. Anal. Appl. (2023), 21:51 

 

The chain is said to be reversible with respect to the target distribution 𝜋(∙) if its transition 

kernel satisfies: 

∀𝑥, 𝑦 ∈ 𝛺, B ∈  ℬ: ∫ 1𝐵(𝑥,𝑦)𝜋(𝑥)𝑃𝑀𝐻 (𝑥, 𝑦)𝑑𝑥
𝛺

= ∫ 1𝐵(𝑥,𝑦)𝜋(𝑦)𝑃𝑀𝐻 (𝑦, 𝑥)𝑑𝑦
𝛺

 

Let 𝑥, 𝑦 ∈  Ω and B ∈ ℬ. Then we have: 

∫ 1𝐵(𝑥,𝑦)𝜋(𝑥)𝑃𝑀𝐻 (𝑥, 𝑦)𝑑𝑥
𝛺

= ∫ 1𝐵(𝑥,𝑦)𝜋(𝑥)𝑎(𝑥, 𝑦)𝑞(𝑥, 𝑦)𝑑𝑥
𝛺

 

+ ∫ 1𝐵(𝑥,𝑦)𝜋(𝑥)
𝛺

[∫ 1 − 𝑎(𝑥, 𝑧)𝑞(𝑥, 𝑧)𝑑𝑧
𝛺

] 𝛿𝑥 (𝑦)𝑑𝑥 

                        = ∫ 1𝐵(𝑥,𝑦)𝜋(𝑥)𝑎(𝑥, 𝑦)𝑞(𝑥, 𝑦)𝑑𝑥
𝛺

+ ∫ 1𝐵(𝑥,𝑥)𝜋(𝑥)
𝛺

[∫ 1 − 𝑎(𝑥, 𝑧)𝑞(𝑥, 𝑧)𝑑𝑧
𝛺

] 𝑑𝑥 

By construction 𝑎 and 𝑞 satisfy,  

𝜋(𝑥)𝑎(𝑥, 𝑦)𝑞(𝑥, 𝑦) = 𝑘𝛽𝑛(𝑥)𝛾−𝑐(𝑥)𝑚𝑖𝑛 {1, 𝛾𝑐(𝑥)−𝑐(𝑦)}𝑞(𝑥, 𝑦) 

                                       = 𝑘 𝑚𝑖𝑛 {𝛽𝑛(𝑥)𝛾−𝑐(𝑥), 𝛽𝑛(𝑥)𝛾−𝑐(𝑦)}𝑞(𝑥, 𝑦)  

And since  𝑛(𝑥) = 𝑛(𝑦) and  𝑞(𝑥, 𝑦) = 𝑞(𝑦, 𝑥). Then we have: 

𝜋(𝑥)𝑎(𝑥, 𝑦)𝑞(𝑥, 𝑦) = 𝑘 min{𝛽𝑛(𝑥)𝛾−𝑐(𝑥), 𝛽𝑛(𝑦)𝛾−𝑐(𝑦)} 𝑞(𝑦, 𝑥) 

                                                                 =𝑘𝛽𝑛(𝑦)𝛾−𝑐(𝑦)𝑚𝑖𝑛 {𝛾𝑐(𝑦)−𝑐(𝑥), 1}𝑞(𝑦, 𝑥) 

 = 𝜋(𝑦)𝑎(𝑦, 𝑥)𝑞(𝑦, 𝑥). 

Finally, 

 ∫ 1𝐵(𝑥,𝑦)𝜋(𝑥)𝑃𝑀𝐻 (𝑥, 𝑦)𝑑𝑥𝛺 = ∫ 1𝐵(𝑥,𝑦)𝜋(𝑦)𝑎(𝑦, 𝑥)𝑞(𝑦, 𝑥)𝑑𝑥𝛺  

                                       + ∫ 1𝐵(𝑦,𝑦)𝜋(𝑦)
𝛺

[∫ 1 − 𝑎(𝑦, 𝑧)𝑞(𝑦, 𝑧)𝑑𝑧
𝛺

] 𝑑𝑦 

           = ∫ 1𝐵(𝑥,𝑦)𝜋(𝑦)𝑃𝑀𝐻 (𝑦, 𝑥)𝑑𝑦.
𝛺

 

This is the condition of 𝜋 − reversibility of the transition mechanism 𝑃𝑀𝐻. 

 A measure 𝜋(∙) is said to be stationary for the transition kernel 𝑃𝑀𝐻 if: 

∀𝑥, 𝑦 ∈ 𝛺;  B, A ∈ ℬ: ∫ 1𝐵(𝑥,𝑦)𝜋(𝑥)𝑃𝑀𝐻 (𝑥, A)𝑑𝑥
𝛺

= ∫ 1𝐵(𝑥,𝑦)𝜋(𝑥)𝑑𝑥
𝛺

 

Let 𝑥 ∈ 𝛺, and B ∈ ℬ. Then for every Borelean A of  ℬ we have:   

 



10 Int. J. Anal. Appl. (2023), 21:51 

 

∫ 1𝐵(𝑥,𝑦)𝜋(𝑥)𝑃𝑀𝐻 (𝑥, A)𝑑𝑥
𝛺

 

         = ∫ 1𝐵(𝑥,𝑦)𝜋(𝑥) [∫ 𝑎(𝑥, 𝑦)𝑞(𝑥, 𝑦)𝑑𝑦
𝛺

] 𝑑𝑥
𝛺

  + ∫ 1𝐵(𝑥,𝑦)𝜋(𝑥)
𝛺

[∫ 1 − 𝑎(𝑥, 𝑧)𝑞(𝑥, 𝑧)𝑑𝑧
𝛺

] 𝛿𝑥 (𝑦)𝑑 

= ∫ 1𝐵(𝑥,𝑦)𝜋(𝑥) [∫ 𝑎(𝑥, 𝑦)𝑞(𝑥, 𝑦)𝑑𝑦
𝛺

] 𝑑𝑥
𝛺

 + ∫ 1𝐵(𝑥,𝑥)𝜋(𝑥)
𝛺

[∫ 1 − 𝑎(𝑥, 𝑧)𝑞(𝑥, 𝑧)𝑑𝑧
𝛺

] 𝑑𝑥 

        = ∫ ∫ 1𝐵(𝑥,𝑦)𝜋(𝑥)𝑎(𝑥, 𝑦)𝑞(𝑥, 𝑦)𝑑𝑦𝑑𝑥
𝛺𝛺

 + ∫ 1𝐵(𝑥,𝑥)𝜋(𝑥)𝑑𝑥
𝛺

− ∫ ∫ 𝜋(𝑥)𝑎(𝑥, 𝑧)𝑞(𝑥, 𝑧)
𝛺

𝑑𝑧𝑑𝑥
𝛺

 

= ∫ 1𝐵(𝑥,𝑥)𝜋(𝑥)𝑑𝑥𝛺 . 

Consequently, the chain assumes 𝜋 as a stationary distribution. 

A measure 𝜋(∙) is said to be irreducible for the transition kernel 𝑃𝑀𝐻 of a Markov chain if:  

∀A ∈ ℬ so that  𝜋(A) > 0  ⇒ ∃𝑡 :  𝑃𝑀𝐻
𝑡 (𝑥, A) > 0 

Let A be a borelean of ℬ, and for t=1 we have: 

∫ 1𝐵(𝑥,A)𝑃𝑀𝐻(𝑥, A)𝑑𝑥
Ω

= ∫ 1𝐵(𝑥,A)𝑎(𝑥, A)𝑞(𝑥, A)𝑑𝑥
𝛺

+ ∫ 1𝐵(𝑥,A)
𝛺

[∫ 1 − 𝑎(𝑥, 𝑧)𝑞(𝑥, 𝑧)𝑑𝑧
𝛺

] 𝛿𝑥 (𝐴)𝑑𝑥 

                            = ∫ 1𝐵(𝑥,A)𝑎(𝑥, A)𝑞(𝑥, A)dx
Ω

+ ∫ 1B(𝑥,𝑥)
Ω

[∫ 1 − 𝑎(𝑥, 𝑧)𝑞(𝑥, 𝑧)𝑑𝑧
Ω

] 𝑑𝑥 

             = ∫ 1𝐵(𝑥,𝐴)𝑎(𝑥, A)𝑞(𝑥, A)𝑑𝑥
Ω

+ 1 − ∫ ∫ 𝑎(𝑥, 𝑧)𝑞(𝑥, 𝑧)𝑑𝑧
𝛺𝛺

𝑑𝑥. 

Since 𝑎(𝑥, A) = 𝑚𝑖𝑛 (1; 𝛾𝑐(𝑥)−𝑐(A)) and 𝑎(𝑥, 𝑧) = 𝑚𝑖𝑛 (1; 𝛾𝑐(𝑥)−𝑐(𝑧)), then four possible cases 

can be distinguished: 

• If 𝑎(𝑥, A) = 1 and 𝑎(𝑥, 𝑧) = 1  then: 

∫ 1𝐵(𝑥,A)𝑃𝑀𝐻 (𝑥, A)𝑑𝑥
𝛺

= ∫ 1𝐵(𝑥,A)𝑞(𝑥, A)𝑑𝑥
𝛺

+ 1 − ∫ ∫ 𝑞(𝑥, 𝑧)𝑑𝑧
𝛺𝛺

𝑑𝑥 

           = ∫ 1𝐵(𝑥,A)𝑞(𝑥, A)𝑑𝑥
𝛺

> 0. 

• If 𝑎(𝑥, A) = 1 and 𝑎(𝑥, 𝑧) = 𝛾𝑐(𝑥)−𝑐(𝑧)  then:  

∫ 1𝐵(𝑥,A)𝑃𝑀𝐻 (𝑥, A)𝑑𝑥
𝛺

= ∫ 1𝐵(𝑥,A)𝑞(𝑥, A)𝑑𝑥
𝛺

+ 1 − ∫ ∫ 𝛾𝑐(𝑥)−𝑐(𝑧)𝑞(𝑥, 𝑧)𝑑𝑧
𝛺𝛺

𝑑𝑥 

         = ∫ 1𝐵(𝑥,A)𝑞(𝑥, A)𝑑𝑥
𝛺

+ 1 − 𝛾𝑐(𝑥)−𝑐(𝑧) > 0. 

• If 𝑎(𝑥, A) = 𝛾𝑐(𝑥)−𝑐(A)  and  𝑎(𝑥, 𝑧) = 1 then: 



11 Int. J. Anal. Appl. (2023), 21:51 

 

∫ 1𝐵(𝑥,A)𝑃𝑀𝐻(𝑥, A)𝑑𝑥
𝛺

= ∫ 1𝐵(𝑥,A)𝛾
𝑐(𝑥)−𝑐(A)𝑞(𝑥, A)𝑑𝑥

𝛺

+ 1 − ∫ ∫ 𝑞(𝑥, 𝑧)𝑑𝑧
𝛺𝛺

𝑑𝑥 

     = 𝛾𝑐(𝑥)−𝑐(A) ∫ 1𝐵(𝑥,A)𝑞(𝑥, A)𝑑𝑥 > 0.
𝛺

 

• If 𝑎(𝑥, A) = 𝛾𝑐(𝑥)−𝑐(𝐴) and 𝑎(𝑥, 𝑧) = 𝛾𝑐(𝑥)−𝑐(𝑧) then : 

∫ 1𝐵(𝑥,A)𝑃𝑀𝐻(𝑥, A)𝑑𝑥
𝛺

= ∫ 1𝐵(𝑥,A)𝛾
𝑐(𝑥)−𝑐(A)𝑞(𝑥, A)𝑑𝑥

𝛺

+ 1 − ∫ ∫ 𝛾𝑐(𝑥)−𝑐(𝑧)𝑞(𝑥, 𝑧)𝑑𝑧
𝛺𝛺

𝑑𝑥 

                                        = 𝛾𝑐(𝑥)−𝑐(A) ∫ 1𝐵(𝑥,𝐴)𝑞(𝑥, A)𝑑𝑥
𝛺

+ 1 − 𝛾𝑐(𝑥)−𝑐(𝑧) ∫ ∫ 𝑞(𝑥, 𝑧)𝑑𝑧
𝛺𝛺

𝑑𝑥 

                                                       = 𝛾𝑐(𝑥)−𝑐(A) ∫ 1𝐵(𝑥,𝐴)𝑞(𝑥, A)𝑑𝑥𝛺 + 1 − 𝛾
𝑐(𝑥)−𝑐(𝑧) > 0. 

So ∫ 1𝐵(𝑥,A)𝑃𝑀𝐻
𝑡 (𝑥, A)𝑑𝑥

𝛺
> 0  ∀ 𝑡 ≥ 0, then 𝑃𝑀𝐻 is 𝜋 −irreducible. 

Since 𝜋  is the invariant distribution of 𝑃𝑀𝐻, then it remains so with 𝑃. As a matter of fact, 

𝜋𝑃𝑀𝐻 = 𝜋, and by recurrence on 𝑛 we get: 

𝜋𝑃𝑀𝐻 = 𝜋𝑃𝑀𝐻
2 = 𝜋𝑃𝑀𝐻

3 = ⋯ = 𝜋𝑃𝑀𝐻
𝑛 = 𝜋 

Yet  𝑃 = 𝑃𝑀𝐻
𝑛 , then we get:  𝜋𝑃 = 𝜋. On the other hand, 𝜋 − reversibility of 𝑃𝑀𝐻 leads to 𝜋 − 

reversibility of 𝑃, i.e: 

𝜋(𝑥)𝑃𝑀𝐻 (𝑥, 𝑦)  =  𝜋(𝑦)𝑃𝑀𝐻 (𝑦, 𝑥) ⇒ 𝜋(𝑥)𝑃(𝑥, 𝑦)  =  𝜋(𝑦)𝑃(𝑦, 𝑥) 

Since  𝜋𝑃𝑀𝐻 = 𝜋𝑃𝑀𝐻
𝑛 = 𝜋, on then we shall get: 

𝜋(𝑥)𝑃𝑀𝐻 (𝑥, 𝑦) = 𝜋(𝑥)𝑃𝑀𝐻
𝑛(𝑥, 𝑦) = 𝜋(𝑦)𝑃𝑀𝐻 (𝑦, 𝑥) = 𝜋(𝑦)𝑃𝑀𝐻

𝑛(𝑦, 𝑥) 

And since  𝑃𝑀𝐻
𝑛 = 𝑃, then we get: 𝜋(𝑥)𝑃(𝑥, 𝑦)  =  𝜋(𝑦)𝑃(𝑦, 𝑥). 

By constructing 𝑃 = 𝑃𝑀𝐻
𝑛 , the  𝜋 −irreducibility of 𝑃𝑀𝐻 leads to 𝜋 −irreducibility of 𝑃. If 𝑃 

is 𝜋 −irreducibile and has an invariant 𝜋 distribution, then 𝑃 is positive recurrent and π is the unique 

invariant distribution of  𝑃 [13] (see, Proposition 1). 

On the other hand, the chain generated by the construction algorithm shall be aperiodic as 

well provided there exists at least a pair of configuration (𝑥, 𝑦) so that 𝑎(𝑥, 𝑦) < 1, and we shall 

ultimately get  𝑃(𝑥, 𝑥)  >  0. We are quick to notice that the chain is aperiodic since the event 

𝑋(𝑁+1)  =  𝑋(𝑁) is likely at practically any moment. As a matter of fact, each state can, then, be 

visited at two subsequent iterations, thus  𝑃1(𝑥, 𝑥) >  0, and their period is then equal to 1. 

Since the chain generated with the algorithm is irreducible and aperiodic, then its kernel with 

transition 𝑃 is primitive (a characterization of primitive Markov Kernel more common in probability 

theory is to say that are irreducible and aperiodic [14]). 



12 Int. J. Anal. Appl. (2023), 21:51 

 

Theorem 4.1. The Markov chain (𝑋𝑘 )𝑘 ≥ 0 obtained from the construction algorithm is uniformly 

ergodic, and its kernel 𝑃 realizes the simulation of the connected component process with 

density 𝜋(𝑥) = 𝑘𝛽𝑛(𝑥)𝛾−𝑐(𝑥), i.e.  𝑣𝑃𝑚   converges to 𝜋 as 𝑚 approaches infinity; where 𝑣 is an initial 

distribution, and we have : 

‖𝑣𝑃𝑚 − 𝜋‖ → 0, 𝑚 → ∞ 

Proof of theorem 4.1 Let 𝑣 be an initial distribution, for every integer 𝑚 and ∀ 𝑥 ∈  𝑁𝑓  we have: 

‖𝑣𝑃𝑚 (𝑥,∙) − 𝜋‖ = ‖𝑣𝑃𝑚 − 𝜋𝑃𝑚 ‖ 

On the other hand,  ‖𝑣𝑃𝑚 − 𝜋𝑃𝑚 ‖ ≤ 2𝐶(𝑃𝑚 ) and 𝐶(𝑃𝑚 ) ≤ (𝐶(𝑃))𝑚 [15] (see lemma 4.2.2, 

P.71), this implies that:  

‖𝑣𝑃𝑚 (𝑥,∙) − 𝜋‖ ≤ 2(𝐶(𝑃))𝑚 

Where, 𝐶(𝑃)is the Dobrushin contraction coefficient of  𝑃 [16]. According to proposition 4.1, 

the kernel  𝑃 is primitive, then 𝐶(𝑃) < 1 [15] (see lemma 4.2.3, P.72), so when 𝑚 approaches infinity, 

‖𝑣𝑃𝑚 − 𝜋‖ tends towards zero. Hence, the chain is uniformly ergodic and converges towards the 

distribution defined in (3.2). 

 

5. SIMULATION STUDY AND DISCUSSION 

In this section, we will compare the quality of point distributions in the proposed designs 

using standard criteria to evaluate the degree of filling of the experimental space and the quality of 

the uniform distribution. 

 

• Criterion of distance (Mindist): it is about maximizing the minimal distance between two 

points of the design [17].  

𝑀𝑖𝑛𝑑𝑖𝑠𝑡 = min
𝑖

min
𝑗≠𝑖

𝑑(𝑥𝑖 , 𝑥𝑗 ) 

Where, 𝑑(𝑥𝑖 , 𝑥𝑗 ) is the Euclidean distance between the points 𝑥𝑖 and 𝑥𝑗. 

• The coverage criterion (Cov): it helps to measure the gap between the points of the 

design and those of a regular grid. This criterion is invalid for a regular grid. Our goal, 

here, is to minimize this gap to get closer to a regular grid, hence ensuring the filling of 

the space without, however, reaching it to comply with a uniform distribution, mainly in 

projection on factorial axes [18]. 



13 Int. J. Anal. Appl. (2023), 21:51 

 

𝑐𝑜𝑣 =
1

𝛿̅
√

1

𝑛
∑(𝛿𝑖 − 𝛿̅)

2
𝑛

𝑖=1

 

Where 𝛿𝑖 = min
𝑖≠𝑗

(𝑑(𝑥𝑖 , 𝑥𝑗 )) and  �̅� =
1

𝑛
 ∑ 𝛿𝑖

𝑛
𝑖=1  . For regular grid, 𝛿1 = 𝛿2 = ⋯ = 𝛿𝑛, then Cov=0. 

• The mesh Ratio (R): is the ratio between maximal and minimal distances separating the 

points of the experimental design. In the case of a regular grid, R=1. Thus, when R is 

close to the value 1 the points are close to those of a regular grid.  

𝑅 =
max
𝑖:1…𝑛

𝛿𝑖

min
𝑖:1…𝑛

𝛿𝑖
 

• The Discrepancy criterion (Disc): Discrepancy measures the difference between the empirical 

distribution function of the design points and that of the uniform distribution. Unlike the 

previous three criteria, Discrepancy does not depend on the distance between points. There 

are different ways to measure Discrepancy, but in this study, we use the L2 norm. [19]. 

𝐷𝑖𝑠𝑐 = (
1

3
)

𝑝

−
21−𝑝

𝑛
∑ ∏ (1 − (𝑥𝑖

𝑗
)

2
) +

1

𝑛2
∑ ∑ ∏ (1 − max(𝑥𝑖

𝑗
− 𝑥𝑘

𝑗
))

𝑝

𝑗=1

𝑛

𝑘=1

𝑛

𝑖=1

𝑝

𝑗=1

𝑛

𝑖=1

 

In table 1, we present the results related to the discrepancy criterion for the proposed designs 

and compare them with sequences with low discrepancy (such as Halton sequence [20], Sobol 

sequence [21], and Faure sequence [22]). It is interesting to note that the proposed designs have 

similarly low discrepancies as those of the low discrepancy sequences. 

Table 1. Discrepancy values for CCD designs: a low discrepancy Halton sequence, a low 
discrepancy Sobol sequence and a low discrepancy Faure sequence. The presented results cover 

three, five, seven, and ten dimensions. 

The underlined values represent the lowest values for each dimension. 

Dimension 
Number of 

Points 
CCD 

Halton 

Sequence 

Sobol 

Sequence 

Faure 

Sequence 

3 

5 

7 

10 

20 

40 

50 

60 

0,001224 

0,000634 

0,000257 

0.0000140 

0,00173285 

0,00040004 

0,00012451 

0,00005601 

0,001382617 

0,000260329 

0,000074156 

0,000007428 

0,00113882 

0,00028458 

0,00012455 

0,00001327 



14 Int. J. Anal. Appl. (2023), 21:51 

 

The designs presented in this article are also compared with commonly used designs in 

computer experiments, excluding sequences with weak discrepancy. To ensure meaningful results, 

the criteria have been tested on 80 designs. The designs considered in this section are typically the 

following: 

• Random Designs (RD). 

• Latin Hypercube Sampling (LHS) [23]. 

• Maximin LHS Designs (mLHS) [24]. 

• Strauss Designs (SD) [6]. 

• Maximal Entropy Designs (Dmax) [25]. 

• Marked Strauss Designs (MSD) ([7], [8]). 

• Connected Component Designs (CCD). 

Figures 4 and 5 present the criteria results in box plots for 80 designs in 10 and 5 dimensions, 

respectively. The plots clearly illustrate the distribution of the data and allow for easy comparison 

between designs. 

  

  

Figure 4. Representations of usual criteria calculated on 80 designs with 60 points for dimension 
10. 



15 Int. J. Anal. Appl. (2023), 21:51 

 

  

  

Figure 5. Representations of usual criteria calculated on 80 designs with 40 points for dimension 
5. 

Some remarks regarding the figures above are worth noting. Maximal entropy designs, Strauss 

designs, Marked Strauss designs, as well as connected component designs all scored well with regard 

to the discrepancy criterion. It is interesting to note that connected component designs are among 

the designs previously mentioned which also have very good distance criteria. The application of the 

recovery criterion helps to show that the proposed designs offer better results regarding this criterion. 

 

6. CONCLUSION 

The use of the connected component process and Markov Chain Monte Carlo (MCMC) 

method is instrumental in creating novel computer designs that are customized based on the 

distribution of a connected component model. This method is highly versatile, as it allows us to 

manipulate the distribution by representing it in a manner that makes it possible to impose certain 

properties such as filling. Furthermore, MCMC methodology presents a fascinating alternative to the 

classical statistical approach that typically involves the independent realization of the same 



16 Int. J. Anal. Appl. (2023), 21:51 

 

distribution. Ultimately, the designs developed in this study were compared to those commonly used 

in computer experiments, and the outcomes were highly satisfactory. 

Conflicts of Interest: The authors declare that there are no conflicts of interest regarding the 

publication of this paper. 

References  

[1] D.J. Daley, D. Vere-Jones, An Introduction to the Theory of Point Processes. Volume I: Elementary 

Theory and Methods, 2nd edn, Springer, Berlin, 1988. 

[2] A.J. Baddeley, J.  Møller, Nearest Neighbour Markov Point Processes and Random Sets. Int. Stat. Rev. 

57 (1989), 89–121. https://doi.org/10.2307/1403381. 

[3] A.J. Baddeley, M.N.M. van Lieshout, Area-Interaction Point Processes, Ann. Inst. Stat. Math. 47 (1995), 

601–619. https://doi.org/10.1007/bf01856536. 

[4] Y.C. Chin, A.J. Baddeley, On Connected Component Markov Point Processes, Adv. Appl. Probab. 31 

(1999), 279–282. https://doi.org/10.1239/aap/1029955135. 

[5] B.D. Ripley, F.P. Kelly, Markov Point Processes, J. Lond. Math. Soc. s2-15 (1977), 188–192.  

https://doi.org/10.1112/jlms/s2-15.1.188. 

[6] J. Franco, Planification d’Expériences Numériques en Phase Exploratoire pour des Codes de Calculs 

Simulant des Phénomènes Complexes. Doctoral Thesis, l’Ecole Nationale Supérieure des Mines de 

SaintEtienne, France, 2008. 

[7] H. Elmossaoui, N. Oukid, F. Hannane, Construction of Computer Experiment Designs Using Marked Point 

Processes, Afr. Mat. 31 (2020), 917–928. https://doi.org/10.1007/s13370-020-00770-9. 

[8] H. Elmossaoui, Contribution à la Méthodologie de la Recherche Expérimentale, Doctoral Thesis, University 

Saad Dahleb, Blida, Algeria, 2020. 

[9] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, E. Teller, Equation of State Calculations 

by Fast Computing Machines, J. Chem. Phys. 21 (1953), 1087–1092. 

https://doi.org/10.1063/1.1699114. 

[10] W.K. Hastings, Monte Carlo Sampling Methods Using Markov Chains and Their Applications, Biometrika. 

57 (1970), 97–109. https://doi.org/10.1093/biomet/57.1.97. 

[11] A. Baddeley, J. Møller, A.G. Pakes, Properties of Residuals For Spatial Point Processes, Ann. Inst. Stat. 

Math. 60 (2007), 627–649. https://doi.org/10.1007/s10463-007-0116-6. 

[12] S. Chib, E. Greenberg, Understanding the Metropolis-Hastings Algorithm, Amer. Stat. 49 (1995), 327-

355. https://doi.org/10.2307/2684568. 

https://doi.org/10.2307/1403381
https://doi.org/10.1007/bf01856536
https://doi.org/10.1239/aap/1029955135
https://doi.org/10.1112/jlms/s2-15.1.188
https://doi.org/10.1007/s13370-020-00770-9
https://doi.org/10.1063/1.1699114
https://doi.org/10.1093/biomet/57.1.97
https://doi.org/10.1007/s10463-007-0116-6
https://doi.org/10.2307/2684568


17 Int. J. Anal. Appl. (2023), 21:51 

 

[13] S. Chib, E. Greenberg, Markov Chain Monte Carlo Simulation Methods in Econometrics, Econ. Theory. 

12 (1996), 409–431. https://doi.org/10.1017/s0266466600006794. 

[14] E. Senata, Non-Negative Matrices and Markov Chains, 2nd edition. Springer, New York Heidelberg Berlin, 

1981. 

[15] G. Winkler, Image Analysis Random Fields and Dynamic Monte Carlo Methods, Springer, Berlin, 1995. 

[16] R.L. Dobrushin, Central Limit Theorem for Nonstationary Markov Chains. I, Theory Probab. Appl. 1 

(1956), 65–80. https://doi.org/10.1137/1101006. 

[17] M. Gunzburger, J. Burkardt, Uniformity Measures for Point Samples in Hypercubes. (2004). 

 https://people.sc.fsu.edu/~jburkardt/publications/gb_2004.pdf. 

[18] M.E. Johnson, L.M. Moore, D. Ylvisaker, Minimax and Maximin Distance Designs, J. Stat. Plan. 

Inference. 26 (1990), 131–148. https://doi.org/10.1016/0378-3758(90)90122-b. 

[19] T.T. Warnock, Computational Investigations of Low-Discrepancy Point Sets II. In: Niederreiter H., and 

Shiue P.JS. (eds) Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing. Lecture Notes 

in Statistics, 106. Springer, New York, 1995. 

[20] J.H. Halton, On the Efficiency of Certain Quasi-Random Sequences of Points in Evaluating Multi-

Dimensional Integrals, Numer. Math. 2 (1960), 84–90. https://doi.org/10.1007/bf01386213. 

[21] I.M. Sobol, Uniformly Distributed Sequences with an Additional Uniform Property, USSR Comput. Math. 

Math. Phys. 16 (1976), 236–242. https://doi.org/10.1016/0041-5553(76)90154-3. 

[22] H. Faure, Discrépance de Suites Associées à un Système de Numération (en Dimension s), Acta Arith. 

41 (1982), 337–351. https://doi.org/10.4064/aa-41-4-337-351. 

[23] W.L. Loh, On Latin Hypercube Sampling, Ann. Stat. 24 (1996), 2058-2080. 

 https://doi.org/10.1214/aos/1069362310. 

[24] M.D. Morris, T.J. Mitchell, Exploratory Designs for Computational Experiments, J. Stat. Plan. Inference. 

43 (1995), 381–402. https://doi.org/10.1016/0378-3758(94)00035-t. 

[25] M.C. Shewry, H.P. Wynn, Maximum Entropy Sampling, J. Appl. Stat. 14 (1987), 165–170. 

 https://doi.org/10.1080/02664768700000020. 

https://doi.org/10.1017/s0266466600006794
https://doi.org/10.1137/1101006
https://people.sc.fsu.edu/~jburkardt/publications/gb_2004.pdf
https://doi.org/10.1016/0378-3758(90)90122-b
https://doi.org/10.1007/bf01386213
https://doi.org/10.1016/0041-5553(76)90154-3
https://doi.org/10.4064/aa-41-4-337-351
https://doi.org/10.1214/aos/1069362310
https://doi.org/10.1016/0378-3758(94)00035-t
https://doi.org/10.1080/02664768700000020