generating random x and y coordinates with a minimum distance

Question

Is there a way in R to generate random coordinates with a minimum distance between them?

E.g. what I'd like to avoid

x <- c(0,3.9,4.1,8)
y <- c(1,4.1,3.9,7)
plot(x~y)

Yodel

How do you define minimal distance? – Dodecahedron

minimum euclidean distance? For example, points are generated with at least a distance of 1 between them – Yodel

Answer 1

This is a classical problem from stochastic geometry. Completely random points in space where the number of points falling in disjoint regions are independent of each other corresponds to a homogeneous Poisson point process (in this case in R^2, but could be in almost any space).

An important feature is that the total number of points has to be random before you can have independence of the counts of points in disjoint regions.

For the Poisson process points can be arbitrarily close together. If you define a process by sampling the Poisson process until you don't have any points that are too close together you have the so-called Gibbs Hardcore process. This has been studied a lot in the literature and there are different ways to simulate it. The R package spatstat has functions to do this. rHardcore is a perfect sampler, but if you want a high intensity of points and a big hard core distance it may not terminate in finite time... The distribution can be obtained as the limit of a Markov chain and rmh.default lets you run a Markov chain with a given Gibbs model as its invariant distribution. This finishes in finite time but only gives a realisation of an approximate distribution.

In rmh.default you can also simulate conditional on a fixed number of points. Note that when you sample in a finite box there is of course an upper limit to how many points you can fit with a given hard core radius, and the closer you are to this limit the more problematic it becomes to sample correctly from the distribution.

Example:

library(spatstat)
beta <- 100; R = 0.1
win <- square(1) # Unit square for simulation
X1 <- rHardcore(beta, R, W = win) # Exact sampling -- beware it may run forever for some par.!
plot(X1, main = paste("Exact sim. of hardcore model; beta =", beta, "and R =", R))

minnndist(X1) # Observed min. nearest neighbour dist.
#> [1] 0.102402

Approximate simulation

model <- rmhmodel(cif="hardcore", par = list(beta=beta, hc=R), w = win)
X2 <- rmh(model)
#> Checking arguments..determining simulation windows...Starting simulation.
#> Initial state...Ready to simulate. Generating proposal points...Running Metropolis-Hastings.
plot(X2, main = paste("Approx. sim. of hardcore model; beta =", beta, "and R =", R))

minnndist(X2) # Observed min. nearest neighbour dist.
#> [1] 0.1005433

Approximate simulation conditional on number of points

X3 <- rmh(model, control = rmhcontrol(p=1), start = list(n.start = 42))
#> Checking arguments..determining simulation windows...Starting simulation.
#> Initial state...Ready to simulate. Generating proposal points...Running Metropolis-Hastings.
plot(X3, main = paste("Approx. sim. given n =", 42))

minnndist(X3) # Observed min. nearest neighbour dist.
#> [1] 0.1018068

Savoirvivre

Thank you for this great answer! Do you know a way that I can achieve a 1D version of this? i.e random points in 1-dimensional space where each point is at least r units away from an adjacent – Kerseymere

A simple approach would be to generate the distance between points sequentially as r+X where X is exponentially distributed. It probably requires a new question to give more details. – Savoirvivre

Answer 2

OK, how about this? You just generate random number pairs without restriction and then remove the onces which are too close. This could be a great start for that:

minimumDistancePairs <- function(x, y, minDistance){
    i <- 1
    repeat{
        distance <- sqrt((x-x[i])^2 + (y-y[i])^2) < minDistance # pythagorean theorem
        distance[i] <- FALSE # distance to oneself is always zero
        if(any(distance)) { # if too close to any other point
            x <- x[-i] # remove element from x
            y <- y[-i] # and remove element from y
        } else { # otherwise...
            i = i + 1 # repeat the procedure with the next element
        }
        if (i > length(x)) break
    }
    data.frame(x,y)
}

minimumDistancePairs(
    c(0,3.9,4.1,8)
    , c(1,4.1,3.9,7)
    , 1
)

will lead to

x   y
1 0.0 1.0
2 4.1 3.9
3 8.0 7.0

Be aware, though, of the fact that these are not random numbers anymore (however you solve problem).

Sudra

Answer 3

You can use rejection sapling https://en.wikipedia.org/wiki/Rejection_sampling

The principle is simple: you resample until you data verify the condition.

> set.seed(1)
> 
> x <- rnorm(2)
> y <- rnorm(2)
> (x[1]-x[2])^2+(y[1]-y[2])^2
[1] 6.565578
> while((x[1]-x[2])^2+(y[1]-y[2])^2 > 1) {
+   x <- rnorm(2)
+   y <- rnorm(2)
+ }
> (x[1]-x[2])^2+(y[1]-y[2])^2
[1] 0.9733252
> 

Unific

Thanks. I think that this is the best method I've seen suggested – Yodel

I don't think so. This method only works on 2 points. It does not check if a point has sufficient distance to the previously generated points. Either I misunderstood the question or this is not the answer. – Sudra

Answer 4

The following is a naive hit-and-miss approach which for some choices of parameters (which were left unspecified in the question) works well. If performance becomes an issue, you could experiment with the package gpuR which has a GPU-accelerated distance matrix calculation.

rand.separated <- function(n,x0,x1,y0,y1,d,trials = 1000){
  for(i in 1:trials){
    nums <- cbind(runif(n,x0,x1),runif(n,y0,y1))
    if(min(dist(nums)) >= d) return(nums)
  }
  return(NA) #no luck
}

This repeatedly draws samples of size n in [x0,x1]x[y0,y1] and then throws the sample away if it doesn't satisfy. As a safety, trials guards against an infinite loop. If solutions are hard to find or n is large you might need to increase or decrease trials.

For example:

> set.seed(2018)
> nums <- rand.separated(25,0,10,0,10,0.2)
> plot(nums)

runs almost instantly and produces:

Tactual

one change to this: instead of seq_along(trials) you need 1:trials. At the moment it only does one replicate – Yodel

@DomBurns embarrassing mistake -- thanks for the correction – Tactual

@JohnColeman, This seems to be a solution that I am interested (in this So Q here)[#65583993. Is it possible to set the co-ordinates 2 two groups , so that they fall left and right ? Any leads would be helpful. – Woolfell

Answer 5

Im not sure what you are asking.

if you want random coordinates here.

c(
runif(1,max=y[1],min=x[1]),
runif(1,max=y[2],min=x[2]),
runif(1,min=y[3],max=x[3]),
runif(1,min=y[4],max=x[4])
)

Epp