¿Cómo remuestrear en R sin repetir las permutaciones?

Question

En R, si establezco.seed(), y luego utilizo la función sample para aleatorizar una lista, ¿puedo garantizar que no generaré la misma permutación?

es decir...

set.seed(25)
limit <- 3
myindex <- seq(0,limit)
for (x in seq(1,factorial(limit))) {
    permutations <- sample(myindex)
    print(permutations)
}

Esto produce

[1] 1 2 0 3
[1] 0 2 1 3
[1] 0 3 2 1
[1] 3 1 2 0
[1] 2 3 0 1
[1] 0 1 3 2

¿todas las permutaciones impresas serán permutaciones únicas? ¿O hay alguna posibilidad, basada en la forma en que está implementada, de que pueda obtener algunas repeticiones?

Quiero ser capaz de hacer esto sin repeticiones, garantizado. ¿Cómo podría hacerlo?

(También quiero evitar tener que usar una función como permn(), que tiene un método muy mecánico para generar todas las permutaciones--no parece aleatorio).

Además, nota al margen, parece que este problema es O((n!)), si no me equivoco.

Answer 1

La pregunta tiene muchas interpretaciones válidas. Los comentarios -especialmente el que indica que se necesitan permutaciones de 15 o más elementos (¡15! = 1307674368000 se hace grande)- sugieren que lo que se quiere es un relativamente pequeño muestra aleatoria, sin reemplazo, de todos los n! = n*(n-1) (n-2) ...*2*1 permutaciones de 1:n. Si esto es cierto, existen soluciones (algo) eficientes.

La siguiente función, rperm acepta dos argumentos n (el tamaño de las permutaciones a muestrear) y m (el número de permutaciones de tamaño n a dibujar). Si m se aproxima o supera a n!, la función tardará mucho tiempo y devolverá muchos valores NA: está pensada para usarse cuando n es relativamente grande (digamos, 8 o más) y m es mucho menor que n! Funciona almacenando en la memoria caché una representación de cadena de las permutaciones encontradas hasta el momento y luego generando nuevas permutaciones (al azar) hasta que se encuentre una nueva. Aprovecha la capacidad de indexación asociativa de R para buscar rápidamente en la lista de permutaciones encontradas anteriormente.

rperm <- function(m, size=2) { # Obtain m unique permutations of 1:size

    # Function to obtain a new permutation.
    newperm <- function() {
        count <- 0                # Protects against infinite loops
        repeat {
            # Generate a permutation and check against previous ones.
            p <- sample(1:size)
            hash.p <- paste(p, collapse="")
            if (is.null(cache[[hash.p]])) break

            # Prepare to try again.
            count <- count+1
            if (count > 1000) {   # 1000 is arbitrary; adjust to taste
                p <- NA           # NA indicates a new permutation wasn't found
                hash.p <- ""
                break
            }
        }
        cache[[hash.p]] <<- TRUE  # Update the list of permutations found
        p                         # Return this (new) permutation
    }

    # Obtain m unique permutations.
    cache <- list()
    replicate(m, newperm())  
} # Returns a `size` by `m` matrix; each column is a permutation of 1:size.

La naturaleza de replicate es devolver las permutaciones como columna vectores; Por ejemplo A continuación se reproduce un ejemplo de la pregunta original, transpuesta :

> set.seed(17)
> rperm(6, size=4)
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    2    4    4    3    4
[2,]    3    4    1    3    1    2
[3,]    4    1    3    2    2    3
[4,]    2    3    2    1    4    1

Los tiempos son excelentes para valores pequeños y moderados de m, hasta unos 10.000, pero se degradan para problemas mayores. Por ejemplo, una muestra de m = 10.000 permutaciones de n = 1000 elementos (una matriz de 10 millones de valores) se obtuvo en 10 segundos; una muestra de m = 20.000 permutaciones de n = 20 elementos requirió 11 segundos, a pesar de que la salida (una matriz de 400.000 entradas) era mucho más pequeña; y el cálculo de la muestra de m = 100.000 permutaciones de n = 20 elementos se abortó después de 260 segundos (no tuve paciencia para esperar a que se completara). Este problema de escala parece estar relacionado con las ineficiencias de escala en el direccionamiento asociativo de R. Se puede solucionar generando muestras en grupos de, digamos, 1000 o así, y luego combinando esas muestras en una muestra grande y eliminando los duplicados. Los expertos en R podrían sugerir soluciones más eficientes o mejores soluciones.

Editar

Podemos lograr un rendimiento asintótico casi lineal dividiendo la caché en una jerarquía de dos cachés, para que R nunca tenga que buscar en una lista grande. Conceptualmente (aunque no como se ha implementado), crea un array indexado por el primer $k$ elementos de una permutación. Las entradas de esta matriz son listas de todas las permutaciones que comparten esos primeros $k$ elementos. Para comprobar si una permutación ha sido vista, utilice su primer $k$ para encontrar su entrada en la caché y luego buscar esa permutación dentro de esa entrada. Podemos elegir $k$ para equilibrar los tamaños previstos de todas las listas. La implementación real no utiliza un $k$ -fold array, lo que sería difícil de programar con suficiente generalidad, pero en su lugar utiliza otra lista.

