8
$\begingroup$

I'm looking at how the expected minimum Euclidean distance between randomly uniform points and the origin changes as we increase the density of random points (points per unit square) around the origin. I have managed to come up with a relationship between the two described as such:

$$\text{Expected Min Distance} =\frac{1}{2\sqrt{\text{Density}}}$$

I came up with this by running some Monte Carlo simulations in R and fitting a curve manually (code below).

My question is: could I have derived this result theoretically rather than through experimentation?

#Stack Overflow example
library(magrittr)
library(ggplot2)


#---------
#FUNCTIONS
#---------
#gen random points within a given radius and given density
gen_circle_points <- function(radius, density) {
  #round radius up then generate points in square with side length = 2*radius
  c_radius <- ceiling(radius)
  coords <- data.frame(
    x = runif((2 * c_radius) ^ 2 * density, -c_radius, c_radius),
    y = runif((2 * c_radius) ^ 2 * density, -c_radius, c_radius)
  )
  return(coords[sqrt(coords$x ^ 2 + coords$y ^ 2) <= radius, ])#filter in circle
}

#Example plot
plot(gen_circle_points(radius = 1,density = 200)) #200 points around origin
points(0,0, col="red",pch=19) #colour origin

enter image description here

#return euclidean distances of points generated by gen_circle_points()
calculate_distances <- function(circle_points) {
  return(sqrt(circle_points$x ^ 2 + circle_points$y ^ 2))
}

#find the smallest distance from output of calculate_distances()
calculate_min_value <- function(distances) {
  return(min(distances))
}


#Try a range of values
density_values <- c(1:100)

expected_min_from_density <- sapply(density_values, function(density) {
  #simulate each density value 1000 times and take an average as estimate for
  #expected minimum distance
  sapply(1:1000, function(i) {
    gen_circle_points(radius=1, density=density) %>%
      calculate_distances() %>%
      calculate_min_value()
  }) %>% mean()
})

results <- data.frame(density_values, expected_min_from_density)

#fit based off exploration
theoretical_fit <- data.frame(density = density_values, 
                              fit = 1 / (sqrt(density_values) * 2))

#plot monte carlo (black) and fit (red dashed)
ggplot(results, aes(x = density_values, y = expected_min_from_density)) +
  geom_line() + 
  geom_line(
    data = theoretical_fit,
    aes(x = density, y = fit),
    color = "red",
    linetype = 2
  )

A graph of density values against expected minimum, both monte carlo and theoretical

$\endgroup$
  • $\begingroup$ The (asymptotic) direct dependence on the inverse root of density follows easily and immediately from considerations of the units of measurement, so the only question concerns why the multiple is $1/2.$ $\endgroup$ – whuber Apr 16 '18 at 13:34
  • $\begingroup$ @whuber Yes I'd noticed the units lined up nicely and yes, the question becomes: where did the 2 come from? $\endgroup$ – Michael Bird Apr 16 '18 at 13:49
  • 1
    $\begingroup$ The $2$ is the width of your square. $\endgroup$ – whuber Apr 16 '18 at 14:47
8
$\begingroup$

Consider the distance to the origin of $n$ independently distributed random variables $(X_i,Y_i)$ that have uniform distributions on the square $[-1,1]^2.$

Writing $R_i^2 = X_i^2+Y_i^2$ for the squared distance, Euclidean geometry shows us that

$$\Pr(R_i \le r \le 1) = \frac{1}{4} \pi\, r^2 $$

while (with a little more work)

$$\Pr(1 \le R_i \le r \le \sqrt{2}) = \frac{1}{4}\left(\pi\, r^2 + 4\sqrt{r^2-1} - 4 r^2 \operatorname{ArcTan}\left(\sqrt{r^2-1}\right)\right).$$

Figure 1: Plot of the distribution function

Together these determine the distribution function $F$ common to all the $R_i.$

Because the $n$ points are independent, so are the distances $R_i,$ whence the survival function of $\min(R_i)$ is

$$S_n(r) = (1 - F(r))^n,$$

implying the mean shortest distance is

