Jarak minimum yang diharapkan dari suatu titik dengan kepadatan yang bervariasi

8

Saya melihat bagaimana jarak Euclidean minimum yang diharapkan antara titik-titik seragam acak dan titik asal berubah saat kami meningkatkan kepadatan titik-titik acak ( titik per unit persegi ) di sekitar titik asal. Saya telah berhasil menemukan hubungan antara keduanya yang digambarkan sebagai berikut:

Expected Min Distance=12Massa jenis

Saya datang dengan ini dengan menjalankan beberapa simulasi Monte Carlo di R dan pas kurva secara manual (kode di bawah).

Pertanyaan saya adalah : dapatkah saya memperoleh hasil ini secara teoritis daripada melalui eksperimen?

#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

masukkan deskripsi gambar di sini

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

Grafik nilai kepadatan terhadap minimum yang diharapkan, baik monte carlo maupun teoretis

Michael Bird
sumber
Ketergantungan langsung (asimptotik) pada akar kebalikan mengikuti dengan mudah dan langsung dari pertimbangan unit pengukuran, sehingga satu-satunya pertanyaan menyangkut mengapa kelipatannya adalah1/2.
Whuber
@whuber Ya saya perhatikan unit berbaris dengan baik dan ya, pertanyaannya menjadi: dari mana 2 berasal?
Michael Bird
1
Angka adalah lebar kotak Anda. 2
Whuber

Jawaban:

8

Pertimbangkan jarak ke asal variabel acak yang didistribusikan secara independen yang memiliki distribusi seragam di alun-alunn(Xsaya,Ysaya)[-1,1]2.

Menulis untuk jarak kuadrat, geometri Euclidean menunjukkan kepada kita bahwaRsaya2=Xsaya2+Ysaya2

Pr(Rsayar1)=14πr2

sementara (dengan sedikit lebih banyak pekerjaan)

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

Gambar 1: Plot fungsi distribusi

Bersama-sama ini menentukan fungsi distribusi umum untuk semuaFRsaya.

Karena titik independen, begitu juga jarak mana fungsi survival adalahnRsaya,min(Rsaya)

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

menyiratkan jarak terpendek rata-rata adalah

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

Untuk hampir semua area dalam integral ini mendekati jadi kami dapat memperkirakannya sebagain1,0,

μsekitar(n)=01Sn(r)dr=01(1-π4r2)ndr.

Kesalahan tidak lebih besar dari bagian integral yang dihilangkan, yang pada gilirannya tidak lebih besar dari

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

yang jelas menurun secara eksponensial dengann.

Kami mungkin akan memperkirakan integrand sebagai

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

Hingga konstanta normalisasi, ini adalah fungsi kerapatan distribusi Normal dengan rata-rata dan varians Konstanta normalisasi yang hilang adalah0σ2=2/(nπ).

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

Oleh karena itu, memperluas integral dari hingga (yang menambahkan kesalahan sebanding dengan ),1e-n

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

Dalam proses memperoleh perkiraan ini tiga kesalahan dibuat. Secara kolektif mereka paling banyak memesan kesalahan terjadi ketika mendekati oleh Gaussian.n-1,Sn(r)

! [Gambar 2: Plot kesalahan simulasi

Gambar ini memplot kali perbedaan antara dan kali jarak terpendek rata-rata yang diamati dalam set data simulasi terpisah untuk setiap Karena berkurang ketika bertambah, ini adalah bukti bahwa kesalahannya adalahn1n105n.nHai(n-1/n)=Hai(n-3/2).

Akhirnya, faktor dalam pertanyaan berasal dari ukuran kuadrat:1/2 kepadatan adalah jumlah titik per satuan luas dan kuadrat memiliki luas , dari manan,[-1,1]24

2Massa jenis=2n/4=n.

Ini adalah Rkode untuk simulasi:

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
sumber
2
Wow! Jawaban yang sangat bagus! Terima kasih banyak, ini bagus. Terima kasih!
Michael Bird
Hai @whuber, saya mencoba mereproduksi dan saya perhatikan persamaan Anda untuk tidak menghasilkan seperti yang ditunjukkan grafik Anda. Ketika saya menghitung Saya mendapat yang memberikan kurva yang Anda berikan. Apakah Anda membuat kesalahan ketik? F(r)F(2)1Pr(1Rsayar2)π/4-r(rArcCos(1/r)-1-1/r2)
Michael Bird
1
@Michael Terima kasih, ada kesalahan ketik - tapi itu bukan yang Anda sarankan: salah satu dari " " saya seharusnya " " Saya telah memperbaiki yang itu. r4
whuber