Bagaimana sampel dari untuk variabel acak, masing-masing dengan fungsi massa yang berbeda, dalam R?

Dalam R, saya memiliki matriks di mana 'th deretan bersesuaian dengan distribusi pada . Pada dasarnya, saya perlu sampel dari setiap baris secara efisien. Implementasi naif adalah:$N \times K$$P$$i$$P$$\{1, ..., K\}$

X = rep(0, N);
for(i in 1:N){
    X[i] = sample(1:K, 1, prob = P[i, ]);
}

Ini terlalu lambat. Pada prinsipnya saya bisa memindahkan ini ke C tapi saya yakin pasti ada cara yang ada untuk melakukan ini. Saya ingin sesuatu dalam semangat kode berikut (yang tidak berfungsi):

X = sample(1:K, N, replace = TRUE, prob = P)

EDIT: Untuk motivasi, ambil dan . Saya punya matriks semua dan saya perlu sampel vektor dari masing-masing.$N = 10000$$K = 100$$P_1, ..., P_{5000}$$N \times K$

Kita dapat melakukan ini dalam beberapa cara sederhana . Yang pertama adalah kode mudah, mudah dimengerti dan cukup cepat. Yang kedua sedikit lebih rumit, tetapi jauh lebih efisien untuk ukuran masalah ini daripada metode pertama atau pendekatan lain yang disebutkan di sini.

Metode 1 : Cepat dan kotor.

Untuk mendapatkan pengamatan tunggal dari distribusi probabilitas setiap baris, kita cukup melakukan hal berikut.

# Q is the cumulative distribution of each row.
Q <- t(apply(P,1,cumsum))

# Get a sample with one observation from the distribution of each row.
X <- rowSums(runif(N) > Q) + 1

Ini menghasilkan distribusi kumulatif dari setiap baris dan kemudian sampel satu pengamatan dari setiap distribusi. Perhatikan bahwa jika kita dapat menggunakan kembali maka kita dapat menghitung sekali dan menyimpannya untuk digunakan nanti. Namun, pertanyaannya membutuhkan sesuatu yang berfungsi untuk berbeda di setiap iterasi.$P$ $P$$Q$$P$

Jika Anda membutuhkan beberapa ( ) pengamatan dari setiap baris, maka ganti baris terakhir dengan yang berikut.$n$

# Returns an N x n matrix
X <- replicate(n, rowSums(runif(N) > Q)+1)

Ini sebenarnya bukan cara yang sangat efisien secara umum untuk melakukan ini, tetapi memang memanfaatkan Rkemampuan vektorisasi, yang biasanya merupakan penentu utama kecepatan eksekusi. Juga mudah dipahami.

Metode 2 : Menggabungkan cdfs.

Misalkan kita memiliki fungsi yang mengambil dua vektor, yang kedua diurutkan dalam urutan nondecreasing monoton dan menemukan indeks dalam vektor kedua dari batas bawah terbesar dari setiap elemen di yang pertama. Kemudian, kita bisa menggunakan fungsi ini dan trik yang licin: Cukup buat jumlah kumulatif dari semua baris cdf. Ini memberikan vektor yang meningkat secara monoton dengan elemen dalam kisaran .$[0,N]$

Ini kodenya.

i <- 0:(N-1)

# Cumulative function of the cdfs of each row of P.
Q <- cumsum(t(P))

# Find the interval and then back adjust
findInterval(runif(N)+i, Q)-i*K+1

Perhatikan apa yang dilakukan baris terakhir, itu menciptakan variabel acak yang didistribusikan dalam dan kemudian memanggil untuk menemukan indeks batas bawah terbesar dari setiap entri. . Jadi, ini memberitahu kita bahwa elemen pertama akan ditemukan antara indeks 1 dan indeks , yang kedua akan ditemukan antara indeks dan , dll, masing-masing sesuai dengan distribusi baris sesuai . Maka kita perlu kembali mentransformasikan untuk mendapatkan masing-masing indeks kembali dalam rentang .$(0,1), (1,2), \dots, (N-1,N)$findIntervalrunif(N)+i$K$$K+1$$2K$$P$$\{1,\ldots,K\}$

Karena findIntervalcepat baik secara algoritmik maupun dari segi implementasi, metode ini ternyata sangat efisien.

Sebuah tolok ukur

Di laptop lama saya (MacBook Pro, 2,66 GHz, 8GB RAM), saya mencoba ini dengan dan dan menghasilkan 5000 sampel ukuran , persis seperti yang disarankan dalam pertanyaan yang diperbarui, untuk total 50 juta varian acak .$N = 10000$$K = 100$$N$

Kode untuk Metode 1 membutuhkan waktu hampir 15 menit untuk dijalankan, atau sekitar 55 ribu variasi acak per detik. Kode untuk Metode 2 membutuhkan waktu sekitar empat setengah menit untuk dijalankan, atau sekitar 183 ribu variasi acak per detik.

Berikut adalah kode demi reproduktifitas. (Perhatikan bahwa, seperti yang ditunjukkan dalam komentar, dihitung ulang untuk masing-masing dari 5000 iterasi untuk mensimulasikan situasi OP.)$Q$

# Benchmark code
N <- 10000
K <- 100

set.seed(17)
P <- matrix(runif(N*K),N,K)
P <- P / rowSums(P)

method.one <- function(P)
{
    Q <- t(apply(P,1,cumsum))
    X <- rowSums(runif(nrow(P)) > Q) + 1
}

method.two <- function(P)
{
    n <- nrow(P)
    i <- 0:(n-1)
    Q <- cumsum(t(P))
    findInterval(runif(n)+i, Q)-i*ncol(P)+1
}

Ini outputnya.

# Method 1: Timing
> system.time(replicate(5e3, method.one(P)))
   user  system elapsed 
691.693 195.812 899.246 

# Method 2: Timing
> system.time(replicate(5e3, method.two(P)))
   user  system elapsed 
182.325  82.430 273.021 

Postscript : Dengan melihat kode untuk findInterval, kita dapat melihat bahwa ia melakukan beberapa pemeriksaan pada input untuk melihat apakah ada NAentri atau jika argumen kedua tidak diurutkan. Karenanya, jika kami ingin memeras lebih banyak kinerja dari ini, kami dapat membuat versi modifikasi kami sendiri findIntervalyang menghapus cek ini yang tidak perlu dalam kasus kami.

Sebuah forlingkaran mungkin sangat lambat di R. Bagaimana dengan vektorisasi sederhana ini sapply?

n <- 10000
k <- 200

S <- 1:k
p <- matrix(rep(1 / k, n * k), nrow = n, ncol = k)
x <- numeric(n)

x <- sapply(1:n, function(i) sample(S, 1, prob = p[i,]))

Tentu saja, p seragam ini hanya untuk pengujian.