¿Qué tiene de malo este algoritmo "ingenuo" de barajado?

Question

Esto es un seguimiento de un Stackoverflow pregunta sobre cómo barajar una matriz al azar .

Existen algoritmos establecidos (como el Baraja Knuth-Fisher-Yates ) que uno debería usar para barajar un array, en lugar de confiar en implementaciones "ingenuas" ad-hoc.

Ahora estoy interesado en probar (o refutar) que mi algoritmo ingenuo está roto (como en: no genera todas las permutaciones posibles con igual probabilidad).

Este es el algoritmo:

Haz un bucle un par de veces (la longitud del array debería bastar), y en cada iteración, obtén dos índices de array al azar e intercambia los dos elementos allí.

Obviamente, esto necesita más números aleatorios que KFY (el doble), pero aparte de eso, ¿funciona correctamente? ¿Y cuál sería el número apropiado de iteraciones (es suficiente la "longitud del array")?

Answer 1

Está roto, aunque si se realizan suficientes barajadas puede ser una excelente aproximación (como han indicado las respuestas anteriores).

Sólo para tener una idea de lo que está pasando, considere la frecuencia con la que su algoritmo generará barajadas de un $k$ matriz de elementos en la que el primer elemento es fijo, $k \ge 2$ . Cuando las permutaciones se generan con igual probabilidad, esto debería ocurrir $1/k$ de la época. Deje que $p_n$ sea la frecuencia relativa de esta ocurrencia después de $n$ baraja con su algoritmo. Seamos generosos, también, y supongamos que realmente está seleccionando distintivo pares de índices uniformemente al azar para sus barajados, de modo que cada par se selecciona con probabilidad $1/{k \choose 2}$ = $2/\left( k (k-1) \right)$ . (Esto significa que no hay barajadas "triviales" desperdiciadas. Por otro lado, rompe totalmente su algoritmo para una matriz de dos elementos, porque usted alterna entre la fijación de los dos elementos y el intercambio de ellos, por lo que si se detiene después de un número predeterminado de pasos, no hay ninguna aleatoriedad en el resultado).

Esta frecuencia satisface una simple recurrencia, porque el primer elemento se encuentra en su lugar original después de $n+1$ baraja de dos maneras disjuntas. Una es que se fijó después de $n$ baraja y la siguiente baraja no mueve el primer elemento. La otra es que se haya movido después de $n$ baraja pero el $n+1^{st}$ La baraja se mueve hacia atrás. La posibilidad de no moviendo el primer elemento igual a ${k-1 \choose 2}/{k \choose 2}$ = $(k-2)/k$ mientras que la posibilidad de desplazar el primer elemento hacia atrás es igual a $1/{k \choose 2}$ = $2/\left( k (k-1) \right)$ . De donde:

$$p_0 =1$$ porque el primer elemento comienza en el lugar que le corresponde;

$$p_{n+1} = \frac{k-2}{k} p_n + \frac{2}{k(k-1)} \left( 1 - p_n \right).$$

La solución es

$$p_n = 1/k + \left( \frac{k-3}{k-1} \right) ^n \frac{k-1}{k}.$$

Restando $1/k$ vemos que la frecuencia es errónea por $\left( \frac{k-3}{k-1} \right) ^n \frac{k-1}{k}$ . Para los grandes $k$ y $n$ una buena aproximación es $\frac{k-1}{k} \exp(-\frac{2n}{k-1})$ . Esto demuestra que el error en esta frecuencia particular disminuirá exponencialmente con el número de intercambios en relación con el tamaño de la matriz ( $n/k$ ), lo que indica que será difícil de detectar con matrices grandes si ha hecho un número relativamente grande de intercambios, pero el error siempre está ahí.

Es difícil ofrecer un análisis exhaustivo de los errores en todas las frecuencias. Sin embargo, es probable que se comporten como éste, que muestra que como mínimo se necesitaría $n$ (el número de intercambios) sea lo suficientemente grande como para que el error sea aceptablemente pequeño. Una solución aproximada es

$$n \gt \frac{1}{2} \left(1 - (k-1) \log(\epsilon) \right)$$

donde $\epsilon$ debe ser muy pequeño en comparación con $1/k$ . Esto implica $n$ debe ser varias veces $k$ incluso para las aproximaciones más burdas ( es decir , donde $\epsilon$ es del orden de $0.01$ veces $1/k$ más o menos).

Todo esto nos lleva a preguntarnos: ¿por qué elegir un algoritmo que no es del todo (sino sólo aproximadamente) correcto, que emplea exactamente las mismas técnicas que otro algoritmo que sí es demostrablemente correcto y que, sin embargo, requiere más cálculos?

