Wie werden Stichproben gleichmäßig zufällig aus mehreren diskreten Variablen generiert, die Einschränkungen unterliegen?

Ich möchte einen Monte-Carlo-Prozess generieren, um eine Urne mit N Kugeln der Farbe I, C [i], zu füllen. Jede Farbe C [i] hat eine minimale und maximale Anzahl von Kugeln, die in die Urne gelegt werden sollten.

Zum Beispiel versuche ich, 100 Kugeln in die Urne zu legen und kann sie mit vier Farben füllen:

• Rot - Minimum 0, Maximum 100 # NB, das tatsächliche Maximum kann nicht realisiert werden.
• Blau - Minimum 50, Maximum 100
• Gelb - Minimum 0, Maximum 50
• Grün - Minimum 25, Maximum 75

Wie kann ich N Stichproben generieren, die eine gleichmäßige Verteilung auf mögliche Ergebnisse gewährleisten?

Ich habe Lösungen für dieses Problem gesehen, bei denen die Kugeln kein Minimum oder Maximum haben oder alternativ die gleichen impliziten Minima und Maxima haben. Siehe zum Beispiel diese Diskussion zu einem etwas anderen Thema:

Gleichmäßig verteilte Gewichte erzeugen, die sich zu einer Einheit summieren?

Aber ich habe Probleme, diese Lösung zu verallgemeinern.

