Hasilkan titik secara efisien antara satuan lingkaran dan satuan bujur sangkar

17

Saya ingin menghasilkan sampel dari wilayah biru yang didefinisikan di sini:

masukkan deskripsi gambar di sini

Solusi naif adalah menggunakan sampel penolakan di unit square, tetapi ini hanya memberikan efisiensi (~ 21,4%).1π/4

Apakah ada cara agar saya dapat mencicipi lebih efisien?

Cam.Davidson.Pilon
sumber
6
Petunjuk : Gunakan simetri untuk menggandakan efisiensi Anda secara sepele.
kardinal
3
Oh suka: jika nilainya (0,0), ini bisa dipetakan ke (1,1)? Saya suka ide itu
Cam.Davidson.Pilon
@ cardinal Tidakkah seharusnya 4x efisiensinya? Anda dapat mengambil sampel dalam dan kemudian mencerminkannya pada sumbu x, sumbu y dan asal. [0,,1]×[0,,1]
Martin Krämer
1
@ Martin: Di keempat wilayah simetris, Anda memiliki tumpang tindih, yang harus Anda hadapi dengan lebih hati-hati.
kardinal
3
@ Martin: Jika saya memahami apa yang Anda gambarkan, itu tidak meningkatkan efisiensi sama sekali . (Anda menemukan satu titik, dan sekarang tahu tiga lainnya --- di daerah empat kali ukuran --- yang melakukan atau tidak terletak di dalam unit disk dengan probabilitas satu sesuai dengan apakah tidak. apakah itu membantu?) Maksud peningkatan efisiensi adalah meningkatkan kemungkinan penerimaan untuk setiap ( x , y ) yang dihasilkan. Mungkin saya orang yang padat? (x,y)(x,y)
kardinal

Jawaban:

10

Apakah dua juta poin per detik akan berhasil?