A continuación se muestran algunos tiempos transcurridos en segundos para un rango de tamaños de permutación y números de permutaciones distintas solicitadas:

 Number Size=10 Size=15 Size=1000 size=10000 size=100000
     10    0.00    0.00      0.02       0.08        1.03
    100    0.01    0.01      0.07       0.64        8.36
   1000    0.08    0.09      0.68       6.38
  10000    0.83    0.87      7.04      65.74
 100000   11.77   10.51     69.33
1000000  195.5   125.5

(El aumento de velocidad aparentemente anómalo del tamaño=10 al tamaño=15 se debe a que el primer nivel de la caché es mayor para el tamaño=15, lo que reduce el número medio de entradas en las listas de segundo nivel, acelerando así la búsqueda asociativa de R. Con un cierto coste de RAM, la ejecución podría ser más rápida aumentando el tamaño de la caché de nivel superior. Sólo con aumentar k.head por 1 (que multiplica el tamaño del nivel superior por 10) aceleró rperm(100000, size=10) de 11,77 segundos a 8,72 segundos, por ejemplo. Al hacer la caché de nivel superior 10 veces más grande no se consiguió una ganancia apreciable, con 8,51 segundos).

Salvo en el caso de 1.000.000 de permutaciones únicas de 10 elementos (¡una parte sustancial de todas las 10! = unos 3,63 millones de permutaciones de este tipo), no se detectó prácticamente ninguna colisión. En este caso excepcional, se produjeron 169.301 colisiones, pero no se produjeron fallos completos (de hecho, se obtuvo un millón de permutaciones únicas).

Obsérvese que con tamaños de permutación grandes (superiores a 20, aproximadamente), la probabilidad de obtener dos permutaciones idénticas, incluso en una muestra tan grande como 1.000.000.000, es prácticamente nula. Por lo tanto, esta solución es aplicable principalmente en situaciones en las que (a) hay un gran número de permutaciones únicas de (b) entre $n=5$ y $n=15$ de elementos, pero aún así, (c) un número sustancialmente menor de todos los $n!$ se necesitan permutaciones.

El código de trabajo es el siguiente.

rperm <- function(m, size=2) { # Obtain m unique permutations of 1:size
    max.failures <- 10

    # Function to index into the upper-level cache.
    prefix <- function(p, k) {    # p is a permutation, k is the prefix size
        sum((p[1:k] - 1) * (size ^ ((1:k)-1))) + 1
    } # Returns a value from 1 through size^k

    # Function to obtain a new permutation.
    newperm <- function() {
        # References cache, k.head, and failures in parent context.
        # Modifies cache and failures.        

        count <- 0                # Protects against infinite loops
        repeat {
            # Generate a permutation and check against previous ones.
            p <- sample(1:size)
            k <- prefix(p, k.head)
            ip <- cache[[k]]
            hash.p <- paste(tail(p,-k.head), collapse="")
            if (is.null(ip[[hash.p]])) break

            # Prepare to try again.
            n.failures <<- n.failures + 1
            count <- count+1
            if (count > max.failures) {  
                p <- NA           # NA indicates a new permutation wasn't found
                hash.p <- ""
                break
            }
        }
        if (count <= max.failures) {
            ip[[hash.p]] <- TRUE      # Update the list of permutations found
            cache[[k]] <<- ip
        }
        p                         # Return this (new) permutation
    }

    # Initialize the cache.
    k.head <- min(size-1, max(1, floor(log(m / log(m)) / log(size))))
    cache <- as.list(1:(size^k.head))
    for (i in 1:(size^k.head)) cache[[i]] <- list()

    # Count failures (for benchmarking and error checking).
    n.failures <- 0

    # Obtain (up to) m unique permutations.
    s <- replicate(m, newperm())
    s[is.na(s)] <- NULL
    list(failures=n.failures, sample=matrix(unlist(s), ncol=size))
} # Returns an m by size matrix; each row is a permutation of 1:size.

Answer 2

Utilizando unique de la manera correcta debería funcionar:

set.seed(2)
limit <- 3
myindex <- seq(0,limit)

endDim<-factorial(limit)
permutations<-sample(myindex)

while(is.null(dim(unique(permutations))) || dim(unique(permutations))[1]!=endDim) {
    permutations <- rbind(permutations,sample(myindex))
}
# Resulting permutations:
unique(permutations)

# Compare to
set.seed(2)
permutations<-sample(myindex)
for(i in 1:endDim)
{
permutations<-rbind(permutations,sample(myindex))
}
permutations
# which contains the same permutation twice

Answer 3

Voy a desviarme un poco de tu primera pregunta, y sugerir que si estás tratando con vectores relativamente cortos, podrías simplemente generar todas las permutaciones usando permn y los ordenan aleatoriamente esos utilizando sample :

x <- combinat:::permn(1:3)
> x[sample(factorial(3),factorial(3),replace = FALSE)]
[[1]]
[1] 1 2 3

[[2]]
[1] 3 2 1

[[3]]
[1] 3 1 2

[[4]]
[1] 2 1 3

[[5]]
[1] 2 3 1

[[6]]
[1] 1 3 2