$$\mu(n) = \int_0^\sqrt{2} S_n(r)\, dr.$$

For $n\gg 1,$ almost all the area in this integral is close to $0,$ so we may approximate it as

$$\mu_\text{approx}(n) = \int_0^1S_n(r)\, dr = \int_0^1\left(1 - \frac{\pi}{4}r^2\right)^n\,dr.$$

The error is not greater than the part of the integral that was omitted, which is in turn no greater than

$$(\sqrt{2}-1)(1-F(1))^n = (\sqrt{2}-1)(1 - \pi/4)^n,$$

which obviously decreases exponentially with $n.$

We may in turn approximate the integrand as

$$\left(1 - \frac{\pi}{4}r^2\right)^n \approx \exp\left(-\frac{1}{2} \frac{r^2}{2/(n\pi)}\right).$$

Up to a normalizing constant, this is the density function of a Normal distribution with mean $0$ and variance $\sigma^2=2/(n\pi).$ The missing normalizing constant is

$$C(n) = \frac{1}{\sqrt{2\pi \sigma^2}} = \frac{1}{\sqrt{2\pi\ 2 / (n\pi)}} = \frac{\sqrt{n}}{2}.$$

Therefore, extending the integral from $1$ to $\infty$ (which adds an error proportional to $e^{-n}$),

$$\mu_\text{approx}(n) \approx \int_0^\infty e^{-t^2/(2\sigma^2)}\,dt = \frac{1}{C(n)} \frac{1}{2} = \frac{1}{\sqrt{n}}.$$

In the process of obtaining this approximation three errors were made. Collectively they are at most of order $n^{-1},$ the error incurred when approximating $S_n(r)$ by the Gaussian.

![Figure 2: Plot of simulation errors

This figure plots $n$ times the difference between $1$ and $\sqrt{n}$ times the mean shortest distance observed in $10^5$ separate simulated datasets for each $n.$ Because they decrease as $n$ grows, this is evidence that the error is $o(n^{-1}/\sqrt{n}) = o(n^{-3/2}).$

Finally, the factor $1/2$ in the question derives from the size of the square: the density is the number of points $n,$ per unit area and the square $[-1,1]^2$ has area $4$, whence

$$2\sqrt{\text{Density}} = 2\sqrt{n/4} = \sqrt{n}.$$


This is the R code for the simulation:

n.sim <- 1e5  # Size of each simulation
d <- 2        # Dimension
n <- 2^(1:11) # Numbers of points in each simulation
#
# Estimate mean distance to the origin for each `n`.
#
y <- sapply(n, function(n.points) {
  x <- array(runif(d*n.points*n.sim, -1, 1), c(d, n.points, n.sim))
  mean(sqrt(apply(colSums(x^2), 2, min)))
})
#
# Plot the errors (normalized) against `n`.
#
library(ggplot2)
ggplot(data.frame(Log2.n = 1:length(n), Error=sqrt(n)* (1 - y * n^(1/d))),
       aes(Log2.n, Error)) + geom_point() + geom_smooth() 
  ylab("Error * n") + ggtitle("Simulation Means")
$\endgroup$
  • 2
    $\begingroup$ Wow! What an answer! Thanks a lot, this is great. Thanks! $\endgroup$ – Michael Bird Apr 16 '18 at 15:50
  • $\begingroup$ Hi @whuber, I was trying to reproduce your $F(r)$ and I noticed your equation for $F(\sqrt{2})$ doesn't return $1$ like your graphs shows. When I calculated $\text{Pr}(1 \leq R_i \leq r \leq \sqrt{2})$ I got $\pi/4 - r (r \text{ArcCos}(1/r) - \sqrt{1-1/r^2})$ which gives the curve you provided. Have you made a typo? $\endgroup$ – Michael Bird May 2 '18 at 15:03
  • 1
    $\begingroup$ @Michael Thank you, there is a typo--but it's not the one you suggest: one of my "$r$" should have been "$4$." I have fixed that one. $\endgroup$ – whuber May 2 '18 at 16:58

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.