Bagaimana cara menghasilkan angka berdasarkan distribusi diskrit arbitrer?

28

Bagaimana cara menghasilkan angka berdasarkan distribusi diskret arbitrer?

Misalnya, saya memiliki satu set angka yang ingin saya hasilkan. Katakanlah mereka diberi label dari 1-3 sebagai berikut.

1: 4%, 2: 50%, 3: 46%

Pada dasarnya, persentase adalah probabilitas bahwa mereka akan muncul di output dari generator nomor acak. Saya memiliki generator nomor acak yang akan menghasilkan distribusi seragam dalam interval [0, 1]. Apakah ada cara untuk melakukan ini?

Tidak ada batasan berapa banyak elemen yang bisa saya miliki, tetapi% akan bertambah hingga 100%.

FurtiveFelon
sumber
2
Saya mungkin menyarankan untuk menentukan "... distribusi diskret arbitrer" dalam judul, jika itu adalah pertanyaan Anda. Kasus kontinu berbeda.
David M Kaplan
3
Cara umum adalah dengan melakukan pencarian biner dalam daftar probabilitas kumulatif, yang dalam contoh ini adalah (0,0.04,0.54,1.0) . Rata-rata ini membutuhkan log(n)/2 probe per peristiwa generasi. Jika tidak ada probabilitas yang sangat kecil, Anda bisa mendapatkan kinerja O(1) dengan membuat vektor nilai-nilai dengan spasi yang sama di [0,1] dan (dalam tahap prakomputasi) yang menetapkan hasil untuk setiap nilai. Misalnya, dalam contoh ini Anda dapat membuat vektor (1,1,1,1,2,,2,3,,3) (dengan50 2 dan46 3). Hasilkan seragam, kalikan dengan 100, dan indeks ke dalam vektor ini: selesai.
whuber
Lihat juga di sini
Glen_b -Reinstate Monica
Tautan "di sini" sebenarnya menautkan ke pertanyaan ini, @Glen_b ... kesalahan salin-tempel?
buruzaemon
@buruzaemon terima kasih ya itu kesalahan; Saya sudah memperbaikinya.
Glen_b -Reinstate Monica

Jawaban:

26

Salah satu algoritma terbaik untuk pengambilan sampel dari distribusi diskrit adalah metode alias .

Metode alias (efisien) mengkompilasi struktur data dua dimensi untuk mempartisi persegi panjang ke dalam area yang proporsional dengan probabilitas.

Figure

Dalam skema ini dari situs direferensikan, persegi panjang dengan tinggi Unit telah dipartisi menjadi empat macam daerah - sebagai dibedakan oleh warna - dalam proporsi , 1 / 3 , 1 / 12 , dan 1 / 12 , di memesan sampel berulang kali dari distribusi diskrit dengan probabilitas ini. Strip vertikal memiliki lebar (unit) konstan. Masing-masing dibagi menjadi satu atau dua potong. Identitas potongan dan lokasi divisi vertikal disimpan dalam tabel yang dapat diakses melalui indeks kolom.1/21/31/121/12

Tabel dapat disampel dalam dua langkah sederhana (satu untuk setiap koordinat) yang membutuhkan hanya menghasilkan dua nilai seragam independen dan perhitungan . Ini meningkatkan pada perhitungan O ( log ( n ) ) yang diperlukan untuk membalikkan CDF diskrit seperti yang dijelaskan dalam balasan lain di sini.O(1)O(log(n))

Lucas
sumber
2
Algoritma ini hanya terbaik jika probabilitasnya murah untuk dihitung. Sebagai contoh jika besar, mungkin lebih baik untuk tidak membangun seluruh pohon. n
probabilityislogic
3
+1 Sejauh ini ini adalah satu - satunya jawaban yang menyarankan dan menjelaskan algoritma yang efisien.
whuber
19

Anda dapat melakukan ini dengan mudah di R, cukup tentukan ukuran yang Anda butuhkan:

sample(x=c(1,2,3), size=1000, replace=TRUE, prob=c(.04,.50,.46))
Dominic Comtois
sumber
3
Secara pribadi, saya lebih suka algoritma (atau suatu tempat untuk mempelajari pengetahuan yang diperlukan), karena saya mencoba untuk memasukkan ini ke dalam aplikasi yang saya bangun :) Terima kasih banyak atas jawaban Anda :)
FurtiveFelon
Hmmm ok ... Mengetahui lebih banyak tentang apa yang ingin Anda lakukan akan membantu kami membimbing Anda. Bisakah Anda memberi tahu kami lebih banyak tentang itu? (Tujuan, konteks, dll.)
Dominic Comtois
Ini untuk pemungutan suara. Misalnya, saya memiliki banyak foto, dan saya hanya dapat menampilkan 6 kepada pengguna pada satu waktu, saya ingin memasukkan "terbaik" ke pengguna pada satu waktu, dan pengguna dapat memilih naik atau turun pada setiap foto . Solusi paling sederhana yang dapat berfungsi saat ini adalah skema yang saya uraikan (setiap angka mewakili foto, setiap suara akan menurunkan probabilitas pada foto itu, dan meningkat pada yang lainnya)
FurtiveFelon
1
@furtivefelon, Anda selalu dapat port kode dari R, o mencari tahu algoritma dari kode dan mengimplementasikannya.
mpiktas
Saya pikir Anda mungkin mendapatkan saran yang bagus (lebih baik) tentang Stack Overflow, karena mungkin ada beberapa solusi terkenal untuk tujuan khusus ini. Saya sarankan juga menyertakan info dari komentar terakhir Anda langsung ke pertanyaan Anda.
Dominic Comtois
19

