Bagaimana cara menghasilkan sampel secara seragam secara acak dari berbagai variabel diskrit yang mengalami kendala?

Saya ingin membuat proses Monte Carlo untuk mengisi guci dengan N bola warna I, C [i]. Setiap warna C [i] memiliki jumlah bola minimum dan maksimum yang harus ditempatkan di dalam guci.

Misalnya, saya mencoba memasukkan 100 bola ke dalam guci, dan dapat mengisinya dengan empat warna:

• Merah - Minimum 0, Maksimum 100 # NB, maksimum yang sebenarnya tidak dapat direalisasikan.
• Biru - Minimum 50, Maksimum 100
• Kuning - Minimum 0, Maksimum 50
• Hijau - Minimum 25, Maksimal 75

Bagaimana saya bisa menghasilkan sampel N yang dipastikan akan didistribusikan secara merata di seluruh hasil yang mungkin?

Saya telah melihat solusi untuk masalah ini di mana bola tidak memiliki minimum atau maksimum, atau sebagai alternatif, memiliki minima dan maxima implisit yang sama. Lihat misalnya, diskusi ini tentang topik yang sedikit berbeda:

Hasilkan bobot yang terdistribusi secara seragam yang menghasilkan kesatuan?

Tetapi saya mengalami masalah dalam menggeneralisasikan solusi ini.

Biarkan menunjukkan jumlah bola warna . Juga, biarkan dan menyatakan jumlah minimum dan maksimum bola warna . Kami ingin mengambil sampel secara seragam pada subjek acak dengan batasan berikut:$n_i$$C_i$$m_i$$M_i$$C_i$$(n_1, \dots, n_I)$

1. $m_i \leq n_i \leq M_i$
2. $\sum_{i=1}^I n_i = N$

Pertama-tama, Anda dapat menghapus batasan batas bawah (yaitu ) dengan memilih bola warna pada awalnya. Ini memodifikasi dua kendala sebagai berikut:$m_i \leq n_i$$m_i$$C_i$

1. $0 \leq n_i \leq b_i = M_i - m_i$
2. $\sum_{i=1}^I n_i = B = N - \sum_{i=1}^I m_i$

Misalkan menunjukkan distribusi seragam yang kami minati. Kita dapat menggunakan aturan rantai dan pemrograman dinamis untuk mengambil sampel dari efisien. Pertama, kami menggunakan aturan rantai untuk menulis sebagai berikut: mana adalah distribusi marjinal . Perhatikan bahwa$P(n_1, \dots, n_I \mid B, b_{1:I})$$P$$P$

Berikut ini adalah implementasi ide Matlab. Kompleksitas algoritma adalah mana . Kode menggunakan dihasilkan secara acak di setiap proses. Akibatnya, beberapa kasus uji yang dihasilkan mungkin tidak memiliki solusi yang valid, dalam hal ini kode mencetak pesan peringatan.$O(I \times B \times K)$$K = \max_{i=1}^I b_i$$m_i$

global dpm b

I = 5; % number of colors
N = 300; % total number of balls

m = randi(50, 1, I)-1; % minimum number of balls from each from each color
M = 99*ones(1, I); % maximum number of balls from each color

% print original constraints
print_constraints(I, N, m, M, 'original constraints');

% remove the lower bound constraints
b = M - m;
B = N - sum(m);
m = zeros(size(m));

% print transformed constraints
print_constraints(I, B, zeros(1, I), b, 'transformed constraints');

