Suppose there is a square box with side length $m$ (measured in pixels). Let there be $n$ points in this box, distributed uniformly within the box (with integer coordinates, aligned to a pixel grid). If we take from each point the Euclidean distance to its nearest neighbor, what would be the expected value of this distance averaged over all points?

My actual problem is about a discrete pixel grid, an $m\times m$ bitmap image, but if that's easier I would be happy with a continuous solution. A more general solution e.g. for a rectangle instead of a square box is welcome, but at this point it is not necessary. I found similar questions about the continous case, but without answers. For me it wouldn't be easy to generalise the case for two points only.

  • 1
    $\begingroup$ A simulation is your best option. $\endgroup$ – quasi Aug 31 '17 at 23:12
  • $\begingroup$ When you say "distributed randomly"... Do you mean uniformly? There are many ways to put points in a box at random. $\endgroup$ – Merkh Aug 31 '17 at 23:29
  • $\begingroup$ Sorry, yes, I did mean uniformly: each possible coordinate pair has an equal chance. $\endgroup$ – Adam Aug 31 '17 at 23:33
  • $\begingroup$ This similar problem was mentioned in an answer earlier: mathoverflow.net/questions/124579/… Unfortunately this answer has since been deleted, probably because stricly speaking it is not an answer to the question. I think it is related and interesting, so worth mentioning anyway. $\endgroup$ – Adam Sep 1 '17 at 12:40
  • $\begingroup$ @Merk Usually, it means uniformly. It's the colloquial way. $\endgroup$ – Felix Marin Jan 8 at 3:52

Empirically it is between $\dfrac{0.5}{\sqrt{n-1}}$ and $\dfrac{0.56}{\sqrt{n-1}}$, tending towards the lower bound as $n$ increases. This is not very different to putting $n-1$ random points in a circle or square with area $1$ and finding the expected minimum distance to the centre.

At the bottom end with $n=2$, the question becomes the expected Euclidean distance between two random points in a square, which an earlier answer gives as $\frac{1}{15}(2+\sqrt{2}+5\log(1+\sqrt{2})) \approx 0.5214$

Simulation results look like this, both directly and multiplying by $\sqrt{n-1}$

enter image description here

with values of

n    expected dist nn     sqrt(n-1) times expected dist nn
2           0.521                   0.521
3           0.389                   0.550
4           0.322                   0.557
5           0.279                   0.558
6           0.250                   0.558
7           0.227                   0.557
8           0.210                   0.555
9           0.196                   0.553
10          0.184                   0.552
20          0.124                   0.540
50          0.075                   0.527
100         0.052                   0.519
200         0.036                   0.514
500         0.023                   0.509
1000        0.016                   0.507

(See comments) The simulation code was