Dalam contoh Anda, katakan Anda menggambar nilai Seragam pseudorandom Anda [0,1] dan menyebutnya U. Kemudian keluaran:

1 jika U <0,04

2 jika U> = 0,04 dan U <0,54

3 jika U> = 0,54

Jika% yang ditentukan adalah a, b, ..., cukup keluaran

nilai 1 jika U

nilai 2 jika U> = a dan U <(a + b)

dll.

Pada dasarnya, kami memetakan% ke dalam himpunan bagian dari [0,1], dan kami tahu probabilitas bahwa nilai acak seragam jatuh ke dalam rentang apa pun hanya panjang rentang itu. Menempatkan rentang dalam urutan tampaknya cara paling sederhana, jika tidak unik, untuk melakukannya. Ini dengan asumsi bahwa Anda hanya bertanya tentang distribusi diskrit; untuk yang berkelanjutan, dapat melakukan sesuatu seperti "sampel penolakan" ( entri Wikipedia ).

David M Kaplan
sumber
8
Algoritma lebih cepat jika Anda mengurutkan kategori dalam urutan penurunan probabilitas. Dengan begitu, Anda melakukan lebih sedikit tes (rata-rata) per nomor acak yang dihasilkan.
jbowman
1
Hanya dengan menambahkan catatan singkat tentang penyortiran - ini akan efektif hanya jika Anda melakukannya sekali pada awal skema pengambilan sampel - sehingga tidak akan berfungsi dengan baik untuk kasus-kasus di mana probabilitas itu sendiri dijadikan sampel sebagai bagian dari keseluruhan skema yang lebih besar ( mis. dan kemudian P r ( Y = j ) = p j ). Dengan menyortir dalam hal ini Anda menambahkan operasi penyortiran ke dalam setiap iterasi pengambilan sampel - yang akan menambahkan O ( n log ( n ) )pjDistPr(Y=j)=pjO(nlog(n))waktu untuk setiap iterasi. Namun, mungkin berguna untuk mengurutkan berdasarkan perkiraan pada ukuran probabilitas pada awal dalam kasus ini.
probabilityislogic
4

Misalkan ada hasil yang mungkin diskrit. Anda membagi interval [ 0 , 1 ] ke dalam sub-terminal berdasarkan pada fungsi massa probabilitas kumulatif, F , untuk memberikan interval yang dipartisi ( 0 , 1 )m[0,1]F(0,1)

I1I2Im

di mana dan F ( 0 ) 0 . Dalam contoh Anda m = 3 danIj=(F(j1),F(j))F(0)0m=3

I1=(0,.04),     I2=(.04,.54),     I3=(.54,1)

karena dan F ( 2 ) = .54 dan F ( 3 ) = 1 .F(1)=.04F(2)=.54F(3)=1

Kemudian Anda dapat menghasilkan dengan distribusi F menggunakan algoritma berikut:XF

(1) menghasilkan UUniform(0,1)

(2) Jika , maka X = j .UIjX=j

  • Langkah ini dapat dilakukan dengan melihat apakah kurang dari masing-masing probabilitas kumulatif , dan melihat di mana titik perubahan (dari menjadi ) terjadi, yang seharusnya merupakan masalah menggunakan operator boolean dalam bahasa pemrograman apa pun yang Anda gunakan dan menemukan di mana yang pertama terjadi dalam vektor.UTRUEFALSEFALSE

Perhatikan bahwa akan berada tepat pada salah satu interval I j karena mereka terpisah dan partisi [ 0 , 1 ] .UIj[0,1]

