Jak ponownie próbkować w R bez powtarzania permutacji?

12

Czy w R, jeśli ustawię set.seed (), a następnie użyję przykładowej funkcji do losowej listy, czy mogę zagwarantować, że nie wygeneruję tej samej permutacji?

to znaczy...

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

To produkuje

[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

czy wszystkie wydrukowane permutacje będą unikatowymi permutacjami? Czy może jest jakaś szansa, w oparciu o sposób, w jaki jest to realizowane, że mogę dostać kilka powtórzeń?

Gwarantuję, że mogę to zrobić bez powtórzeń. Jak mam to zrobić?

(Chcę również uniknąć konieczności używania funkcji takiej jak permn (), która ma bardzo mechanistyczną metodę generowania wszystkich permutacji --- nie wygląda losowo.)

Ponadto, sidenote --- wygląda na to, że ten problem to O ((n!)!), Jeśli się nie mylę.

Mittenchops
źródło
Domyślnie argument „replace” dla „sample” jest ustawiony na FALSE.
ocram
Dzięki ocram, ale to działa w obrębie konkretnej próbki. To gwarantuje, że 0,1,2 i 3 nie powtórzą się podczas losowania (więc nie mogę losować 0,1,2,2), ale nie wiem, czy to gwarantuje, że druga próbka, Nie mogę ponownie narysować tej samej sekwencji 0123. Właśnie to zastanawiam się pod względem implementacji, czy ustawienie zarodka ma jakikolwiek wpływ na to powtórzenie.
Mittenchops,
Tak, to w końcu zrozumiałem, czytając odpowiedzi ;-)
ocram
1
Jeśli limitprzekroczy 12, prawdopodobnie zabraknie pamięci RAM, gdy R spróbuje przydzielić miejsce seq(1,factorial(limit)). (12! Wymaga około 2 GB, więc 13! Potrzebuje około 25 GB, 14! Około 350 GB itp.)
whuber
2
Istnieje szybkie, kompaktowe i eleganckie rozwiązanie do generowania losowych sekwencji wszystkich permutacji 1: n, pod warunkiem, że możesz wygodnie przechowywać n! liczby całkowite z zakresu 0: (n!). Łączy reprezentację permutacji w tabeli inwersji z silną reprezentacją liczbową.
whuber

Odpowiedzi:

9

Pytanie ma wiele ważnych interpretacji. Komentarze - szczególnie te wskazujące na konieczność permutacji 15 lub więcej elementów (15! = 1307674368000 stają się duże) - sugerują, że to, czego potrzebujemy, to stosunkowo niewielka losowa próbka, bez zamiany, wszystkich n! = n * (n-1) (n-2) ... * 2 * 1 permutacje 1: n. Jeśli to prawda, istnieją (nieco) wydajne rozwiązania.

Poniższa funkcja rpermakceptuje dwa argumenty n(rozmiar permutacji do próbki) i m(liczba permutacji o rozmiarze n do narysowania). Jeśli m zbliży się lub przekroczy n !, funkcja zajmie dużo czasu i zwróci wiele wartości NA: jest przeznaczona do użycia, gdy n jest względnie duże (powiedzmy 8 lub więcej), a m jest znacznie mniejsze niż n !. Działa poprzez buforowanie reprezentacji łańcuchowej znalezionych dotychczas permutacji, a następnie generowanie nowych permutacji (losowo), aż do znalezienia nowej. Wykorzystuje zdolność asocjatywnego indeksowania list R do szybkiego przeszukiwania listy znalezionych wcześniej permutacji.

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.

Naturą replicatejest zwracanie permutacji jako wektorów kolumnowych ; np . następujący tekst odtwarza przykład z pierwotnego pytania, transponowanego :

> 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

Czasy są doskonałe dla małych do umiarkowanych wartości m, do około 10 000, ale zmniejszają się w przypadku większych problemów. Na przykład próbkę m = 10 000 permutacji n = 1000 elementów (macierz 10 milionów wartości) uzyskano w 10 sekund; próbka m = 20 000 permutacji n = 20 elementów wymagała 11 sekund, mimo że wynik (macierz 400 000 wpisów) był znacznie mniejszy; a próbka obliczeniowa m = 100 000 permutacji n = 20 elementów została przerwana po 260 sekundach (nie miałem czasu czekać na zakończenie). Ten problem skalowania wydaje się być związany z nieefektywnością skalowania w adresowaniu asocjacyjnym R. Można obejść ten problem, generując próbki w grupach, powiedzmy, około 1000, a następnie łącząc je w dużą próbkę i usuwając duplikaty.

Edytować

kkkkk-fold tablica, która byłaby trudna do zaprogramowania w wystarczającej ogólności, ale zamiast tego używa innej listy.

