O que há de errado com esse algoritmo de embaralhamento “ingênuo”?

23

Este é um seguimento de uma pergunta do Stackoverflow sobre embaralhar uma matriz aleatoriamente .

Existem algoritmos estabelecidos (como o Knuth-Fisher-Yates Shuffle ) que se deve usar para embaralhar uma matriz, em vez de confiar em implementações ad-hoc "ingênuas".

Agora estou interessado em provar (ou refutar) que meu algoritmo ingênuo está quebrado (como em: não gera todas as permutações possíveis com igual probabilidade).

Aqui está o algoritmo:

Faça um loop algumas vezes (o comprimento da matriz deve funcionar) e, a cada iteração, obtenha dois índices aleatórios da matriz e troque os dois elementos.

Obviamente, isso precisa de mais números aleatórios que o KFY (o dobro), mas, além disso, ele funciona corretamente? E qual seria o número apropriado de iterações (o "comprimento da matriz" é suficiente)?

Thilo
fonte
4
Eu simplesmente não consigo entender por que as pessoas pensam que essa troca é 'mais simples' ou 'mais ingênua' que o EF ... Quando eu estava resolvendo esse problema pela primeira vez, acabei de implementar o AF (sem saber que ele tem nome) , apenas porque parecia a maneira mais simples de fazer isso por mim.
1
@mbq: pessoalmente, acho-os igualmente fáceis, embora eu concorde que o EF pareça mais "natural" para mim.
Nico
3
Quando pesquisei algoritmos de embaralhamento depois de escrever o meu próprio (uma prática que abandonei desde então), eu estava "merda, isso foi feito e tem um nome !"
JM não é estatístico

Respostas:

12

Está quebrado, embora se você executar baralhamento suficiente, pode ser uma excelente aproximação (como as respostas anteriores indicaram).

