Simulasikan konstrain normal pada batas bawah atau atas dalam R

Saya ingin menghasilkan data acak dari distribusi normal terbatas menggunakan R.

Sebagai contoh saya mungkin ingin mensimulasikan variabel dari distribusi normal dengan mean=3, sd= 2dan nilai apa pun yang lebih besar dari 5 diresampled dari distribusi normal yang sama.

Jadi, untuk fungsi umum, saya bisa melakukan hal berikut.

rnorm(n=100, mean=3, sd=2)

Saya kemudian memiliki beberapa pemikiran:

• Iterasikan ifelsefungsi dengan loop yang berulang hingga semua nilai dibatasi untuk berada di dalam batas.
• Mensimulasikan lebih banyak nilai daripada yang dibutuhkan dan mengambil yang pertama nyang memenuhi batasan.
• Hindari simulator variabel normal yang di-vektor-kan dan gunakan perulangan for dengan do sementara di dalam untuk mensimulasikan setiap pengamatan satu per satu dan perulangan jika diperlukan.

Semua hal di atas tampak agak kikuk.

Pertanyaan

• Apa cara sederhana untuk mensimulasikan variabel normal acak terbatas dalam R dari normal dengan mean = 3, sd = 2 dan max = 5?
• Lebih umum apa yang merupakan cara umum yang baik untuk menggabungkan kendala ke dalam variabel simulasi dalam R?

Ini disebut distribusi normal terpotong:

http://en.wikipedia.org/wiki/Truncated_normal_distribution

Christian Robert menulis tentang pendekatan untuk melakukannya untuk berbagai situasi (menggunakan berbeda tergantung di mana titik pemotongan berada) di sini:

Robert, CP (1995) "Simulasi variabel normal terpotong",
Statistik dan Komputasi, Volume 5, Edisi 2, Juni, hal 121-125

Kertas tersedia di http://arxiv.org/abs/0907.4010

Ini membahas sejumlah ide berbeda untuk titik pemotongan yang berbeda. Ini bukan satu-satunya cara untuk mendekati ini dengan cara apa pun tetapi memiliki kinerja yang cukup baik. Jika Anda ingin melakukan banyak norma terpotong yang berbeda dengan berbagai titik pemotongan, itu akan menjadi pendekatan yang masuk akal. Seperti yang Anda catat, msm::tnormdidasarkan pada pendekatan Robert, sementara truncnorm::truncnormmengimplementasikan sampler accept-reject Geweke (1991); ini terkait dengan pendekatan dalam makalah Robert. Catatan yang msm::tnormmencakup fungsi densitas, cdf, dan kuantil (invers cdf) dengan Rcara biasa .

Referensi yang lebih tua dengan pendekatan adalah buku Luc Devroye ; sejak keluar dari cetak, dia mendapatkan kembali hak cipta dan membuatnya tersedia sebagai unduhan.

Contoh khusus Anda sama dengan sampling standar normal terpotong pada 1 (jika $t$ adalah titik pemotongan, $(t-\mu)/\sigma = (5-3)/2 = 1$), dan kemudian scaling hasilnya (kalikan dengan $\sigma$ dan tambahkan $\mu$).

Dalam kasus khusus itu, Robert menyarankan bahwa ide Anda (dalam inkarnasi kedua atau ketiga) cukup masuk akal. Anda mendapatkan nilai yang dapat diterima sekitar 84% dari waktu dan menghasilkan sekitar$1.19 n$ normals rata-rata (Anda dapat bekerja di luar batas sehingga Anda menghasilkan nilai yang cukup menggunakan algoritma vektor mengatakan 99,5% dari waktu, dan kemudian sesekali menghasilkan beberapa yang terakhir kurang efisien - bahkan satu per satu).

Ada juga diskusi tentang implementasi dalam kode R di sini (dan di Rccp di jawaban lain untuk pertanyaan yang sama, tetapi kode R sebenarnya ada lebih cepat). Kode R biasa di sana menghasilkan 50.000 normal terpotong dalam 6 milidetik, meskipun normal terpotong tertentu hanya memotong ekor ekstrim, jadi pemotongan yang lebih substantif akan berarti hasilnya lebih lambat. Ini mengimplementasikan ide menghasilkan "terlalu banyak" dengan menghitung berapa banyak yang harus dihasilkan untuk hampir pasti mendapatkan cukup.

Jika saya perlu hanya satu jenis tertentu yang terpotong normal berkali-kali, saya mungkin akan melihat mengadaptasi versi metode ziggurat, atau yang serupa, dengan masalah.

Sebenarnya sepertinya Nicolas Chopin sudah melakukannya, jadi saya bukan satu-satunya orang yang terpikir:

http://arxiv.org/abs/1201.6140

Dia membahas beberapa algoritma lain dan membandingkan waktu untuk 3 versi algoritmanya dengan algoritma lain untuk menghasilkan 10 ^ 8 normals acak untuk berbagai titik pemotongan.

Mungkin tidak mengejutkan, algoritmanya ternyata relatif cepat.

Dari grafik di kertas, bahkan yang paling lambat dari algoritma yang dia bandingkan dengan di (untuk mereka) titik pemotongan terburuk menghasilkan $10^8$ nilai dalam waktu sekitar 3 detik - yang menunjukkan bahwa salah satu algoritma yang dibahas di sana mungkin dapat diterima jika diterapkan dengan cukup baik.

Sunting: Salah satu yang saya tidak yakin disebutkan di sini (tapi mungkin itu di salah satu tautan) adalah untuk mengubah (melalui inverse normal cdf) seragam terpotong - tetapi seragam dapat dipotong dengan hanya menghasilkan seragam dalam batas pemotongan . Jika invers cdf normal cepat, ini cepat dan mudah dan berfungsi dengan baik untuk berbagai titik pemotongan.

Sebagai lanjutan dari referensi @ glen_b dan memfokuskan secara eksklusif pada implementasi R. Ada beberapa fungsi yang dirancang untuk mengambil sampel dari distribusi normal terpotong:

• rtruncnorm(100, a=-Inf, b=5, mean=3, sd=2)dalam truncnormpaket
• rtnorm(100, 3, 2, upper=5)dalam msmpaket