Makro
sumber
Bukankah seharusnya semua interval itu setengah tertutup? Kalau tidak, batas antara interval tidak termasuk .. yaitu. {[0,0.04), [0.04,0.54), [0.54,1]}
naught101
1
untuk setiap titik u (yaitu ukuran Lebesgue dari interval setengah terbuka adalah sama dengan yang dari interval terbuka) jadi saya tidak berpikir itu penting. P(U=u)=0u
Makro
1
Pada mesin digital presisi-terbatas, mungkin suatu hari nanti sebelum akhir jagat raya ...
jbowman
1
Cukup adil, @whuber, lihat edit saya.
Makro
1
OK, itu adalah sebuah algoritma. BTW, mengapa Anda tidak mengembalikan sesuatu seperti min(which(u < cp))? Sebaiknya hindari menghitung ulang jumlah kumulatif pada setiap panggilan juga. Dengan perhitungan itu, seluruh algoritma dikurangi menjadi min(which(runif(1) < cp)). Atau lebih baik, karena OP meminta untuk menghasilkan angka ( jamak ), vektorkan sebagai n<-10; apply(matrix(runif(n),1), 2, function(u) min(which(u < cp))).
whuber
2

Salah satu algoritma sederhana adalah mulai dengan nomor acak seragam Anda dan dalam satu lingkaran pertama kurangi probabilitas pertama, jika hasilnya negatif maka Anda mengembalikan nilai pertama, jika masih positif maka Anda pergi ke iterasi berikutnya dan kurangi probabilitas berikutnya , periksa apakah negatif, dll.

Ini bagus karena jumlah nilai / probabilitas bisa tak terbatas tetapi Anda hanya perlu menghitung probabilitas ketika Anda mendekati angka-angka itu (untuk sesuatu seperti menghasilkan dari Poisson atau distribusi binomial negatif).

Jika Anda memiliki serangkaian probabilitas terbatas, tetapi akan menghasilkan banyak angka dari mereka maka akan lebih efisien untuk mengurutkan probabilitas sehingga Anda mengurangi yang pertama, kemudian yang ke-2 berikutnya dan seterusnya.

Greg Snow
sumber
2

First of all, let me draw your attention to a python library with ready-to-use classes for either integer or floating point random number generation that follow arbitrary distribution.

Generally speaking there are several approaches to this problem. Some are linear in time, but require large memory storage, some run in O(n log(n)) time. Some are optimized for integer numbers and some are defined for circular histograms (for example: generating random time spots during a day). In the above mentioned library I used this paper for integer number cases and this recipe for floating point numbers. It (still) lacks circular histogram support and is generally messy, but it works well.

Boris Gorelik
sumber
2

I had the same problem. Given a set where each item has a probability and whose items' probabilities sum up to one, I wanted to draw a sample efficiently, i.e. without sorting anything and without repeatedly iterating over the set.

The following function draws the lowest of N uniformly distributed random numbers within the interval [a,1). Let r be a random number from [0,1).

next(N,a)=1(1a)rN

You can use this function to draw an ascending series (ai) of N uniformly distributed random numbers in [0,1). Here is an example with N=10:

a0=next(10,0)
a1=next(9,a0)
a2=next(8,a1)

a9=next(1,a8)

While drawing that ascending series (ai) of uniformly distributed numbers, iterate over the set of probabilities P which represents your arbitraty (yet finite) distribution. Let 0k<|P| be the iterator and pkP. After drawing ai, increment k zero or more times until p0pk>ai. Then add pk to your sample and move on with drawing ai+1.


Example with the op's set {(1,0.04),(2,0.5),(3,0.46)} and sample size N=10:

i  a_i    k  Sum   Draw
0  0.031  0  0.04  1
1  0.200  1  0.54  2
2  0.236  1  0.54  2
3  0.402  1  0.54  2
4  0.488  1  0.54  2
5  0.589  2  1.0   3
6  0.625  2  1.0   3
7  0.638  2  1.0   3
8  0.738  2  1.0   3
9  0.942  2  1.0   3

Sample: (1,2,2,2,2,3,3,3,3,3)


If you wonder about the next function: It is the inverse of the probability that one of N uniformly distributed random numbers lies within the interval [a,x) with x1.

casi
sumber
It appears the problem you are addressing abruptly changed in the second paragraph from one that samples from an arbitrary discrete distribution to sampling from a uniform distribution. Its solution appears not to be relevant to the question that was asked here.
whuber
I clarified the last part.
casi
Your answer still seems unrelated to the question. Could you perhaps provide a small but nontrivial worked example of your algorithm? Show us how it would generate a single draw from the set {1,2,3} according to the probabilities given in the question.
whuber
I added an example. My answer has something in common with David M Kaplan's answer (stats.stackexchange.com/a/26860/93386), but requires just one instead of N (= sample size) iterations over the set, at the expense of drawing N N-th roots. I profiled both procedures, and mine was much faster.
casi
Thank you for the clarification (+1). It may be of interest to many readers that this isn't a simple random sample, because the outcomes appear in a predetermined, fixed order: a random permutation would have to be applied to the results in order to create a simple random sample. You might also be interested in a parallelizable version of this algorithm in which
aj=i=1jlog(ui)i=1N+1log(ui)
where u1,,uN+1 is a simple random sample of Uniform(0,1] variates.
whuber