Este es un seguimiento de una pregunta de Stackoverflow acerca de mezclar aleatoriamente una matriz .

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

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

Aquí está el algoritmo:

Haga un bucle un par de veces (la longitud de la matriz debería hacerlo), y en cada iteración, obtenga dos índices de matriz aleatorios e intercambie 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 de la matriz")?

Está roto, aunque si realiza suficientes barajaduras puede ser una aproximación excelente (como lo han indicado las respuestas anteriores).

Solo para tener una idea de lo que está sucediendo, considere con qué frecuencia su algoritmo generará barajaduras de una matriz de elementos en la que se arregla el primer elemento, k ≥ 2 . Cuando las permutaciones se generan con la misma probabilidad, esto debería suceder 1 / k del tiempo. Deje que p n sea ​​la frecuencia relativa de esta ocurrencia después de que n baraje con su algoritmo. Seamos generosos, también, y supongamos que en realidad estás seleccionando pares distintos de índices de manera uniforme al azar para tus barajas, de modo que cada par se seleccione con probabilidad 1 / ($k$$k \ge 2$$1/k$$p_n$$n$ =2/(k(k-1)). (Esto significa que no hay desperdicios "triviales" desperdiciados. Por otro lado, rompe totalmente su algoritmo para una matriz de dos elementos, porque alterna entre arreglar los dos elementos e intercambiarlos, así que si se detiene después de un número predeterminado de pasos, no hay aleatoriedad para el resultado en absoluto!)$1/{k \choose 2}$$2/\left( k (k-1) \right)$

Esta frecuencia satisface una recurrencia simple, porque el primer elemento se encuentra en su lugar original después de que baraja de dos maneras distintas. Una es que se solucionó después de n aleatorios y el siguiente aleatorio no mueve el primer elemento. La otra es que se movió después de n barajar, pero el n + 1 s t barajar lo mueve hacia atrás. La posibilidad de no mover el primer elemento es igual a ( k - $n+1$$n$$n$$n+1^{st}$ =(k-2)/k, mientras que la posibilidad de mover el primer elemento hacia atrás es igual a1/ (${k-1 \choose 2}/{k \choose 2}$$(k-2)/k$ =2/(k(k-1)). De dónde:$1/{k \choose 2}$$2/\left( k (k-1) \right)$

porque el primer elemento comienza en el lugar que le corresponde;

La solucion es

Restando , vemos que la frecuencia es incorrecta por ( k - $1/k$ . Parakyngrandes, una buena aproximación esk-$\left( \frac{k-3}{k-1} \right) ^n \frac{k-1}{k}$$k$$n$. Esto muestra que el erroren esta frecuencia particulardisminuirá exponencialmente con el número de intercambios en relación con el tamaño del conjunto (n/k), lo que indica que será difícil de detectar con conjuntos grandes si ha realizado un número relativamente grande de intercambios --pero el error siempre está ahí.$\frac{k-1}{k} \exp(-\frac{2n}{k-1})$$n/k$

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

donde debe ser muy pequeño en comparación con 1 / k . Esto implica que n debería ser varias veces k incluso para aproximaciones crudas ( es decir , donde ϵ está en el orden de 0.01 veces 1 / k más o menos).$\epsilon$$1/k$$n$$k$$\epsilon$$0.01$$1/k$

Todo esto plantea la pregunta: ¿por qué elegiría usar un algoritmo que no es del todo (pero solo aproximadamente) correcto, emplea exactamente las mismas técnicas que otro algoritmo que es demostrablemente correcto y, sin embargo, requiere más cómputo?

Editar

El comentario de Thilo es acertado (y esperaba que nadie lo señalara, ¡así podría ahorrarme este trabajo extra!). Déjame explicarte la lógica.

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

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

• Observe que el tamaño del error disminuye monotónicamente con el número de intercambios distintos. Por lo tanto, realizar menos intercambios en promedio también aumenta el error, en promedio. Pero este es un precio que debe estar dispuesto a pagar para superar el problema descrito en la primera viñeta. 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 excepción aparente interesante: una mirada cercana a la fórmula de error sugiere que no hay error en el caso . Esto no es un error: es correcto. Sin embargo, aquí he examinado solo una estadística relacionada con la distribución uniforme de 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 uniformemente. De hecho, después de 2 n intercambios reales, las únicas permutaciones posibles que se pueden generar son ( 123 ) , $k=3$$k=3$$2n$$(123)$ , y la identidad. Solo el último fija una posición determinada, de modo que exactamente un tercio de las permutaciones arreglan una posición. ¡Pero faltan la mitad de las permutaciones! En el otro caso, después de 2 n + 1 intercambios reales, las únicas permutaciones posibles son ( 12 ) , ( 23 ) y ( 13 ) . Una vez más, exactamente uno de estos arreglará cualquier posición, por lo que nuevamente obtenemos la frecuencia correcta de permutaciones que fijan esa posición, pero nuevamente obtenemos solo la mitad de las permutaciones posibles.$(321)$$2n+1$$(12)$$(23)$$(13)$

Este pequeño ejemplo ayuda a revelar las principales líneas del argumento: al ser "generosos" subestimamos conservadoramente la tasa de error para una estadística particular. Debido a que esa tasa de error no es cero para todos los , vemos que el algoritmo está roto. Además, al analizar la disminución en la tasa de error para esta estadística , establecemos un límite inferior en el número de iteraciones del algoritmo necesarias para tener alguna esperanza de aproximarnos a una distribución uniforme de permutaciones.$k \ge 4$

Creo que su algoritmo simple barajará las cartas correctamente ya que el número baraja tiende al infinito.

Supongamos que tiene tres cartas: {A, B, C}. Suponga que sus tarjetas comienzan en el siguiente orden: A, B, C. Luego, después de una mezcla, 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.

Usando 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).

Al codificar esto en R (ver el código a continuación), da la probabilidad de que la tarjeta A esté en la posición {1,2,3} como {0.33334, 0.33333, 0.33333} después de diez barajas.

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

$1/n!$$t$$A/n^{2t}$$A$$1/n! = A/n^{2t}$$n^{2t}/n! = A$$n \ge 3$$n$$n^{2t}/n!$$n!$$n=52$$1/52!$$3, 5, 7, ..., 47$$1/52^{2t}$$A/52^{2t}$$1/52!$

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

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

$\frac 12 n \log n$$(1-\epsilon) \frac12 n \log n$$(1+\epsilon) \frac 12 n \log n$$L^2$$7$

Tenga en cuenta que no soy estadístico, pero pondré mis 2 centavos.

Hice una pequeña prueba en R (cuidado, es muy lento para alto numTrials, el código probablemente se puede optimizar):

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 swapscon numTrials+1filas (una por prueba + la original) y numElementscolumnas (una por cada elemento vectorial). Si el método es correcto, la distribución de cada columna (es decir, de los valores de cada elemento en los ensayos) no debe ser diferente de la distribución de los datos originales.

Debido a que nuestros datos originales se distribuyeron normalmente, esperaríamos que todas las columnas no se desviaran de eso.

Si corremos

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

Obtenemos:

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

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

Lo que nos da p = 0.9926

Si verificamos 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")

obtenemos:

Entonces, 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

Tenga en cuenta que, obviamente, con un número menor de pruebas, la situación no es tan buena:

50 ensayos

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

100 ensayos

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

500 ensayos

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

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, a saber, 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. Llame a una de estas listas un rastro del programa.

Eso significa que hay $length ^ {2 \times length \times num\_passes}$tales rastros, y cada rastro es igualmente probable. También podemos asociar con cada traza una permutación de la matriz. A saber, la permutación al final de la ejecución asociada con la traza.

Existen $length !$ posibles permutaciones. $length ! < length ^ {2 \times length \times num\_passes}$ entonces, en general, una permutación dada está asociada con más de un rastro.

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

Elige cualquier prima $p$ tal que $p < length$, pero tal que $p \nmid length$, que puedes hacer por cualquier $length > 2$. Luego$p \bigm| length!$ pero no divide $length ^ {2 \times length \times num\_passes}$. Resulta que$length ! \nmid length ^ {2 \times length \times num\_passes}$ y entonces todas las permutaciones no pueden ser igualmente probables si $length > 2$.

¿Existe tal prima? Sí. Si$length$ eran divisibles por todos los números primos $p < length$, luego $length-1$ debe ser primo, pero luego $length-1$ sería una prima que es menor pero no divide $length$.

Compare esto con Fisher-Yates. En la primera iteración, usted elige entre$length$opciones. La segunda iteración tiene$length-1$opciones, y así sucesivamente. En otras palabras, tienes$length !$ rastros, y $length! \bigm| length!$. No es difícil demostrar que cada traza resulta en una permutación diferente, y desde allí es fácil ver que Fisher-Yates genera cada permutación con la misma probabilidad.