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
, j
y 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)