% initialize the dynamic programming matrix (dpm)
% if dpm(i, n) <> -1, it denotes the value of the following marginal probability
% \sum_{k=1}^{i-1} P(n_1, ..., n_i |
dpm = -ones(I, B);

% sample the number of balls of each color, one at a time, using chain rule
running_B = B;  % we change the value of "running_B" on the fly, as we sample balls of different colors
for i = I : -1 : 1
    % compute marginal distribution P(n_i)

    % instead of P(n_i) we compute q(n_i) which is un-normalized.
    q_ni = zeros(1, b(i) + 1); % possibilities for ni are 0, 1, ..., b(i)
    for ni = 0 : b(i)
        q_ni(ni+1) = dpfunc(i-1, running_B-ni);
    end
    if(sum(q_ni) == 0)
        fprintf('Impossible!!! constraints can not be satisfied!\n');
        return;
    end
    P_ni = q_ni / sum(q_ni);
    ni = discretesample(P_ni, 1) - 1;
    fprintf('n_%d=%d\n', i, ni);
    running_B = running_B - ni;
end

di mana fungsi print_constraints adalah

function [] = print_constraints(I, N, m, M, note)
    fprintf('\n------ %s ------ \n', note);
    fprintf('%d <= n_%d <= %d\n', [m; [1:I]; M]);
    fprintf('========================\n');
    fprintf('sum_{i=1}^%d n_i = %d\n', I, N);
end

dan fungsi dpfunc melakukan perhitungan pemrograman dinamis sebagai berikut:

function res = dpfunc(i, n)
    global dpm b

    % check boundary cases
    if(n == 0)
        res = 1;
        return;
    end
    if(i == 0) % gets here only if n <> 0
        res = 0;
        return;
    end

    if(n < 0)
        res = 0;
        return;
    end

    if(dpm(i, n) == -1) % if <> -1, it has been compute before, so, just use it!
        % compute the value of dpm(i, n) = \sum_{n_1, ..., n_i} valid(n, n_1, ..., n_i)
        % where "valid" return 1 if \sum_{j=1}^i n_i = n and 0 <= n_i <= b_i, for all i
        % and 0 otherwise.
        dpm(i, n) = 0;
        for ni = 0 : b(i)
            dpm(i, n) = dpm(i, n) + dpfunc(i-1, n-ni);
        end
    end
    res = dpm(i, n);
end

dan akhirnya, fungsi discretesample (p, 1) mengambil sampel acak dari distribusi diskrit . Anda dapat menemukan satu implementasi dari fungsi ini di sini .$p$

Mari kita pertimbangkan generalisasi masalah ini. Ada kaleng cat warna berbeda dan bola. Bisakah menyimpan hingga bola. Anda ingin membuat konfigurasi bola di kaleng yang memiliki setidaknya bola di can untuk setiap , setiap konfigurasi dengan probabilitas yang sama.$m=4$$m=4$$n^{(0)}=100$$i$$a^{(0)}_i = (100, 100, 50, 75)$$b_i = (0, 50, 0, 25)$$i$$i$

Konfigurasi seperti ini berada dalam korespondensi satu-ke-satu dengan konfigurasi yang diperoleh setelah mengeluarkan bola dari can , membatasi bola tersisa paling banyak per kaleng. Karena itu saya hanya akan menghasilkan ini dan membiarkan Anda menyesuaikannya setelah itu (dengan meletakkan bola kembali ke can untuk setiap ).$b_i$$i$$n = n^{(0)} - \sum_i b_i = 100 - (0+50+0+25)=25$$a_i = a^{(0)}_i - b_i = (100, 50, 50, 50)$$b_i$$i$$i$

Untuk menghitung konfigurasi ini, perbaiki semua kecuali dua indeks, misalkan dan . Misalkan ada bola sudah di dapat untuk setiap berbeda dari dan . Itu menyisakan bola . Bersyarat di tempat bola bola ini didistribusikan secara seragam di dalam kaleng dan . Kemungkinan konfigurasi yang di nomor (lihat komentar), mulai dari menempatkan sebanyak bola di dapat$i$$j$$s_k$$k$$k$$i$$j$$s_i+s_j$$n - (s_i+s_j)$$i$$j$$1 + \min(a_i + a_j - s_i - s_j, s_i+s_j)$$i$mungkin semua jalan melalui menempatkan sebanyak bola di dapat mungkin.$j$

Jika mau, Anda bisa menghitung jumlah total konfigurasi dengan menerapkan argumen ini secara rekursif ke kaleng tersisa . Namun, untuk mendapatkan sampel kita bahkan tidak perlu tahu jumlah ini. Yang perlu kita lakukan adalah berulang kali mengunjungi semua pasangan kaleng yang tidak berurutan kaleng dan secara acak (dan seragam) mengubah distribusi bola di dalam kedua kaleng itu. Ini adalah rantai Markov dengan distribusi probabilitas terbatas yang seragam pada semua kemungkinan kondisi (seperti yang ditunjukkan dengan menggunakan metode standar). Oleh karena itu sudah cukup untuk memulai di setiap$m-2$$\{i,j\}$nyatakan, jalankan rantai cukup lama untuk mencapai distribusi membatasi, dan kemudian melacak keadaan yang dikunjungi oleh prosedur ini. Seperti biasa, untuk menghindari korelasi serial, urutan keadaan ini harus "ditipiskan" dengan melewatinya (atau ditinjau kembali secara acak). Menipis dengan faktor sekitar setengah jumlah kaleng cenderung bekerja dengan baik, karena setelah itu banyak langkah rata-rata setiap kaleng dapat terpengaruh, menghasilkan konfigurasi yang benar-benar baru.

Algoritma ini membutuhkan biaya untuk menghasilkan setiap konfigurasi acak rata-rata. Meskipun ada algoritma lain, ini memiliki keunggulan karena tidak perlu melakukan perhitungan kombinatorial sebelumnya.$O(m)$$O(m)$

Sebagai contoh, mari kita selesaikan situasi yang lebih kecil secara manual. Misalnya dan . Ada 15 konfigurasi yang valid, yang dapat ditulis sebagai string nomor hunian. Misalnya, tempatkan dua bola ke kaleng kedua dan satu bola di kaleng keempat. Meniru argumen, mari kita pertimbangkan total hunian dari dua kaleng pertama. Ketika itu adalah bola, tidak ada bola yang tersisa untuk dua kaleng terakhir. Itu memberi negara$a=(4,3,2,1)$$n=3$0201$s_1+s_2=3$

30**, 21**, 12**, 03**

di mana **mewakili semua nomor hunian yang mungkin untuk dua kaleng terakhir: yaitu 00,. Ketika , adalah$s_1+s_2=2$

20**, 11**, 02**

di mana sekarang **dapat berupa 10atau 01. Itu memberi lebih banyak status. Ketika , adalah$3\times 2=6$$s_1+s_2=1$

10**, 01**

di mana sekarang **bisa 20, 11tetapi tidak 02 (karena batas dari satu bola di dapat terakhir). Itu memberi lebih banyak status. Akhirnya, ketika , semua bola berada di dua kaleng terakhir, yang harus penuh hingga batas dan . The kemungkinan sama oleh karena itu negara-negara yang$2\times 2=4$$s_1+s_2=0$$2$$1$$4+6+4+1=15$

3000, 2100, 1200, 0300; 2010, 2001, 1110, 1101, 0210, 0201; 1020, 1011, 0120, 0111; 0021.

Dengan menggunakan kode di bawah ini, urutan konfigurasi seperti itu dihasilkan dan ditipiskan ke setiap sepertiga, menciptakan konfigurasi dari negara. Frekuensi mereka adalah sebagai berikut:$10,009$$3337$$15$

State: 3000 2100 1200 0300 2010 1110 0210 1020 0120 2001 1101 0201 1011 0111 0021 
Count:  202  227  232  218  216  208  238  227  237  209  239  222  243  211  208 

Sebuah uji keseragaman memberikan nilai , ( derajat kebebasan): bahwa kesepakatan yang indah dengan hipotesis bahwa prosedur ini menghasilkan negara kemungkinan sama.$\chi^2$$\chi^2$$11.2$$p=0.67$$14$

Ini Rkode diatur untuk menangani situasi dalam pertanyaan. Ubah adan nuntuk bekerja dengan situasi lain. Atur Nagar cukup besar untuk menghasilkan jumlah realisasi yang Anda butuhkan setelah penjarangan .

Kode ini menipu sedikit dengan bersepeda secara sistematis melalui semua pasangan. Jika Anda ingin menjadi ketat tentang menjalankan rantai Markov, menghasilkan , dan secara acak, seperti yang diberikan dalam kode berkomentar. $(i,j)$ijij

#
# Gibbs-like sampler.
#
# `a` is an array of maximum numbers of balls of each type.  Its values should
#     all be integers greater than zero.
# `n` is the total number of balls.
#------------------------------------------------------------------------------#
g <- function(j, state, a) {
  #
  # `state` contains the occupancy numbers.
  # `a`     is the array of maximum occupancy numbers.
  # `j`     is a pair of indexes into `a` to "rotate".
  #
  k <- sum(state[j]) # Total occupancy.
  x <- floor(runif(1, max(0, k - a[j[2]]), min(k, a[j[1]]) + 1))
  state[j] <- c(x, k-x)
  return(state)
}
#
# Set up the problem.
#
a <- c(100, 50, 50, 50)
n <- 25

# a <- 4:1
# n <- 3
#
# Initialize the state.
#
state <- round(n * a / sum(a))
i <- 1
while (sum(state) < n) {
  if (state[i] < a[i]) state[i] <- state[i] + 1
  i <- i+1
}
while (sum(state) > n) {
  i <- i-1
  if (state[i] > 0) state[i] <- state[i] - 1
}
#
# Conduct a sequence of random changes.
#
set.seed(17)
N <- 1e5
sim <- matrix(state, ncol=1)
u <- ceiling(N / choose(length(state), 2) / 2)
i <- rep(rep(1:length(state), each=length(state)-1), u)
j <- rep(rep(length(state):1, length(state)-1), u)
ij <- rbind(i, j)
#
# Alternatively, generate `ij` randomly:
#   i <- sample.int(length(state), N, replace=TRUE)
#   j <- sample.int(length(state)-1, N, replace=TRUE)
#   ij <- rbind(i, ((i+j-1) %% length(state))+1)
#
sim <- cbind(sim, apply(ij, 2, function(j) {state <<- g(j, state, a); state}))
rownames(sim) <- paste("Can", 1:nrow(sim))
#
# Thin them for use.  Each column is a state.
#
thin <- function(x, stride=1, start=1) {
  i <- round(seq(start, ncol(x), by=stride))
  x[, i]
}
#
# Make a scatterplot of the results, to illustrate.
#
par(mfrow=c(1,1))
s <- thin(sim, stride=max(1, N/1e4))
pairs(t(s) + runif(length(s), -1/2, 1/2), cex=1/2, col="#00000005", pch=16)