Lassen die Anzahl der Kugeln der Farbe bezeichnen . Außerdem bezeichnen und die minimale bzw. maximale Anzahl von Kugeln der Farbe . Wir wollen gleichmäßig zufällig folgende Einschränkungen gelten:C i m i M i C i ( n 1 , … , n I$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$

Zunächst können Sie die Einschränkungen der unteren Grenze (dh ) entfernen, indem Sie zunächst Kugeln der Farbe . Dadurch werden die beiden Einschränkungen wie folgt geändert:m i C $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$

Lassen bezeichnen die einheitliche Verteilung , dass wir interessiert sind , können wir Kettenregel und dynamische Programmierung Probe verwenden. effizient. Zuerst verwenden wir die Kettenregel , um wie folgt zu schreiben : wobei ist die Randverteilung von . Man beachte, dassP P P ( n 1 , … , n I ∣ B , b 1 : I$P(n_1, \dots, n_I \mid B, b_{1:I})$$P$$P$ P(nI|B,b1:I)=∑ n 1 ,…, n I - 1 P(n1,…,nI|B,b1:I)nIP(nI|B,b1:I)nIB

$$\begin{align} P(n_1, \dots, n_I \mid B, b_{1:I}) &= P(n_I \mid B, b_{1:I}) P(n_1, \dots, n_{I-1} \mid n_I, B, b_{1:I}) \\ &= P(n_I \mid B, b_{1:I}) P(n_1, \dots, n_{I-1} \mid B-n_I, b_{1:I-1}) \quad (1) \end{align}$$ $P(n_I | B, b_{1:I}) = \sum_{n_1, \dots, n_{I-1}} P(n_1, \dots, n_I | B, b_{1:I})$$n_I$$P(n_I | B, b_{1:I})$ist eine diskrete Verteilung und kann mithilfe dynamischer Programmierung effizient berechnet werden. Beachten Sie auch, dass der zweite Term in (1) rekursiv berechnet werden kann. Wir probieren in der ersten Runde der Rekursion, aktualisieren die Gesamtzahl der Bälle auf und kehren in der nächsten Runde zu Probe zurück.$n_I$n I - $B - n_I$$n_{I-1}$

Das Folgende ist eine Matlab-Implementierung der Idee. Die Komplexität des Algorithmus ist wobei . Der Code verwendet zufällig generierte s in jedem Lauf. Infolgedessen haben einige der generierten Testfälle möglicherweise keine gültige Lösung. In diesem Fall druckt der Code eine Warnmeldung aus.$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

wo die Funktion print_constraints ist

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

und die Funktion dpfunc führt die dynamische Programmierberechnung wie folgt durch:

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

und schließlich zieht die Funktion diskrete Stichprobe (p, 1) eine Zufallsstichprobe aus der diskreten Verteilung . Eine Implementierung dieser Funktion finden Sie hier .$p$

Betrachten wir eine Verallgemeinerung dieses Problems. Es gibt Farbdosen mit verschiedenen Farben und Kugeln. Kann bis zu Bälle halten? Sie möchten Konfigurationen von Kugeln in den Dosen mit mindestens Kugeln in der Dose für jedes erzeugen , wobei jede Konfiguration mit gleicher Wahrscheinlichkeit erfolgt.$m=4$$m=4$$n^{(0)}=100$$i$$a^{(0)}_i = (100, 100, 50, 75)$$b_i = (0, 50, 0, 25)$$i$$i$

Solche Konfigurationen stimmen eins zu eins mit den Konfigurationen die nach dem Entfernen von Kugeln aus der Dose , wodurch verbleibende Bälle zu höchstens pro Dose. Ich werde diese daher einfach generieren und Sie sie anschließend anpassen lassen (indem Sie Bälle für jedes wieder in can legen ). i n = n ( 0 ) - ∑ i b i = 100 - ( 0 + 50 + 0 + 25 ) = 25 a i = a ( 0 ) i - b i = ( 100 , 50 , 50 , 50 ) b i i $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$

Um diese Konfigurationen hochzuzählen, korrigieren Sie alle bis auf zwei Indizes, z. B. und . Angenommen, es gibt bereits Bälle in Dose für jedes sich von und . Das lässt Bälle. Abhängig davon, wo sich die anderen Kugeln befinden, sind diese gleichmäßig in den Dosen und . Die möglichen Konfigurationen sind in der Anzahl (siehe die Kommentare) und reichen von der Platzierung so vieler Bälle in Dose$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$wie möglich den ganzen Weg so viele Bälle in Dose durch Platzieren wie möglich.$j$

Wenn Sie möchten, können Sie die Gesamtzahl der Konfigurationen zählen, indem Sie dieses Argument rekursiv auf die verbleibenden Dosen anwenden . Um Proben zu erhalten, müssen wir diese Anzahl jedoch nicht einmal kennen. Alles, was wir tun müssen, ist, wiederholt alle möglichen ungeordneten Paare von Dosen zu besuchen und die Verteilung der Kugeln innerhalb dieser beiden Dosen zufällig (und gleichmäßig) zu ändern. Dies ist eine Markov-Kette mit einer begrenzten Wahrscheinlichkeitsverteilung, die über alle möglichen Zustände gleichmäßig ist (wie mit Standardmethoden leicht gezeigt werden kann). Daher reicht es aus, in einem zu beginnen$m-2$$\{i,j\}$Führen Sie die Kette so lange aus, bis die Grenzverteilung erreicht ist, und verfolgen Sie dann die von diesem Verfahren besuchten Zustände. Wie üblich sollte diese Folge von Zuständen "verdünnt" werden, um eine serielle Korrelation zu vermeiden, indem sie übersprungen (oder zufällig wiederholt werden). Das Ausdünnen um den Faktor etwa der Hälfte der Anzahl der Dosen funktioniert in der Regel gut, da nach durchschnittlich vielen Schritten jede Dose betroffen ist und eine wirklich neue Konfiguration entsteht.

Dieser Algorithmus kostet Aufwand, um durchschnittlich jede zufällige Konfiguration zu generieren. Obwohl andere -Algorithmen existieren, hat dieser den Vorteil, dass die kombinatorischen Berechnungen nicht vorher durchgeführt werden müssen.$O(m)$$O(m)$

---

Lassen Sie uns als Beispiel eine kleinere Situation manuell herausarbeiten. Sei zum Beispiel und . Es gibt 15 gültige Konfigurationen, die als Zeichenfolgen von Belegungsnummern geschrieben werden können. Zum Beispiel werden zwei Kugeln in die zweite Dose und eine Kugel in der vierten Dose. Betrachten wir das Argument und betrachten wir die Gesamtbelegung der ersten beiden Dosen. Wenn das Bälle ist, bleiben für die letzten beiden Dosen keine Bälle übrig. Das gibt die Staaten$a=(4,3,2,1)$$n=3$0201$s_1+s_2=3$

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

wobei **alle möglichen Belegungsnummern für die letzten beiden Dosen darstellt: nämlich 00. Wenn , sind die Zustände$s_1+s_2=2$

20**, 11**, 02**

wo **kann jetzt entweder 10oder sein 01. Das ergibt weitere Zustände. Wenn , sind die Zustände$3\times 2=6$$s_1+s_2=1$

10**, 01**

wo jetzt **sein kann 20, 11aber nicht 02 (wegen der Grenze von einer Kugel in der letzten Dose). Das ergibt weitere Zustände. Wenn schließlich , befinden sich alle Kugeln in den letzten beiden Dosen, die bis zu ihren Grenzen von und voll sein müssen . Die gleich wahrscheinlichen Zustände sind daher$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.

Unter Verwendung des folgenden Codes wurde eine Folge von solcher Konfigurationen erzeugt und auf jede dritte verdünnt, wodurch Konfigurationen der Zustände erzeugt wurden. Ihre Frequenzen waren die folgenden:$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 

Ein Test der Gleichmäßigkeit ergibt einen Wert von , ( Freiheitsgrade): Das ist eine schöne Übereinstimmung mit der Hypothese, dass dieses Verfahren gleich wahrscheinliche Zustände erzeugt.$\chi^2$$\chi^2$$11.2$$p=0.67$$14$

---

Dieser RCode ist so eingerichtet, dass er die jeweilige Situation behandelt. Ändern Sie sich aund narbeiten Sie mit anderen Situationen. Stellen NSie den Wert so ein, dass die Anzahl der nach dem Ausdünnen erforderlichen Realisierungen generiert wird .

Dieser Code betrügt ein wenig, indem er systematisch alle Paare durchläuft . Wenn Sie über das Ausführen der Markow - Kette zu streng sein wollen, erzeugen , und zufällig, in dem kommentierten Code als gegeben. $(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)