Apenas para entender o que está acontecendo, considere com que frequência seu algoritmo irá gerar embaralhamento de uma matriz de elementos na qual o primeiro elemento é fixo, . Quando permutações são geradas com igual probabilidade, isso deve ocorrer do tempo. Seja a frequência relativa dessa ocorrência após embaralhar com seu algoritmo. Sejamos generosos também, e suponha que você esteja realmente selecionando pares distintos de índices uniformemente aleatoriamente para seus embaralhamentos, de modo que cada par seja selecionado com probabilidade =k 2 1 / k p n n 1 / ( kkk21/kpnn 2/(k(k-1))1/(k2)2/(k(k1)). (Isso significa que não há desperdícios "triviais" desperdiçados. Por outro lado, ele interrompe totalmente seu algoritmo para uma matriz de dois elementos, porque você alterna entre fixar os dois elementos e trocá-los; portanto, se você parar após um número predeterminado de etapas, não há aleatoriedade para o resultado!)

Essa frequência satisfaz uma recorrência simples, porque o primeiro elemento é encontrado em seu lugar original após embaralhar de duas maneiras disjuntas. Uma é que ele foi corrigido após shuffles e o próximo shuffle não move o primeiro elemento. A outra é que ele foi movido após shuffles, mas o move para trás. A chance de não mover o primeiro elemento é igual a = , enquanto a chance de mover o primeiro elemento para trás é igual a = . De onde:n n n + 1 s t ( k - 1n+1nnn+1st (k-2)/k1/ ( k(k12)/(k2)(k2)/k 2/(k(k-1))1/(k2)2/(k(k1))

p0=1
porque o primeiro elemento começa em seu devido lugar;

pn+1=k-2kpn+2k(k-1)(1-pn).

A solução é

pn=1/k+(k-3k-1)nk-1k.

Subtraindo , vemos que a frequência está errada por . Para e grandes , uma boa aproximação é . Isso mostra que o erro nessa frequência específica diminuirá exponencialmente com o número de trocas em relação ao tamanho da matriz ( ), indicando que será difícil detectar com matrizes grandes se você tiver feito um número relativamente grande de trocas. - mas o erro está sempre lá.( k - 31/k knk-1(k-3k-1)nk-1kknn/kk-1kexp(-2nk-1)n/k

É difícil fornecer uma análise abrangente dos erros em todas as frequências. É provável que eles se comportem como este, o que mostra que, no mínimo, você precisaria de (o número de trocas) para ser grande o suficiente para tornar o erro aceitávelmente pequeno. Uma solução aproximada én

n>12(1(k1)log(ϵ))

onde deve ser muito pequeno comparado a . Isso implica que deve ser várias vezes para aproximações grosseiras ( ou seja , onde é da ordem de vezes ou mais).1 / k n k ϵ 0,01 1 / kϵ1/knkϵ0.011/k

Tudo isso levanta a questão: por que você escolheria usar um algoritmo que não é muito (mas apenas aproximadamente) correto, emprega exatamente as mesmas técnicas que outro algoritmo que é comprovadamente correto e, no entanto, que requer mais computação?

Editar

O comentário de Thilo é adequado (e eu esperava que ninguém apontasse isso, para que eu pudesse ser poupada desse trabalho extra!). Deixe-me explicar a lógica.

  • Se você gerar trocas reais a cada vez, estará totalmente ferrado. O problema que apontei para o caso se estende a todas as matrizes. Apenas metade de todas as permutações possíveis pode ser obtida aplicando um número par de swaps; a outra metade é obtida aplicando um número ímpar de swaps. Portanto, nessa situação, você nunca pode gerar em lugar algum uma distribuição uniforme de permutações (mas há tantas possíveis que um estudo de simulação para qualquer considerável não será capaz de detectar o problema). Isso é muito ruim.kk=2k

  • Portanto, é aconselhável gerar swaps aleatoriamente, gerando as duas posições independentemente, aleatoriamente. Isso significa que há uma chance de cada vez que um elemento é trocado; isto é, de não fazer nada. Esse processo efetivamente diminui um pouco o algoritmo: após etapas, esperamos que apenas cerca de trocas verdadeiras ocorram.n k - 11/knk1kN<N

  • Observe que o tamanho do erro diminui monotonicamente com o número de trocas distintas. Portanto, realizar menos swaps em média também aumenta o erro, em média. Mas este é um preço que você deve estar disposto a pagar para superar o problema descrito no primeiro item. Consequentemente, minha estimativa de erro é conservadoramente baixa, aproximadamente por um fator de .(k-1)/k

Eu também queria destacar uma exceção aparente interessante: uma análise mais detalhada da fórmula do erro sugere que não erro no caso . Isso não é um erro: está correto. No entanto, aqui examinei apenas uma estatística relacionada à distribuição uniforme de permutações. O fato de o algoritmo poder reproduzir esta estatística quando (ou seja, obter a frequência certa de permutações que fixam qualquer posição) não garante que as permutações tenham sido realmente distribuídas uniformemente. De fato, após swaps reais, as únicas permutações possíveis que podem ser geradas são ,k = 3 2 n ( 123 ) ( 321 ) 2 n + 1 ( 12 ) ( 23 ) ( 13 )k=3k=32n(123)(321)e a identidade. Somente o último fixa uma determinada posição; portanto, exatamente um terço das permutações fixa uma posição. Mas metade das permutações está faltando! No outro caso, após swaps reais, as únicas permutações possíveis são , e . Novamente, exatamente um deles fixará qualquer posição, então obteremos a frequência correta de permutações que fixam essa posição, mas novamente obteremos apenas metade das permutações possíveis.2n+1(12)(23)(13)

Este pequeno exemplo ajuda a revelar as principais linhas do argumento: por ser "generoso", subestimamos conservadoramente a taxa de erro de uma estatística específica. Como essa taxa de erro é diferente de zero para todos os , vemos que o algoritmo está quebrado. Além disso, analisando o decaimento na taxa de erro dessa estatística , estabelecemos um limite mais baixo para o número de iterações do algoritmo necessário para ter alguma esperança de aproximar uma distribuição uniforme de permutações.k4

whuber
fonte
1
"Vamos ser generosos também, e suponha que você esteja realmente selecionando pares distintos de índices de maneira uniforme e aleatória para seus shuffles". Não entendo por que essa suposição pode ser feita e como é generosa. Parece descartar possíveis permutações, resultando em uma distribuição ainda menos aleatória.
Thilo
1
@ Thilo: Obrigado. Seu comentário merece uma resposta estendida, então eu o coloquei na própria resposta. Deixe-me salientar aqui que ser "generoso" na verdade não descarta nenhuma permutação: apenas elimina etapas no algoritmo que de outra forma não fariam nada.
whuber
2
Esse problema pode ser analisado completamente como uma cadeia de Markov no gráfico de Cayley do grupo de permutação. Os cálculos numéricos para k = 1 a 7 (uma matriz 5040 por 5040!) Confirmam que os maiores valores próprios de tamanho (após 1 e -1) são exatamente . Isso implica que, depois de lidar com o problema de alternar o sinal da permutação (correspondente ao valor próprio de -1), os erros em todas as probabilidades decaem na taxa ou Mais rápido. Suspeito que isso continue valendo para todos os maiores . ( 1 - 2 / ( k - 1 ) ) n k(k3)/(k1)=12/(k1)(12/(k1))nk
whuber
1
Você pode fazer muito melhor que pois as probabilidades são invariantes nas classes de conjugação, e existem apenas partições de para que você possa analisar uma matriz . 15 7 15 × 155040×504015715×15
Douglas Zare
8

Acho que seu algoritmo simples embaralha as cartas corretamente, pois o número de embaralhamentos tende ao infinito.

Suponha que você tenha três cartas: {A, B, C}. Suponha que suas cartas comecem na seguinte ordem: A, B, C. Depois de um shuffle, você tem as seguintes combinações:

{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}

Portanto, a probabilidade de a carta A estar na posição {1,2,3} é {5/9, 2/9, 2/9}.

Se embaralharmos as cartas uma segunda vez, então:

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) 

Isso dá 0,407.

Usando a mesma idéia, podemos formar um relacionamento de recorrência, ou seja:

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

Codificar isso em R (veja o código abaixo) fornece a probabilidade de o cartão A estar na posição {1,2,3} como {0,33333, 0,333333, 0,3333} após dez shuffles.

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
csgillespie
fonte
1
+1. Isso demonstra que a probabilidade de um determinado cartão terminar em uma determinada posição se aproxima da proporção esperada à medida que o número de shuffles aumenta. No entanto, o mesmo também se aplica a um algoritmo que apenas gira o array uma vez de forma aleatória: todas as cartas têm uma probabilidade igual de terminar em todas as posições, mas ainda não há aleatoriedade (o array permanece classificado).
Thilo
@ Thilo: Desculpe, eu não sigo o seu comentário. Um "algoritmo gira de forma aleatória", mas ainda não existe "aleatoriedade"? Você poderia explicar mais?
csgillespie
Se você "embaralha" uma matriz de elementos N girando-a entre as posições 0 e N-1 (aleatoriamente), cada cartão tem exatamente a mesma probabilidade de terminar em qualquer uma das posições N, mas 2 ainda está sempre localizado entre 1 e 3.
Thilo
1
@ Thio: Ah, entendi seu ponto. Bem, você pode calcular a probabilidade (usando exatamente a mesma idéia acima), para o Pr (A na posição 2) e Pr (A na posição 3) - aqui para os cartões B e C. Você verá que todas as probabilidades tendem a 1/3. Nota: minha resposta fornece apenas um caso específico, enquanto que a @whuber nice response fornece o caso geral.
Csgillespie
4

Uma maneira de ver que você não terá uma distribuição perfeitamente uniforme é pela divisibilidade. Na distribuição uniforme, a probabilidade de cada permutação é de. Ao gerar uma sequência de transposições aleatórios, e em seguida, recolher por sequências seu produto, as probabilidades chegar são da forma para algum número inteiro . Se , então . Pelo Postulado de Bertrand (um teorema), para existem números primos que ocorrem no denominador e que não dividem , entãonão é um número inteiro e não há como dividir as transposições uniformemente emt A / n 2 t A 1 / n ! = A / n 2 t n 2 t / n ! = A n 3 n n 2 t / n ! n ! n = 52 1 / 52 ! 3 , 5 , 7 , . . . , 47 1 /1/n!tUMA/n2tUMA1/n!=UMA/n2tn2t/n!=UMAn3nn2t/n!n!permutações. Por exemplo, se , então o denominador deé divisível por enquanto o denominador de não é, portanto não pode ser reduzido para.n=521/52!3,5,7,...,471/522tUMA/522t1/52!

Quantos você precisa para aproximar bem uma permutação aleatória? A geração de uma permutação aleatória por transposições aleatórias foi analisada por Diaconis e Shahshahani usando a teoria da representação do grupo simétrico em

Diaconis, P., Shahshahani, M. (1981): "Gerando uma permutação aleatória com transposições aleatórias". Z. Wahrsch. Verw. Geb. 57, 159-179.

Uma conclusão foi que são necessárias transposições no sentido de que após as permutações estão longe de serem aleatórias, mas após o resultado é quase aleatório, tanto no sentido da variação total quanto da distância . Esse tipo de fenômeno de corte é comum em caminhadas aleatórias em grupos e está relacionado ao famoso resultado de que você precisa de embaralhamento de rifles antes que um baralho se torne quase aleatório.12nregistron(1-ϵ)12nregistron(1+ϵ)12nregistroneu27

Douglas Zare
fonte
2

Tenha em mente que eu não sou um estatístico, mas vou colocar meus 2 centavos.

Fiz um pequeno teste em R (cuidado, é muito lento para alto numTrials, o código provavelmente pode ser otimizado):

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

Isso gerará uma matriz swapscom numTrials+1linhas (uma por tentativa + o original) e numElementscolunas (uma por cada elemento do vetor). Se o método estiver correto, a distribuição de cada coluna (ou seja, dos valores de cada elemento nas tentativas) não deve ser diferente da distribuição dos dados originais.

Como nossos dados originais eram normalmente distribuídos, esperaríamos que todas as colunas não se desviassem disso.

Se corrermos

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

Nós temos:

Histogramas de ensaios aleatórios

o que parece muito promissor. Agora, se queremos confirmar estatisticamente que as distribuições não se desviam do original, acho que poderíamos usar um teste de Kolmogorov-Smirnov (por favor, algum estatístico pode confirmar que isso está certo?) E, por exemplo,

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

O que nos dá p = 0,9926

Se verificarmos todas as colunas:

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

E nós corremos

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

Nós temos:

Histograma dos valores de p de Kolmogorov-Smirnov

Portanto, para a grande maioria dos elementos da matriz, seu método de troca deu um bom resultado, como você também pode ver olhando os quartis.

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

Observe que, obviamente, com um número menor de tentativas, a situação não é tão boa:

50 tentativas

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

100 ensaios

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

500 ensaios

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

Aqui está como eu estou interpretando seu algoritmo, em pseudo-có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]

2×euength×nvocêm_pumasses[0 0,euength-1]euength

euength2×euength×nvocêm_pumasses

euength!euength!<euength2×euength×nvocêm_pumasses

euength!|euength2×euength×nvocêm_pumasses

pp<euengthpeuengtheuength>2p|euength!euength2×euength×nvocêm_pumasseseuength!euength2×euength×nvocêm_pumasseseuength>2

euengthp<euengtheuength-1euength-1euength

euengtheuength-1euength!euength!|euength!. Não é difícil mostrar que cada traço resulta em uma permutação diferente e, a partir daí, é fácil ver que Fisher-Yates gera cada permutação com igual probabilidade.

tzs
fonte