Distribusinya simetris: kita hanya perlu mengerjakan distribusi untuk seperdelapan dari lingkaran penuh dan kemudian menyalinnya di sekitar oktan lainnya. Dalam koordinat kutub , distribusi kumulatif sudut Θ untuk lokasi acak ( X , Y ) pada nilai θ diberikan oleh area antara segitiga ( 0 , 0 ) , ( 1 , 0 ) , ( 1 , tan θ ) dan lengkungan lingkaran memanjang dari ((r,θ)Θ(X,Y)θ(0,0),(1,0),(1,tanθ) hingga ( cos θ , sin θ ) . Dengan demikian sebanding dengan(1,0)(cosθ,sinθ)

FΘ(θ)=Pr(Θθ)12tan(θ)θ2,

dari mana kepadatannya

fΘ(θ)=ddθFΘ(θ)tan2(θ).

Kami dapat mengambil sampel dari kerapatan ini menggunakan, katakanlah, metode penolakan (yang memiliki efisiensi ).8/π254.6479%

Densitas kondisional dari koordinat radial sebanding dengan r d r antara r = 1 dan r = sec θ . Itu dapat disampel dengan inversi CDF yang mudah.Rrdrr=1r=secθ

Jika kita menghasilkan sampel independen , konversi kembali ke koordinat Cartesian ( x i , y i ) sampel oktan ini. Karena sampel independen, swapping koordinat secara acak menghasilkan sampel acak independen dari kuadran pertama, seperti yang diinginkan. (Swap acak memerlukan hanya menghasilkan variabel Binomial tunggal untuk menentukan berapa banyak realisasi untuk bertukar.)(ri,θi)(xi,yi)

Setiap realisasi membutuhkan rata-rata satu varian seragam (untuk R ) ditambah 1 / ( 8 π - 2 ) kali dua varian seragam (untuk Θ ) dan sejumlah kecil perhitungan (cepat). Itu 4 / ( π - 4 ) 4,66 variasi per titik (yang, tentu saja, memiliki dua koordinat). Rincian lengkap ada pada contoh kode di bawah ini. Angka ini memplot 10.000 dari lebih dari setengah juta poin yang dihasilkan.(X,Y)R1/(8π2)Θ4/(π4)4.66

Angka

Berikut adalah Rkode yang menghasilkan simulasi ini dan menghitung waktunya.

n.sim <- 1e6
x.time <- system.time({
  # Generate trial angles `theta`
  theta <- sqrt(runif(n.sim)) * pi/4
  # Rejection step.
  theta <- theta[runif(n.sim) * 4 * theta <= pi * tan(theta)^2]
  # Generate radial coordinates `r`.
  n <- length(theta)
  r <- sqrt(1 + runif(n) * tan(theta)^2)
  # Convert to Cartesian coordinates.
  # (The products will generate a full circle)
  x <- r * cos(theta) #* c(1,1,-1,-1)
  y <- r * sin(theta) #* c(1,-1,1,-1)
  # Swap approximately half the coordinates.
  k <- rbinom(1, n, 1/2)
  if (k > 0) {
    z <- y[1:k]
    y[1:k] <- x[1:k]
    x[1:k] <- z
  }
})
message(signif(x.time[3] * 1e6/n, 2), " seconds per million points.")
#
# Plot the result to confirm.
#
plot(c(0,1), c(0,1), type="n", bty="n", asp=1, xlab="x", ylab="y")
rect(-1, -1, 1, 1, col="White", border="#00000040")
m <- sample.int(n, min(n, 1e4))
points(x[m],y[m], pch=19, cex=1/2, col="#0000e010")
whuber
sumber
1
Saya tidak mengerti kalimat ini: "Karena sampel independen, secara sistematis menukar koordinat setiap sampel kedua menghasilkan sampel acak independen dari kuadran pertama, seperti yang diinginkan." Sepertinya saya bahwa secara sistematis menukar koordinat setiap sampel kedua menghasilkan sampel yang sangat tergantung. Misalnya, menurut saya implementasi Anda dalam kode menghasilkan setengah juta sampel berturut-turut dari oktan yang sama?
A. Rex
7
Sebenarnya, pendekatan ini tidak cukup berhasil (untuk titik awal) karena menghasilkan jumlah sampel yang identik dalam dua oktan: Titik sampel, karenanya, tergantung. Sekarang, jika Anda membalik koin yang tidak bias untuk menentukan oktan untuk setiap sampel ...
kardinal
1
@ Cardinal Anda benar; Saya akan memperbaikinya - tanpa (asimtotik) meningkatkan jumlah varian acak untuk menghasilkan!
whuber
2
Secara tegas (dan, sekali lagi, hanya dalam pengertian teoretis paling murni), dalam kasus sampel terbatas, modifikasi Anda tidak memerlukan variasi acak seragam yang seragam. Intinya: Dari varian acak seragam yang pertama, buat urutan flipping dari bit pertama . Kemudian, gunakan sisa (kali 2 n ) sebagai koordinat pertama yang dihasilkan. n2n
kardinal
2
@ Xi'an Saya tidak dapat memperoleh invers yang dapat dihitung dengan mudah. Saya dapat melakukan sedikit lebih baik dengan penolakan sampel dari distribusi dengan kepadatan sebanding dengan (efisiensinya adalah ( 4 - π ) / ( π - 2 ) 75 % ), dengan biaya harus menghitung arcsine . 2sin(θ)2(4π)/(π2)75%
whuber
13

Saya mengusulkan solusi berikut, yang seharusnya lebih sederhana, lebih efisien dan / atau lebih murah secara komputasi daripada sesi lainnya oleh @ cardinal, @whuber dan @ stephan-kolassa sejauh ini.

Ini melibatkan langkah-langkah sederhana berikut:

u1Unif(0,1)u2Unif(0,1).

min{u1,u2},max{u1,u2}

[xy]=[11]+[2212210][min{u1,u2}max{u1,u2}].

xyu1>u2

x2+y2<1.

Intuisi di balik algoritma ini ditunjukkan pada gambar. masukkan deskripsi gambar di sini

Langkah 2a dan 2b dapat digabungkan menjadi satu langkah:

x=1+22min(u1,u2)u2y=1+22min(u1,u2)u1

Kode berikut mengimplementasikan algoritma di atas (dan mengujinya menggunakan kode @ whuber).

n.sim <- 1e6
x.time <- system.time({
    # Draw two standard uniform samples
    u_1 <- runif(n.sim)
    u_2 <- runif(n.sim)
    # Apply shear transformation and swap
    tmp <- 1 + sqrt(2)/2 * pmin(u_1, u_2)
    x <- tmp - u_2
    y <- tmp - u_1
    # Reject if inside circle
    accept <- x^2 + y^2 > 1
    x <- x[accept]
    y <- y[accept]
    n <- length(x)
})
message(signif(x.time[3] * 1e6/n, 2), " seconds per million points.")
#
# Plot the result to confirm.
#
plot(c(0,1), c(0,1), type="n", bty="n", asp=1, xlab="x", ylab="y")
rect(-1, -1, 1, 1, col="White", border="#00000040")
m <- sample.int(n, min(n, 1e4))
points(x[m],y[m], pch=19, cex=1/2, col="#0000e010")

Beberapa tes cepat menghasilkan hasil berikut.

Algoritma /stats//a/258349 . Terbaik 3: 0,33 detik per juta poin.

Algoritma ini. Terbaik 3: 0,18 detik per juta poin.

Luca Citi
sumber
3
+1 Sangat bagus! Terima kasih telah berbagi solusi yang bijaksana, cerdas, dan sederhana.
whuber
Ide yang hebat! Saya sedang berpikir tentang pemetaan dari unit persegi ke bagian ini, tetapi tidak memikirkan pemetaan yang tidak sempurna dan kemudian skema penolakan. Terima kasih telah memperluas pikiran saya!
Cam.Davidson.Pilon
7

Ya, lebih efisien bisa dilakukan, tapi saya harap Anda tidak mencari lebih cepat .

xx

f(x)=11x2.

Wolfram membantu Anda untuk mengintegrasikan itu :

0xf(y)dy=12x1x2+x12arcsinx.

So the cumulative distribution function F would be this expression, scaled to integrate to 1 (i.e., divided by 01f(y)dy).

Now, to generate your x value, pick a random number t, uniformly distributed between 0 and 1. Then find x such that F(x)=t. That is, we need to invert the CDF (inverse transform sampling). This can be done, but it's not easy. Nor fast.

Finally, given x, pick a random y that is uniformly distributed between 1-x2 dan 1.

Di bawah ini adalah kode R. Perhatikan bahwa saya sedang mengevaluasi CDF di kisix nilai-nilai, dan itupun dibutuhkan beberapa menit.

Anda mungkin dapat mempercepat inversi CDF sedikit jika Anda berinvestasi beberapa pemikiran. Kemudian lagi, berpikir itu menyakitkan. Saya pribadi akan mengambil sampel penolakan, yang lebih cepat dan jauh lebih sedikit kesalahan, kecuali saya punya alasan yang sangat bagus untuk tidak melakukannya.

epsilon <- 1e-6
xx <- seq(0,1,by=epsilon)
x.cdf <- function(x) x-(x*sqrt(1-x^2)+asin(x))/2
xx.cdf <- x.cdf(xx)/x.cdf(1)

nn <- 1e4
rr <- matrix(nrow=nn,ncol=2)
set.seed(1)
pb <- winProgressBar(max=nn)
for ( ii in 1:nn ) {
    setWinProgressBar(pb,ii,paste(ii,"of",nn))
    x <- max(xx[xx.cdf<runif(1)])
    y <- runif(1,sqrt(1-x^2),1)
    rr[ii,] <- c(x,y)
}
close(pb)

plot(rr,pch=19,cex=.3,xlab="",ylab="")

randoms

S. Kolassa - Reinstate Monica
sumber
Saya ingin tahu apakah menggunakan polinomial Chebyshev untuk mendekati CDF akan meningkatkan kecepatan evaluasi.
Sycorax berkata Reinstate Monica
@ Scorax, bukan tanpa modifikasi; lihat misalnya pengobatan chebfun untuk singularitas aljabar di titik akhir.
JM bukan ahli statistik