8 votos

¿Cómo generar muestras uniformemente al azar de múltiples variables discretas sujeta a limitaciones?

Me gustaría generar un Monte Carlo proceso de llenar una urna con N bolas de colores me, C[i]. Cada color C[i] tiene un mínimo y un máximo número de bolas que deben ser colocados en la urna.

Por ejemplo, estoy tratando de poner 100 bolas en la urna, y se puede llenar con cuatro colores:

  • Rojo - Mínimo 0, Máximo 100 # NB, la máxima real no puede ser observado.
  • Azul - Un Mínimo De 50, Máximo 100
  • Amarillo - Mínimo 0, Máximo 50
  • Verde - Con Un Mínimo De 25, Máximo 75

¿Cómo puedo generar un N muestras de que está garantizada para ser distribuidos de manera uniforme a través de los resultados posibles?

He visto soluciones a este problema en el que las bolas no tienen un mínimo o máximo, o, alternativamente, tiene el mismo implícita de mínimos y máximos. Véase, por ejemplo, esta discusión en un poco diferente:

Generar distribuida uniformemente en pesos que se suma a la unidad?

Pero estoy teniendo problemas para generalizar esta solución.

4voto

Thomas Puntos 1

Deje $n_i$ denotar el número de bolas del color $C_i$. También, vamos a $m_i$ $M_i$ denotar el mínimo y el máximo número de bolas del color $C_i$, respectivamente. Queremos muestra $(n_1, \dots, n_I)$ uniformemente al azar sujeción a las siguientes limitaciones:

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

Primero de todos, usted puede quitar el límite inferior de restricciones (es decir,$m_i \leq n_i$) mediante la selección de $m_i$ bolas de color $C_i$ inicialmente. Esto modifica las dos restricciones de la siguiente manera:

  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$

Deje $P(n_1, \dots, n_I \mid B, b_{1:I})$ denotar el uniforme de la distribución que nos interesa. Podemos utilizar la regla de la cadena y la dinámica de la programación de la muestra de $P$ eficiente. En primer lugar, utilizamos la regla de la cadena para escribir $P$ como sigue: $$ \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} $$ donde $P(n_I | B, b_{1:I}) = \sum_{n_1, \dots, n_{I-1}} P(n_1, \dots, n_I | B, b_{1:I})$ es la distribución marginal de $n_I$. Tenga en cuenta que $P(n_I | B, b_{1:I})$ es una distribución discreta y puede ser calculada de manera eficiente utilizando programación dinámica. También, tenga en cuenta que el segundo término de (1) puede ser calculada de forma recursiva. Nos muestra $n_I$ en la primera ronda de la recursividad, la actualización, el número total de bolas de a $B - n_I$ y recurse de la muestra $n_{I-1}$ en la siguiente ronda.

La siguiente es una implementación en Matlab de la idea. La complejidad del algoritmo es $O(I \times B \times K)$ donde $K = \max_{i=1}^I b_i$. Utiliza el código generado aleatoriamente $m_i$s en cada carrera. Como resultado, algunos de los que se generan casos de prueba, no puede tener ninguna solución válida, en cuyo caso el código se imprime un mensaje de advertencia.

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

donde la función print_constraints es

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

y la función dpfunc realiza la programación dinámica de cálculo de la siguiente manera:

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

y por último, la función discretesample(p, 1) dibuja una muestra aleatoria de la distribución discreta $p$. Usted puede encontrar una implementación de esta función aquí.

2voto

jldugger Puntos 7490

Vamos a considerar una generalización de este problema. Hay $m=4$ latas de pintura de $m=4$ distintos colores y $n^{(0)}=100$ bolas. Puede $i$ puede albergar hasta a $a^{(0)}_i = (100, 100, 50, 75)$ bolas. Se desea generar configuraciones de bolas en las latas de tener al menos $b_i = (0, 50, 0, 25)$ bolas en can $i$ por cada $i$, cada configuración con igual probabilidad.

Estas configuraciones están en una correspondencia uno a uno con las configuraciones que se obtiene después de eliminar a $b_i$ bolas de can $i$, limitando el $n = n^{(0)} - \sum_i b_i = 100 - (0+50+0+25)=25$ restante de bolas en la mayoría de las $a_i = a^{(0)}_i - b_i = (100, 50, 50, 50)$ puede. Por lo tanto, sólo generar estos y permiten ajustar después (por poner $b_i$ bolas de nuevo en can $i$ por cada $i$).

