straussDesign {DiceDesign}

Designs based on Strauss process

Description

Space-Filling Designs based on Strauss process

Usage

straussDesign(n,dimension, RND, alpha=0.5, repulsion=0.001, NMC=1000,
              constraints1D=0, repulsion1D=0.0001, seed=NULL)

Arguments

n

the number of experiments

dimension

the number of input variables

RND

a real number which represents the radius of interaction

alpha

the potential power (default, fixed at 0.5)

repulsion

the repulsion parameter in the unit cube (gamma)

NMC

the number of McMC iterations (this number must be large to converge)

constraints1D

1 to impose 1D projection constraints, 0 otherwise

repulsion1D

the repulsion parameter in 1D

seed

seed for the uniform generation of number

Details

Strauss designs are Space-Filling designs initially defined from Strauss process:

\pi (X) = k \gamma^{s(X)}

where s(X) is is the number of pairs of points (x^{i}, x^{j}) of the design X = \left( x^{1}, \ldots, x^{n} \right) that are separated by a distance no greater than the radius of interaction RND, k is the normalizing constant and \gamma is the repulsion parameter. This distribution corresponds to the particular case alpha=0.

For the general case, a stochastic simulation is used to construct a Markov chain which converges to a spatial density of points \pi(X) described by the Strauss-Gibbs potential. In practice, the Metropolis-Hastings algorithm is implemented to simulate a distribution of points which converges to the stationary law:

\pi (X) \propto exp(-U(X))

with a potentiel U defined by:

U(X) = \beta \sum_{1 \leq i < j \leq n} \varphi \left( \| x^{i}-x^{j} \| \right)

where \beta = - \ln \gamma, \varphi (h) = \left( 1 - \frac{h}{RND} \right) ^{\alpha} if h \leq RND and 0 otherwise.

The input parameters of straussDesign function can be interpreted as follows:

- RND is used to compute the number of pairs of points of the design separated by a distance no more than RND. A point is said "in interaction" with another if the spheres of radius RND/2 centered on these points intersect.

- alpha is the potential power \alpha. The case alpha=0 corresponds to Strauss process (0-1 potential).

- repulsion is equal to the \gamma parameter of the Strauss process. Note that \gamma belongs to ]0,1].

- constraints1D allows to specify some constraints into the margin. If constraints1D==1, two repulsion parameters are needed: one for the all space (repulsion) and the other for the 1D projection (repulsion1D). Default values are repulsion=0.001 and repulsion1D=0.001. Note that the value of the radius of interaction in the one-dimensional axis is not an input parameter and is automatically fixed at 0.75/n.

Value

A list containing:

n

the number of experiments

dimension

the number d of variables

design_init

the initial distribution of n points [0,1]^{d}

radius

the radius of interaction

alpha

the potential power alpha

repulsion

the repulsion parameter \gamma

NMC

the number of iterations McMC

constraints1D

an integer indicating if constraints on the factorial axis are imposed. If its value is different from zero, a component repulsion1D containing the value of the repulsion parameter \gamma in dimension 1 is added at the list.

design

the design of experiments in [0,1]^{d}

seed

the seed corresponding to the design

Author(s)

J. Franco

References

J. Franco, X. Bay, B. Corre and D. Dupuy (2008) Planification d'experiences numeriques a partir du processus ponctuel de Strauss, https://hal.science/hal-00260701/fr/.

Examples

## Strauss-Gibbs designs in dimension 2 (n=20 points)
S1 <- straussDesign(n=20, dimension=2, RND=0.2)
plot(S1$design, xlim=c(0,1), ylim=c(0,1))
theta <- seq(0,2*pi, by=2*pi/(100 - 1))
for(i in 1:S1$n){
   lines(S1$design[i,1]+S1$radius/2*cos(theta),
	     S1$design[i,2]+S1$radius/2*sin(theta), col='red')
}

## 2D-Strauss design
S2 <- straussDesign(n=20, dimension=2, RND=0.2, NMC=200,
	         constraints1D=0, alpha=0, repulsion=0.01)
plot(S2$design,xlim=c(0,1),ylim=c(0,1))

## 2D-Strauss designs with constraints on the axis
S3 <- straussDesign(n=20, dimension=2, RND=0.18, NMC=200, 
	         constraints1D=1, alpha=0.5, repulsion=0.1, repulsion1D=0.01)
plot(S3$design, xlim=c(0,1),ylim=c(0,1))
rug(S3$design[,1], side=1)
rug(S3$design[,2], side=2)

## Change the dimnames, adjust to range (-10, 10) and round to 2 digits
xDRDN(S3, letter="T", dgts=2, range=c(-10, 10))