Comment générer uniformément des échantillons au hasard à partir de plusieurs variables discrètes soumises à des contraintes?

Je voudrais générer un processus Monte Carlo pour remplir une urne avec N boules de couleurs I, C [i]. Chaque couleur C [i] a un nombre minimum et maximum de boules qui doivent être placées dans l'urne.

Par exemple, j'essaie de mettre 100 balles dans l'urne et je peux la remplir de quatre couleurs:

• Rouge - Minimum 0, maximum 100 # NB, le maximum réel ne peut pas être atteint.
• Bleu - Minimum 50, maximum 100
• Jaune - Minimum 0, maximum 50
• Vert - Minimum 25, maximum 75

Comment puis-je générer un N échantillons garantissant une distribution uniforme entre les résultats possibles?

J'ai vu des solutions à ce problème où les boules n'ont pas de minimum ou maximum, ou alternativement, ont les mêmes minima et maxima implicites. Voir par exemple, cette discussion sur un sujet légèrement différent:

Générer des poids uniformément répartis qui correspondent à l'unité?

Mais j'ai des problèmes pour généraliser cette solution.

Soit le nombre de boules de couleur . Soit également et le nombre minimal et maximal de boules de couleur , respectivement. Nous voulons échantillonner uniformément au hasard sous réserve des contraintes suivantes:$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$

Tout d'abord, vous pouvez supprimer les contraintes de limite inférieure (ie ) en choisissant initialement boules de couleur . Cela modifie les deux contraintes comme suit:$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$

Soit la distribution uniforme qui nous intéresse. Nous pouvons utiliser la règle de chaîne et la programmation dynamique pour échantillonner efficacement à partir deTout d'abord, nous utilisons la règle de chaîne pour écrire comme suit: où est la distribution marginale de . Notez que$P(n_1, \dots, n_I \mid B, b_{1:I})$$P$$P$

$$\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})$est une distribution discrète et peut être calculée efficacement en utilisant la programmation dynamique. Notez également que le deuxième terme de (1) peut être calculé récursivement. Nous échantillonnons dans le premier tour de la récursivité, à jour le nombre total de boules à et pour échantillonner au tour suivant.$n_I$$B - n_I$$n_{I-1}$

Ce qui suit est une implémentation Matlab de l'idée. La complexité de l'algorithme est où . Le code utilise des générés aléatoirement à chaque exécution. Par conséquent, certains des cas de test générés peuvent ne pas avoir de solution valide, auquel cas le code imprime un message d'avertissement.$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

où la fonction print_constraints est

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

et la fonction dpfunc effectue le calcul de programmation dynamique comme suit:

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

et enfin, la fonction discreteample (p, 1) tire un échantillon aléatoire de la distribution discrète . Vous pouvez trouver une implémentation de cette fonction ici .$p$

Considérons une généralisation de ce problème. Il y a pots de peinture de couleurs distinctes et boules. Puis- tenir jusqu'à balles. Vous souhaitez générer des configurations de billes dans les canettes ayant au moins billes en can pour chaque , chaque configuration avec une probabilité égale.$m=4$$m=4$$n^{(0)}=100$$i$$a^{(0)}_i = (100, 100, 50, 75)$$b_i = (0, 50, 0, 25)$$i$$i$

De telles configurations sont en correspondance biunivoque avec les configurations obtenues après avoir boules de la boîte , limitant les balles restantes à au plus par boîte. Je vais donc simplement les générer et vous laisser les ajuster ensuite (en remettant les boules dans le can pour chaque ).$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$

Pour compter ces configurations, corrigez tous les indices sauf deux, disons et . Supposons qu'il y ait déjà balles dans la boîte pour chaque différent de et . Cela laisse des boules . À condition que se trouvent les autres boules, celles-ci sont réparties uniformément dans les boîtes et . Les configurations possibles sont en nombre (voir les commentaires), allant de la mise en place d'autant de balles dans le can$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$autant que possible tout en plaçant le plus de balles possible dans la boîte .$j$

Si vous le souhaitez, vous pouvez compter le nombre total de configurations en appliquant cet argument de manière récursive aux boîtes restantes . Cependant, pour obtenir des échantillons, nous n'avons même pas besoin de connaître ce nombre. Tout ce que nous devons faire est de visiter à plusieurs reprises toutes les paires possibles non ordonnées de boîtes et de modifier de manière aléatoire (et uniforme) la distribution des boules à l'intérieur de ces deux boîtes. Il s'agit d'une chaîne de Markov avec une distribution de probabilité limite uniforme sur tous les états possibles (comme cela est facilement montré en utilisant des méthodes standard). Par conséquent , il suffit de commencer dans une$m-2$$\{i,j\}$, exécutez la chaîne suffisamment longtemps pour atteindre la distribution limite, puis suivez les états visités par cette procédure. Comme d'habitude, pour éviter une corrélation sérielle, cette séquence d'états doit être "amincie" en les sautant (ou revisitée au hasard). L'amincissement d'un facteur d'environ la moitié du nombre de boîtes a tendance à bien fonctionner, car après cela, de nombreuses étapes en moyenne, chaque boîte a été affectée, produisant une configuration véritablement nouvelle.

Cet algorithme coûte effort pour générer chaque configuration aléatoire en moyenne. Bien qu'il existe d'autres algorithmes , celui-ci présente l'avantage de ne pas avoir à effectuer au préalable les calculs combinatoires.$O(m)$$O(m)$

---

Par exemple, travaillons manuellement sur une situation plus petite. Soit et , par exemple. Il existe 15 configurations valides, qui peuvent être écrites sous forme de chaînes de nombres d'occupation. Par exemple, place deux balles dans la deuxième boîte et une balle dans la quatrième boîte. En émulant l'argument, considérons l'occupation totale des deux premières boîtes. Lorsque c'est balles, aucune balle n'est laissée pour les deux dernières boîtes. Cela donne aux États$a=(4,3,2,1)$$n=3$0201$s_1+s_2=3$

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

où **représente tous les numéros d'occupation possibles pour les deux dernières boîtes: à savoir 00. Lorsque , les états sont$s_1+s_2=2$

20**, 11**, 02**

où maintenant **peut être soit 10ou 01. Cela donne états supplémentaires. Lorsque , les états sont$3\times 2=6$$s_1+s_2=1$

10**, 01**

où maintenant **peut être 20, 11mais pas 02 (en raison de la limite d'une balle dans la dernière boîte). Cela donne états supplémentaires. Enfin, lorsque , toutes les boules sont dans les deux dernières boîtes, qui doivent être pleines jusqu'à leurs limites de et . Les états également probables sont donc$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.

En utilisant le code ci-dessous, une séquence de ces configurations a été générée et réduite à un tiers, créant configurations des états. Leurs fréquences étaient les suivantes:$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 

A tests d'uniformité donne une valeur de , ( degrés de liberté): ce qui est beau accord avec l'hypothèse que cette procédure produit également des états probables.$\chi^2$$\chi^2$$11.2$$p=0.67$$14$

---

Ce Rcode est configuré pour gérer la situation dans la question. Changer aet ntravailler avec d'autres situations. Défini Npour être suffisamment grand pour générer le nombre de réalisations dont vous avez besoin après l'amincissement .

Ce code triche un peu en parcourant systématiquement toutes les paires . Si vous voulez être stricte sur l' exécution de la chaîne de Markov, générer , et au hasard, comme indiqué dans le code commenté. $(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)