Expected minimum distance from a point with varying density 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


#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: 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).$$

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.

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")