Oto kilka upływających czasów w sekundach dla zakresu rozmiarów permutacji i liczby różnych żądanych kombinacji:

 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

(Pozornie nietypowe przyspieszenie od rozmiaru = 10 do rozmiaru = 15 jest spowodowane tym, że pierwszy poziom pamięci podręcznej jest większy dla rozmiaru = 15, co zmniejsza średnią liczbę wpisów na listach drugiego poziomu, przyspieszając w ten sposób wyszukiwanie asocjacyjne R. koszt pamięci RAM, wykonanie można przyspieszyć, zwiększając rozmiar pamięci podręcznej wyższego poziomu. Po prostu zwiększenie k.heado 1 (co zwielokrotnia rozmiar górnego poziomu przez 10) przyspieszyło np. rperm(100000, size=10)z 11,77 sekundy do 8,72 sekundy. pamięć podręczna 10-krotnie większa, ale nie osiągnęła żadnego znaczącego wzmocnienia, taktując 8,51 sekundy.)

Z wyjątkiem przypadku 1 000 000 unikalnych permutacji 10 elementów (znaczna część wszystkich 10! = Około 3,63 miliona takich permutacji) praktycznie nigdy nie wykryto kolizji. W tym wyjątkowym przypadku doszło do 169 301 kolizji, ale nie doszło do całkowitej awarii (faktycznie uzyskano milion unikalnych kombinacji).

n=5n=15n!

Następuje kod roboczy.

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.
Whuber
źródło
Jest blisko, ale zauważam, że dostaję błędy, takie jak 1, 2 i 4, ale myślę, że rozumiem, co masz na myśli i powinienem być w stanie z tym pracować. Dzięki! > rperm(6,3) $failures [1] 9 $sample [,1] [,2] [,3] [1,] 3 1 3 [2,] 2 2 1 [3,] 1 3 2 [4,] 1 2 2 [5,] 3 3 1 [6,] 2 1 3
Mittenchops,
3

uniqueWłaściwe użycie powinno załatwić sprawę:

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
MånsT
źródło
Przepraszamy za nieprawidłowe wyjaśnienie kodu. Jestem teraz w pośpiechu, ale chętnie odpowiem na pytania, które masz później. Poza tym nie mam pojęcia o szybkości powyższego kodu ...
MånsT
1
Sfunkcjonalizowałem to, co mi dałeś w ten sposób: `myperm <- function (limit) {myindex <- seq (0, limit) endDim <-factorial (limit) permutations <-sample (myindex) while (is.null (dim (unikalny (permutacje))) || dim (unique (permutations)) [1]! = endDim) {permutations <- rbind (permutations, sample (myindex))} return (unique (permutations))}} Działa, ale podczas gdy ja mogę zrobić limit = 6, limit = 7 powoduje przegrzanie mojego komputera. = PI uważa, że ​​musi istnieć sposób na
podpróbowanie
@Mittenchops, dlaczego mówicie, że musimy używać unikatowego do ponownego próbkowania w R bez powtarzania permutacji? Dziękuję Ci.
Frank
2

Przejdę trochę do twojego pierwszego pytania i zasugeruję, że jeśli masz do czynienia ze stosunkowo krótkimi wektorami, możesz po prostu wygenerować wszystkie permutacje za pomocą permni losowo uporządkować te za pomocą 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
joran
źródło
Bardzo mi się to podoba i jestem pewien, że to właściwe myślenie. Ale mój problem polega na tym, że używam sekwencji zwiększającej się do 10. Permn () był znacznie wolniejszy między silnią (7) i silnią (8), więc myślę, że 9 i 10 będą zbyt duże.
Mittenchops,
@Mittenchops Prawda, ale nadal jest możliwe, że naprawdę musisz je obliczyć tylko raz, prawda? Zapisz je do pliku, a następnie załaduj je, gdy ich potrzebujesz, i „próbkuj” z wcześniej zdefiniowanej listy. Możesz więc wykonać powolne obliczenia permn(10)lub cokolwiek innego tylko raz.
joran
Zgadza się, ale jeśli przechowuję gdzieś wszystkie permutacje, nawet to rozkłada się wokół silni (15) --- po prostu za dużo miejsca do przechowywania. Dlatego zastanawiam się, czy ustawienie zarodka pozwoli mi na zbiorowe próbkowanie permutacji --- a jeśli nie, czy istnieje algorytm do tego.
Mittenchops
@Mittenchops Ustawienie zarodka nie wpłynie na wydajność, gwarantuje tylko taki sam start za każdym razem, gdy zadzwonisz do PRNG.
Roman Luštrik
1
@Mitten Zobacz pomoc dla set.seed: opisuje, jak zapisać stan RNG i przywrócić go później.
whuber