Apa keuntungan generator acak eksponensial menggunakan metode Ahrens dan Dieter (1972) daripada dengan inverse transform?

Pertanyaan saya terinspirasi oleh generator nomor acak eksponensial bawaan R , fungsinya rexp(). Ketika mencoba untuk menghasilkan angka acak yang didistribusikan secara eksponensial, banyak buku teks merekomendasikan metode transformasi terbalik seperti yang diuraikan dalam halaman Wikipedia ini . Saya sadar bahwa ada metode lain untuk menyelesaikan tugas ini. Secara khusus, kode sumber R menggunakan algoritma yang digariskan dalam sebuah makalah oleh Ahrens & Dieter (1972) .

Saya telah meyakinkan diri saya sendiri bahwa metode Ahrens-Dieter (AD) benar. Namun, saya tidak melihat manfaat dari menggunakan metode mereka dibandingkan dengan metode inverse transform (IT). AD tidak hanya lebih kompleks untuk diterapkan daripada IT. Sepertinya tidak ada manfaat kecepatan juga. Berikut adalah kode R saya untuk membandingkan kedua metode yang diikuti oleh hasilnya.

invTrans <- function(n)
    -log(runif(n))
print("For the inverse transform:")
print(system.time(invTrans(1e8)))
print("For the Ahrens-Dieter algorithm:")
print(system.time(rexp(1e8)))

Hasil:

[1] "For the inverse transform:" 
user     system     elapsed
4.227    0.266      4.597 
[1] "For the Ahrens-Dieter algorithm:"
user     system     elapsed
4.919    0.265      5.213

Membandingkan kode untuk dua metode, AD menggambar setidaknya dua angka acak seragam (dengan fungsi Cunif_rand() ) untuk mendapatkan satu nomor acak eksponensial. TI hanya membutuhkan satu nomor acak yang seragam. Agaknya tim inti R memutuskan untuk tidak mengimplementasikan TI karena diasumsikan bahwa mengambil logaritma mungkin lebih lambat daripada menghasilkan angka acak yang lebih seragam. Saya mengerti bahwa kecepatan mengambil logaritma dapat bergantung pada mesin, tetapi setidaknya bagi saya yang sebaliknya adalah benar. Mungkin ada masalah di sekitar ketepatan numerik IT yang berkaitan dengan singularitas logaritma pada 0? Tapi kemudian, kode sumber R sexp.cmengungkapkan bahwa implementasi AD juga kehilangan beberapa ketepatan numerik karena bagian berikut dari kode C menghilangkan bit-bit utama dari nomor acak seragam u .

double u = unif_rand();
while(u <= 0. || u >= 1.) u = unif_rand();
for (;;) {
    u += u;
    if (u > 1.)
        break;
    a += q[0];
}
u -= 1.;

u kemudian didaur ulang sebagai nomor acak seragam dalam sisa sexp.c . Sejauh ini, tampak seolah-olah

IT lebih mudah untuk dikodekan,
TI lebih cepat, dan
IT dan AD mungkin kehilangan akurasi numerik.

Saya akan sangat menghargai jika seseorang dapat menjelaskan mengapa R masih mengimplementasikan AD sebagai satu-satunya pilihan yang tersedia untuk rexp().

r simulation random-generation exponential inverse-cdf Pratyush Lebih
sumber

Dengan generator angka acak, "kode lebih mudah" tidak benar-benar menjadi pertimbangan kecuali Anda yang melakukannya! Kecepatan dan ketepatan adalah dua pertimbangan saja. (Untuk generator yang seragam, ada juga periode generator.) Di masa lalu, AD lebih cepat. Pada kotak Linux saya, AD berjalan sekitar 1/2 waktu fungsi invTrans Anda, dan pada laptop saya sekitar 2/3 waktu. Anda mungkin ingin menggunakan microbenchmark untuk pengaturan waktu yang lebih komprehensif juga.

jbowman

Saya menyarankan agar kita tidak memigrasikannya. Ini sepertinya menjadi topik bagi saya.

Amoeba berkata Reinstate Monica

