Saya ingin menghasilkan sampel dari wilayah biru yang didefinisikan di sini:
Solusi naif adalah menggunakan sampel penolakan di unit square, tetapi ini hanya memberikan efisiensi (~ 21,4%).
Apakah ada cara agar saya dapat mencicipi lebih efisien?
probability
sampling
monte-carlo
random-generation
Cam.Davidson.Pilon
sumber
sumber
Jawaban:
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θ)
dari mana kepadatannya
Kami dapat mengambil sampel dari kerapatan ini menggunakan, katakanlah, metode penolakan (yang memiliki efisiensi ).8/π−2≈54.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.R rdr r=1 r=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) R 1/(8π−2) Θ 4/(π−4)≈4.66
Berikut adalah
R
kode yang menghasilkan simulasi ini dan menghitung waktunya.sumber
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:
Intuisi di balik algoritma ini ditunjukkan pada gambar.
Langkah 2a dan 2b dapat digabungkan menjadi satu langkah:
Kode berikut mengimplementasikan algoritma di atas (dan mengujinya menggunakan kode @ whuber).
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.
sumber
Ya, lebih efisien bisa dilakukan, tapi saya harap Anda tidak mencari lebih cepat .
Wolfram membantu Anda untuk mengintegrasikan itu :
So the cumulative distribution functionF would be this expression, scaled to integrate to 1 (i.e., divided by ∫10f(y)dy ).
Now, to generate yourx 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, givenx , 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.
sumber