set.seed(1)          # to reproduce
maxn <- 1024         # top number for n-1 neighbours 
cases <- 10^6        # simulations 
x0 <- runif(cases)   # base points  
y0 <- runif(cases)
mindistsq <- rep(Inf, cases)
emindist <- numeric(maxn)
for (i in 1:maxn){
  x <- runif(cases)  # possible neighbour
  y <- runif(cases)
  distsq <- (x0-x)^2 + (y0-y)^2
  mindistsq <- ifelse(distsq < mindistsq, distsq, mindistsq)
  emindist[i] <- mean(sqrt(mindistsq))
emindistsqrtn1 <- emindist * sqrt(1:maxn)   # multiplying by sqrt(n-1)
plot(2:(maxn+1), emindist, xlab="number of points in square", 
                     main="Expected distance of nearest neighbour in square",
                     ylab="Expected dist nearest neighbour",
               ylim = c(0,0.6))
plot(2:(maxn+1), emindistsqrtn1 , xlab="number of points in square", 
               ylab="sqrt(n-1) times expct dist nn",
               main="sqrt(n-1) times Expected distance of nearest neighbour",
               ylim = c(0.5,0.56))
exn1 <- c(1:9, 19, 49, 99, 199, 499, 999)
tab <- cbind(round(emindist[exn1],3), round(emindistsqrtn1[exn1],3))  
colnames(tab) <- c("expected min","sqrt(n-1) times exp min")
rownames(tab) <- (exn1+1) 
  • $\begingroup$ So I guess these results are for a unit square ($m=1$). Did you use floating point numbers as the coordinates of the points? I suppose you used a pseudorandom number generator to get these coordinates. How many times did you calculate these values for each $n$? $\endgroup$ – Adam Sep 1 '17 at 12:59
  • $\begingroup$ @Adam. As you guessed, this is based on a unit square but you can just scale distances by the side of the square. I used R, which has double precision (so nearly 16 decimal digits), and I simulated $1$ million points with their nearest neighbours for each $n$. So I expect the rounded numbers in my table are accurate to $\pm0.001$ or better. I will add the simulation code $\endgroup$ – Henry Sep 1 '17 at 13:33

For the continuous version, assume that you pick $n$ points randomly from $[0, 1]^2$. You are trying to find $$E\left[ \frac{1}{n}\sum_{k=1}^n \min_{i, 1 \le i \le n, i \not=k}||X_k-X_i|| \right] \tag 1$$

where each of $X_i$ are randomly and uniformly chosen from $[0, 1]^2$.

Because the expected value of a sum of random variables is equal to the sum of the expected values, this simplifies to $$ \frac{1}{n}\sum_{k=1}^n E\left[\min_{i, 1 \le i \le n, i \not=k}||X_k-X_i|| \right] \tag 2$$

By symmetry, this is $$E\left[\min_{i, 1 \le i \le n-1}||X_n-X_i|| \right] \tag 3$$ Let $F_n(x) = P\left(\min_{i, 1 \le i \le n-1}||X_n-X_i|| > x\right)$, which would be $1$ for $x < 0$ and $0$ for $x > \sqrt{2}$. Then $(3)$ is $\int_0^{\sqrt{2}}F_n(x) dx$. $F_n(x)$ is equal to $$P\left(||X_n-X_1|| > x\right)^{n-1} \tag 4$$

The probability $P\left(||X_n-X_1|| > x\right)$ is given using formula $(6)$ from this paper as $\int_x^{\sqrt{2}}g(t)dt$ where $$g(x) = \begin{cases} 2x^{3}-8x^{2}+2\pi x & 0 \leq x < 1 \\ -2x^{3}-\left(2\pi+4\right)x+8x\operatorname{arccsc}\left(x\right)+8x\sqrt{x^{2}-1} & 1\leq x < \sqrt{2} \end{cases} \tag 5$$

Therefore, $P\left(||X_n-X_1|| > x\right)$ is given by $$\begin{cases} -\frac{x^{4}}{2}+\frac{8x^{3}}{3}-\pi x^{2}+1 & 0 \leq x < 1 \\ \frac{2}{3}-\frac{4}{3}\left(2x^{2}+1\right)\sqrt{x^{2}-1}+\frac{1}{2}x^{4}+\left(2+\pi\right)x^{2}-4x^{2}\operatorname{arccsc}\left(x\right) & 1 \leq x < \sqrt{2} \end{cases} \tag 6$$

The integral $\int_0^{\sqrt{2}}P\left(||X_n-X_1|| > x\right)^{n-1}dx$ is unlikely to have a closed form (or if it does, it would probably be very complicated). Instead, for an approximation, split it up as $$\int_0^{1}P\left(||X_n-X_1|| > x\right)^{n-1}dx + \int_1^{\sqrt{2}}P\left(||X_n-X_1|| > x\right)^{n-1}dx$$

The second integral has an upper bound of $\left(\frac{19}{6}-\pi\right)^{n-1}\left(\sqrt{2}-1\right)$, obtained because each value in that interval has an upper bound of $P\left(||X_n-X_1|| > 1\right)^{n-1} = \left(\frac{19}{6}-\pi\right)^{n-1}$. The first integral is given by $$\sum_{k=0}^{n-1}\sum_{m=0}^{n-k-1}\sum_{r=0}^{n-k-1}\frac{\binom{n-1}{k, m, r}}{\left(1+2k+3m+4r\right)}\left(-\pi\right)^{k}\left(\frac{8}{3}\right)^{m}\left(-\frac{1}{2}\right)^{r}$$

Numerically, this approaches $\frac{1}{2\sqrt{n-1}}$, but I'm having trouble proving it.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.