Mengingat bahwa saya tidak dapat memikirkan satu skenario di mana rexp(n)akan menjadi hambatan, perbedaan dalam kecepatan bukanlah argumen yang kuat untuk perubahan (setidaknya bagi saya). Saya mungkin lebih peduli tentang akurasi numerik, meskipun tidak jelas bagi saya yang mana yang lebih dapat diandalkan secara numerik.

Cliff AB

@amoeba Saya berpikir bahwa "Apa yang akan menjadi keuntungan dari ..." akan menjadi pengulangan kata yang akan jelas-topik di sini, dan tidak akan mempengaruhi jawaban yang ada. Saya kira "Mengapa orang-orang yang membuat R memutuskan untuk melakukan ..." sebenarnya adalah (a) pertanyaan khusus perangkat lunak, (b) memerlukan bukti dalam dokumentasi atau telepati, sehingga bisa dibilang di luar topik di sini. Secara pribadi saya lebih suka pertanyaan itu diulang untuk membuatnya lebih jelas dalam lingkup situs, tetapi saya tidak melihat ini sebagai alasan yang cukup kuat untuk menutupnya.

Silverfish

@amoeba Saya sudah mencobanya. Tidak yakin judul yang saya sarankan baru terutama tata bahasa, dan mungkin beberapa bagian lain dari teks pertanyaan bisa dilakukan dengan perubahan. Tapi saya harap ini lebih jelas pada topik, setidaknya, dan saya tidak berpikir itu membatalkan atau memerlukan perubahan untuk salah satu jawaban.

Silverfish

Jawaban:

Di komputer saya (maafkan bahasa Prancis saya!):

> print(system.time(rexp(1e8)))
utilisateur     système      écoulé 
      4.617       0.320       4.935 
> print(system.time(rexp(1e8)))
utilisateur     système      écoulé 
      4.589       2.045       6.629 
> print(system.time(-log(runif(1e8))))
utilisateur     système      écoulé 
      7.455       1.080       8.528 
> print(system.time(-log(runif(1e8))))
utilisateur     système      écoulé 
      9.140       1.489      10.623

transformasi terbalik tidak lebih buruk. Tetapi Anda harus berhati-hati terhadap variabilitas. Memperkenalkan parameter laju mengarah ke variabilitas yang lebih besar untuk transformasi terbalik:

> print(system.time(rexp(1e8,rate=.01)))
utilisateur     système      écoulé 
      4.594       0.456       5.047 
> print(system.time(rexp(1e8,rate=.01)))
utilisateur     système      écoulé 
      4.661       1.319       5.976 
> print(system.time(-log(runif(1e8))/.01))
utilisateur     système      écoulé 
     15.675       2.139      17.803 
> print(system.time(-log(runif(1e8))/.01))
utilisateur     système      écoulé 
      7.863       1.122       8.977 
> print(system.time(rexp(1e8,rate=101.01)))
utilisateur     système      écoulé 
      4.610       0.220       4.826 
> print(system.time(rexp(1e8,rate=101.01)))
utilisateur     système      écoulé 
      4.621       0.156       4.774 
> print(system.time(-log(runif(1e8))/101.01))
utilisateur     système      écoulé 
      7.858       0.965       8.819 > 
> print(system.time(-log(runif(1e8))/101.01))
utilisateur     système      écoulé 
     13.924       1.345      15.262

Berikut perbandingan yang digunakan rbenchmark:

> benchmark(x=rexp(1e6,rate=101.01))
  elapsed user.self sys.self
  4.617     4.564    0.056
> benchmark(x=-log(runif(1e6))/101.01)
  elapsed user.self sys.self
  14.749   14.571    0.184
> benchmark(x=rgamma(1e6,shape=1,rate=101.01))
  elapsed user.self sys.self
  14.421   14.362    0.063
> benchmark(x=rexp(1e6,rate=.01))
  elapsed user.self sys.self
  9.414     9.281    0.136
> benchmark(x=-log(runif(1e6))/.01)
  elapsed user.self sys.self
  7.953     7.866    0.092