Para contar estas configuraciones, reparar todos, pero dos de los índices, de decir $i$$j$. Supongamos que hay $s_k$ bolas ya en can $k$ por cada $k$ diferentes de$i$$j$. Que deja a $s_i+s_j$ bolas. Condicional en caso de que el otro $n - (s_i+s_j)$ bolas son, estos se distribuyen de manera uniforme dentro de las latas $i$$j$. Las configuraciones posibles son $1 + \min(a_i + a_j - s_i - s_j, s_i+s_j)$ en número (ver los comentarios), que van desde la colocación de tantas bolas en can $i$ como sea posible, todo el camino a través de la colocación de tantas bolas en can $j$ como sea posible.

Si usted desea, usted puede contar el número total de configuraciones mediante la aplicación de este argumento de forma recursiva para el resto de $m-2$ latas. Sin embargo, para obtener muestras ni siquiera necesitamos saber este recuento. Todo lo que tenemos que hacer es repetidamente visitar todos los posibles pares no ordenados $\{i,j\}$ de las latas y de forma aleatoria (y de manera uniforme) cambiar la distribución de las bolas dentro de los dos latas. Esta es una cadena de Markov con un limitante de la distribución de probabilidad que es uniforme en todos los estados posibles (como es fácilmente demostrado con los métodos estándar). Por lo tanto, es suficiente para iniciar en cualquier estado, ejecución de la cadena lo suficientemente larga para llegar a la limitación de la distribución y, a continuación, seguir la pista de los estados visitados por este procedimiento. Como de costumbre, para evitar la correlación serial, esta secuencia de estados debe ser "adelgazada" por saltar a través de ellos (o revisado al azar). Adelgazamiento por un factor de alrededor de la mitad el número de latas tiende a funcionar bien, porque después de muchos pasos, en promedio, cada uno puede ha sido afectada, produciendo una verdadera nueva configuración.

Este algoritmo de costos $O(m)$ esfuerzo para generar cada aleatoria de configuración en promedio. Aunque otros $O(m)$ algoritmos de existir, este tiene la ventaja de no necesitar hacer la combinatoria de los cálculos de antemano.


Como ejemplo, veamos un pequeño situación manualmente. Deje $a=(4,3,2,1)$$n=3$, por ejemplo. Hay 15 configuraciones válidas, el cual puede ser escrito como cadenas de números de ocupación. Por ejemplo, 0201 coloca dos bolas en el segundo y una bola en el cuarto puede. Emulando el argumento, vamos a considerar el total de la ocupación de las dos primeras latas. Cuando es $s_1+s_2=3$ pelotas, las pelotas, se dejan para los últimos dos latas. Que da a los estados

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

donde ** representa todas las posibilidades de ocupación de los números de los últimos dos latas: a saber, 00. Al $s_1+s_2=2$, de los estados

20**, 11**, 02**

donde ahora ** puede ser 10 o 01. Que da $3\times 2=6$ más estados. Al $s_1+s_2=1$, de los estados

10**, 01**

donde ahora ** puede 20, 11, pero no 02 (debido a que el límite de una pelota en el pasado). Que da $2\times 2=4$ más estados. Finalmente, cuando se $s_1+s_2=0$, todas las bolas en los últimos dos latas, que debe ser completo a sus límites de $2$$1$. El $4+6+4+1=15$ igualmente probables por consiguiente, los estados son

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

Usando el código a continuación, una secuencia de $10,009$ dichas configuraciones se generó y diluido a cada tercero, la creación de $3337$ configuraciones de la $15$ estados. Sus frecuencias fueron los siguientes:

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 

Un $\chi^2$ prueba de homogeneidad da un $\chi^2$ valor de $11.2$, $p=0.67$ ($14$ grados de libertad): que es hermoso acuerdo con la hipótesis de que este procedimiento produce igualmente probables de los estados.


Este R código está configurado para manejar la situación en la pregunta. Cambio a y n trabajar con otras situaciones. Set N ser lo suficientemente grande como para generar el número de realizaciones que usted necesita después de adelgazamiento.

Este código de trucos un poco de ciclismo de forma sistemática a través de todos los $(i,j)$ pares. Si se quiere ser estricto sobre la ejecución de la cadena de Markov, generar i, jy ij al azar, como se indica en el código comentado.

#
# 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)

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by:

X