Distance minimale attendue d'un point de densité variable

8

J'examine comment la distance euclidienne minimale attendue entre des points uniformément aléatoires et l'origine change à mesure que nous augmentons la densité de points aléatoires ( points par unité de carré ) autour de l'origine. J'ai réussi à trouver une relation entre les deux décrits comme tels:

Distance minimale attendue=12Densité

J'ai trouvé cela en exécutant des simulations de Monte Carlo en R et en ajustant une courbe manuellement (code ci-dessous).

Ma question est : aurais-je pu dériver ce résultat théoriquement plutôt que par l'expérimentation?

#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

entrez la description de l'image ici

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

Un graphique des valeurs de densité par rapport au minimum attendu, à la fois à Monte Carlo et théorique

Michael Bird
la source
La dépendance directe (asymptotique) à la racine inverse de la densité découle facilement et immédiatement des considérations sur les unités de mesure, de sorte que la seule question concerne la raison pour laquelle le multiple est 1/2.
whuber
@whuber Oui, j'avais remarqué que les unités étaient bien alignées et oui, la question devient: d'où viennent les 2?
Michael Bird
1
le 2est la largeur de votre carré.
whuber

Réponses:

8

Considérez la distance à l'origine de n variables aléatoires distribuées indépendamment (Xje,Ouije) qui ont des distributions uniformes sur le carré [-1,1]2.

L'écriture Rje2=Xje2+Ouije2 pour la distance au carré, la géométrie euclidienne nous montre que

Pr(Rjer1)=14πr2

tout (avec un peu plus de travail)

Pr(1Rjer2)=14(πr2+4r2-1-4r2ArcTan(r2-1)).

Figure 1: Tracé de la fonction de distribution

Ensemble, ceux-ci déterminent la fonction de distribution commune à tous lesFRje.

Parce que les points sont indépendants, les distances le sont où la fonction de survie de estnRje,min(Rje)

Sn(r)=(1-F(r))n,

impliquant la distance la plus courte moyenne est

μ(n)=02Sn(r)r.

Pour presque toute l'aire de cette intégrale est proche de nous pouvons donc l'approcher commen1,0,

μenviron(n)=01Sn(r)r=01(1-π4r2)nr.

L'erreur n'est pas supérieure à la partie de l'intégrale qui a été omise, qui n'est à son tour pas supérieure à

(2-1)(1-F(1))n=(2-1)(1-π/4)n,

ce qui diminue évidemment exponentiellement avecn.

On peut à son tour approximer l'intégrale comme

(1-π4r2)nexp(-12r22/(nπ)).

Jusqu'à une constante de normalisation, il s'agit de la fonction de densité d'une distribution normale avec une moyenne de et une variance La constante de normalisation manquante est0σ2=2/(nπ).

C(n)=12πσ2=12π 2/(nπ)=n2.

Par conséquent, l'extension de l'intégrale de1 à (ce qui ajoute une erreur proportionnelle à ),e-n

μenviron(n)0e-t2/(2σ2)t=1C(n)12=1n.

Dans le processus d'obtention de cette approximation, trois erreurs ont été commises. Collectivement, ils sont au plus d'ordren-1, l'erreur encourue lors de l'approximation de par le gaussien.Sn(r)

! [Figure 2: Tracé des erreurs de simulation

Cette figure représente fois la différence entre et fois la distance moyenne la plus courte observée dans ensembles de données simulés distincts pour chaquen1ndix5n. Parce qu'ils diminuent à mesure que croît, cela prouve que l'erreur estno(n-1/n)=o(n-3/2).

Enfin, le facteur dans la question découle de la taille du carré:1/2 la densité est le nombre de points par unité de surface et le carrén,[-1,1]2 a la surface , d'où4

2Densité=2n/4=n.

Voici le Rcode de la 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")
whuber
la source
2
Hou la la! Quelle réponse! Merci beaucoup, c'est super. Merci!
Michael Bird
Bonjour @whuber, j'essayais de reproduire votre et j'ai remarqué que votre équation pour ne renvoie pas comme le montrent vos graphiques. Lorsque j'ai calculé j'ai obtenu qui donne la courbe que vous avez fournie. Avez-vous fait une faute de frappe? F(r)F(2)1Pr(1Rjer2)π/4-r(rArcCos(1/r)-1-1/r2)
Michael Bird
1
@Michael Merci, il y a une faute de frappe - mais ce n'est pas celle que vous suggérez: un de mes " " aurait dû être " ". J'ai corrigé celui-là. r4
whuber