Editar

El comentario de Thilo es acertado (y esperaba que nadie lo señalara, para ahorrarme este trabajo extra). Permítanme explicar la lógica.

• Si te aseguras de generar intercambios reales cada vez, estás totalmente jodido. El problema que señalé para el caso $k=2$ se extiende a todas las matrices. Sólo la mitad de todas las permutaciones posibles pueden obtenerse aplicando un número par de permutas; la otra mitad se obtiene aplicando un número impar de permutas. Por lo tanto, en esta situación, nunca se puede generar ni de lejos una distribución uniforme de permutaciones (pero hay tantas posibles que un estudio de simulación para cualquier $k$ no podrá detectar el problema). Eso es muy malo.

• Por lo tanto, es prudente generar intercambios al azar generando las dos posiciones de forma independiente al azar. Esto significa que hay un $1/k$ probabilidad cada vez de intercambiar un elemento consigo mismo; es decir, de no hacer nada. Este proceso ralentiza un poco el algoritmo: después de $n$ pasos, esperamos sólo unos $\frac{k-1}{k} N \lt N$ verdaderos intercambios que se han producido.

• Obsérvese que el tamaño del error disminuye monótonamente con el número de intercambios distintos. Por lo tanto, la realización de menos intercambios en promedio también aumenta el error, por término medio. Pero este es un precio que debería estar dispuesto a pagar para superar el problema descrito en el primer punto. En consecuencia, mi estimación de error es conservadoramente baja, aproximadamente por un factor de $(k-1)/k$ .

También quería señalar una interesante excepción aparente: un examen detallado de la fórmula de error sugiere que hay no error en el caso $k=3$ . Esto no es un error: es correcto. Sin embargo, aquí sólo he examinado una estadística relacionada con la distribución uniforme de las permutaciones. El hecho de que el algoritmo pueda reproducir esta estadística cuando $k=3$ (es decir, obtener la frecuencia correcta de permutaciones que fijan una posición determinada) no garantiza que las permutaciones se hayan distribuido efectivamente de manera uniforme. En efecto, después de $2n$ intercambios reales, las únicas permutaciones posibles que se pueden generar son $(123)$ , $(321)$ y la identidad. Sólo esta última fija una posición determinada, por lo que exactamente un tercio de las permutaciones fijan una posición. Pero falta la mitad de las permutaciones. En el otro caso, después de $2n+1$ intercambios reales, las únicas permutaciones posibles son $(12)$ , $(23)$ et $(13)$ . De nuevo, exactamente una de ellas fijará una posición determinada, por lo que de nuevo obtenemos la frecuencia correcta de permutaciones que fijan esa posición, pero de nuevo obtenemos sólo la mitad de las permutaciones posibles.

Este pequeño ejemplo ayuda a revelar las principales líneas del argumento: al ser "generosos" conservamos subestimar la tasa de error de una estadística concreta. Como esa tasa de error es distinta de cero para todos los $k \ge 4$ vemos que el algoritmo está roto. Además, analizando el decaimiento de la tasa de error para esta estadística establecemos un límite inferior sobre el número de iteraciones del algoritmo necesarias para tener alguna esperanza de aproximarse a una distribución uniforme de permutaciones.

Answer 2

Creo que tu sencillo algoritmo barajará las cartas correctamente ya que el número de barajadas tiende al infinito.

Supongamos que tienes tres cartas: {A,B,C}. Suponga que sus cartas comienzan en el siguiente orden: A,B,C. Entonces, después de una barajada tienes las siguientes combinaciones:

{A,B,C}, {A,B,C}, {A,B,C} #You get this if choose the same RN twice.
{A,C,B}, {A,C,B}
{C,B,A}, {C,B,A}
{B,A,C}, {B,A,C}

Por lo tanto, la probabilidad de que la carta A esté en la posición {1,2,3} es {5/9, 2/9, 2/9}.

Si barajamos las cartas por segunda vez, entonces:

Pr(A in position 1 after 2 shuffles) = 5/9*Pr(A in position 1 after 1 shuffle) 
                                     + 2/9*Pr(A in position 2 after 1 shuffle) 
                                     + 2/9*Pr(A in position 3 after 1 shuffle) 

Esto da 0,407.

Utilizando la misma idea, podemos formar una relación de recurrencia, es decir

Pr(A in position 1 after n shuffles) = 5/9*Pr(A in position 1 after (n-1) shuffles) 
                                     + 2/9*Pr(A in position 2 after (n-1) shuffles) 
                                     + 2/9*Pr(A in position 3 after (n-1) shuffles).