> benchmark(x=rgamma(1e6,shape=1,rate=.01))
  elapsed user.self sys.self
  26.69    26.649    0.056

Jadi jarak tempuh masih bervariasi, tergantung pada skalanya!

Xi'an
sumber

Di laptop saya, waktunya cocok dengan OP sangat dekat sehingga saya curiga kami memiliki mesin yang sama (atau setidaknya prosesor yang sama). Tapi saya pikir poin Anda di sini adalah keunggulan kecepatan yang diamati tergantung pada platform, dan mengingat perbedaan minimal, tidak ada keuntungan yang jelas bagi satu sama lain dalam hal kecepatan.

Cliff AB

Bisakah Anda melakukan yang microbenchmarksebaliknya?

Firebug

Waktu sistem tampaknya mengukur overhead yang sangat bervariasi, mungkin karena gangguan dan paging memori. Sangat menarik, seperti dicatat oleh @Cliff, untuk melihat perbedaan relatif yang cukup besar dalam kinerja antar sistem. Pada Xeon dengan banyak RAM, misalnya, saya mengeluarkan hampir tidak ada waktu sistem (0,05 hingga 0,32 detik), sekitar 12% lebih lama waktu pengguna untuk rexp, 3% waktu pengguna lebih pendek untuk -log(runif()), dan tidak ada variabilitas dengan parameter laju ( total detik). Kita semua secara implisit berasumsi mencapai waktu untuk dan sebanding dengan apa yang akan diperoleh seseorang dengan subrutin Fortran.

5.27 \pm 0.02

$5.27\pm0.02$ Rlogrunif

Whuber

Ini hanya mengutip artikel di bawah bagian "Algoritma LG: (Metode logaritma)":

Dalam FORTRAN algoritma paling baik diprogram secara langsung sebagai menghindari panggilan subprogram. Kinerja ini 504 sec per sampel. Dari waktu ini 361 diambil oleh rutin logaritma pabrikan dan 105 mtk oleh generator untuk REGOL dari variabel terdistribusi seragam . Oleh karena itu tidak ada gunanya mencoba meningkatkan kecepatan dengan menulis fungsi assembler yang harus menggunakan dua subprogram yang sama. $X=-ALOG(REGOL(IR))$ $\mu$ $\mu$ $\mu$ $u$

Jadi sepertinya penulis memilih metode lain untuk menghindari pembatasan "pabrikan" logaritma lambat ini. Mungkin pertanyaan ini kemudian lebih baik dipindahkan ke stackoverflow di mana seseorang dengan pengetahuan tentang nyali R dapat berkomentar.

Alex R.
sumber

Hanya menjalankan ini dengan microbenchmark; pada mesin saya, pendekatan asli R lebih cepat:

library(microbenchmark)
microbenchmark(times = 10L,
               R_native = rexp(1e8),
               dir_inv = -log(runif(1e8)))
# Unit: seconds
#      expr      min       lq     mean   median       uq      max neval
#  R_native 3.643980 3.655015 3.687062 3.677351 3.699971 3.783529    10
#   dir_inv 5.780103 5.783707 5.888088 5.912384 5.946964 6.050098    10

Demi kebaruan, inilah yang memastikan itu tidak sepenuhnya karena memiliki : $\lambda = 1$

lambdas = seq(0, 10, length.out = 25L)[-1L]
png("~/Desktop/micro.png")
matplot(lambdas, 
        ts <- 
          t(sapply(lambdas, function(ll)
            print(microbenchmark(times = 50L,
                                 R_native = rexp(5e5, rate = ll),
                                 dir_inv = -log(runif(5e5))/ll),
                  unit = "relative")[ , "median"])),
        type = "l", lwd = 3L, xlab = expression(lambda),
        ylab = "Relative Timing", lty = 1L,
        col = c("black", "red"), las = 1L,
        main = paste0("Direct Computation of Exponential Variates\n",
                      "vs. R Native Generator (Ahrens-Dieter)"))
text(lambdas[1L], ts[1L, ], c("A-D", "Direct"), pos = 3L)
dev.off()

MichaelChirico
sumber