Codificando esto en R (ver código abajo), da la probabilidad de que la carta A esté en la posición {1,2,3} como {0,33334, 0,33333, 0,33333} después de diez barajadas.

Código R

## m is the probability matrix of card position
## Row is position
## Col is card A, B, C
m = matrix(0, nrow=3, ncol=3)
m[1,1] = 1; m[2,2] = 1; m[3,3] = 1

## Transition matrix
m_trans = matrix(2/9, nrow=3, ncol=3)
m_trans[1,1] = 5/9; m_trans[2,2] = 5/9; m_trans[3,3] = 5/9

for(i in 1:10){
  old_m = m
  m[1,1] = sum(m_trans[,1]*old_m[,1])
  m[2,1] = sum(m_trans[,2]*old_m[,1])
  m[3,1] = sum(m_trans[,3]*old_m[,1])

  m[1,2] = sum(m_trans[,1]*old_m[,2])
  m[2,2] = sum(m_trans[,2]*old_m[,2])
  m[3,2] = sum(m_trans[,3]*old_m[,2])

  m[1,3] = sum(m_trans[,1]*old_m[,3])
  m[2,3] = sum(m_trans[,2]*old_m[,3])
  m[3,3] = sum(m_trans[,3]*old_m[,3])
}  
m

Answer 3

Una forma de ver que no se obtiene una distribución perfectamente uniforme es por la divisibilidad. En la distribución uniforme, la probabilidad de cada permutación es $1/n!$ . Cuando se genera una secuencia de $t$ transposiciones aleatorias, y luego recoger las secuencias por su producto, las probabilidades que se obtienen son de la forma $A/n^{2t}$ para algún número entero $A$ . Si $1/n! = A/n^{2t}$ entonces $n^{2t}/n! = A$ . Por el postulado de Bertrand (un teorema), para $n \ge 3$ hay primos que aparecen en el denominador y que no dividen $n$ Así que $n^{2t}/n!$ no es un número entero, y no hay manera de dividir las transposiciones uniformemente en $n!$ permutaciones. Por ejemplo, si $n=52$ entonces el denominador de $1/52!$ es divisible por $3, 5, 7, ..., 47$ mientras que el denominador de $1/52^{2t}$ no lo es, así que $A/52^{2t}$ no puede reducirse a $1/52!$ .

¿Cuántos se necesitan para aproximar bien una permutación aleatoria? La generación de una permutación aleatoria mediante transposiciones aleatorias fue analizada por Diaconis y Shahshahani utilizando la teoría de la representación del grupo simétrico en

Diaconis, P., Shahshahani, M. (1981): "Generación de una permutación aleatoria con transposiciones aleatorias". Z. Wahrsch. Verw. Geb. 57, 159-179.

Una de las conclusiones fue que se necesita $\frac 12 n \log n$ transposiciones en el sentido de que después de $(1-\epsilon) \frac12 n \log n$ las permutaciones están lejos de ser aleatorias, pero después $(1+\epsilon) \frac 12 n \log n$ el resultado es casi aleatorio, tanto en el sentido de la variación total como $L^2$ distancia. Este tipo de fenómeno de corte es común en los paseos aleatorios sobre grupos, y está relacionado con el famoso resultado de que se necesita $7$ barajar antes de que una baraja se acerque al azar.

Answer 4

Tengan en cuenta que no soy un estadístico, pero pondré mis dos centavos.

He hecho una pequeña prueba en R (cuidado, es muy lento para altas numTrials El código puede ser optimizado):

numElements <- 1000
numTrials <- 5000

swapVec <- function()
    {
    vec.swp <- vec

    for (i in 1:numElements)
        {
        i <- sample(1:numElements)
        j <- sample(1:numElements)

        tmp <- vec.swp[i]
        vec.swp[i] <- vec.swp[j]
        vec.swp[j] <- tmp
        }

    return (vec.swp)
    }

# Create a normally distributed array of numElements length
vec <- rnorm(numElements)

# Do several "swapping trials" so we can make some stats on them
swaps <- vec
prog <- txtProgressBar(0, numTrials, style=3)

for (t in 1:numTrials)
    {
    swaps <- rbind(swaps, swapVec())
    setTxtProgressBar(prog, t)
    }

Esto generará una matriz swaps con numTrials+1 filas (una por ensayo + el original) y numElements columnas (una por cada elemento del vector). Si el método es correcto, la distribución de cada columna (es decir, de los valores de cada elemento a lo largo de los ensayos) no debería ser diferente de la distribución de los datos originales.

Dado que nuestros datos originales se distribuyen normalmente, es de esperar que todas las columnas no se desvíen de ella.

Si ejecutamos

par(mfrow= c(2,2))
# Our original data
hist(swaps[1,], 100, col="black", freq=FALSE, main="Original")
# Three "randomly" chosen columns
hist(swaps[,1], 100, col="black", freq=FALSE, main="Trial # 1") 
hist(swaps[,257], 100, col="black", freq=FALSE, main="Trial # 257")
hist(swaps[,844], 100, col="black", freq=FALSE, main="Trial # 844")

Lo conseguimos:

que parece muy prometedor. Ahora bien, si queremos confirmar estadísticamente que las distribuciones no se desvían del original creo que podríamos utilizar una prueba de Kolmogorov-Smirnov (por favor, ¿puede algún estadístico confirmar que esto es correcto?) y hacer, por ejemplo

ks.test(swaps[1, ], swaps[, 234])

Lo que nos da p=0.9926

Si comprobamos todas las columnas:

ks.results <- apply(swaps, 2, function(col){ks.test(swaps[1,], col)})
p.values <- unlist(lapply(ks.results, function(x){x$p.value})

Y corremos

hist(p.values, 100, col="black")

nos encontramos con que:

Por lo tanto, para la gran mayoría de los elementos de la matriz, su método de intercambio ha dado un buen resultado, como también puede ver mirando los cuartiles.

1> quantile(p.values)
       0%       25%       50%       75%      100% 
0.6819832 0.9963731 0.9999188 0.9999996 1.0000000

Obsérvese que, obviamente, con un número menor de ensayos la situación no es tan buena:

50 pruebas

1> quantile(p.values)
          0%          25%          50%          75%         100% 
0.0003399635 0.2920976389 0.5583204486 0.8103852744 0.9999165730

100 pruebas

          0%         25%         50%         75%        100% 
 0.001434198 0.327553996 0.596603804 0.828037097 0.999999591 

500 pruebas

         0%         25%         50%         75%        100% 
0.007834701 0.504698404 0.764231550 0.934223503 0.999995887

Answer 5

Así es como estoy interpretando su algoritmo, en pseudocódigo:

void shuffle(array, length, num_passes)
  for (pass = 0; pass < num_passes; ++pass) 
    for (n = 0; n < length; ++)
      i = random_in(0, length-1)
      j = random_in(0, lenght-1)
      swap(array[i], array[j]

Podemos asociar una ejecución de este algoritmo con una lista de $2 \times length \times num\_passes$ enteros, es decir, los enteros devueltos por random_in() mientras se ejecuta el programa. Cada uno de estos enteros está en $[0, length-1]$ y también lo ha hecho $length$ valores posibles. Llama a una de estas listas un rastro del programa.

Eso significa que hay $length ^ {2 \times length \times num\_passes}$ tales rastros, y cada uno de ellos es igualmente probable. También podemos asociar a cada traza una permutación de la matriz. Es decir, la permutación al final de la carrera asociada a la traza.

Hay $length !$ posibles permutaciones. $length ! < length ^ {2 \times length \times num\_passes}$ por lo que, en general, una determinada permutación está asociada a más de una traza.

Recuerda que las trazas son todas igual de probables, por lo que para que todas las permutaciones sean igual de probables, dos permutaciones cualesquiera deben estar asociadas al mismo número de trazas. Si esto es cierto, entonces debemos tener $length ! \bigm| length ^ {2 \times length \times num\_passes}$ .

Escoge cualquier primo $p$ tal que $p < length$ pero de manera que $p \nmid length$ que se puede hacer para cualquier $length > 2$ . Entonces $p \bigm| length!$ pero no divide $length ^ {2 \times length \times num\_passes}$ . De ello se desprende que $length ! \nmid length ^ {2 \times length \times num\_passes}$ y por lo tanto todas las permutaciones no pueden ser igualmente probables si $length > 2$ .

¿Existe tal prima? Sí. Si $length$ fueran divisibles por todos los primos $p < length$ entonces $length-1$ debe ser primo, pero entonces $length-1$ sería un primo tal que es menor que pero no divide a $length$ .

Compara esto con Fisher-Yates. En la primera iteración, se elige entre $length$ opciones. La segunda iteración tiene $length-1$ opciones, etc. En otras palabras, tiene $length !$ trazos, y $length! \bigm| length!$ . No es difícil demostrar que cada traza da lugar a una permutación diferente, y a partir de ahí es fácil ver que Fisher-Yates genera cada permutación con igual probabilidad.