R: Komenda sample() podlegająca ograniczeniu

Question

Próbuję losowo próbkować 7 liczb od 0 do 7 (z zamianą), ale z zastrzeżeniem, że wybrane liczby sumują się do 7. Tak na przykład, wyjście 0 1 1 2 3 0 0 jest w porządku, ale wynik 1 2 3 4 5 6 7 nie jest. Czy istnieje sposób użycia przykładowego polecenia z dodanymi ograniczeniami?

Zamierzam użyć funkcji replicate() z poleceniem sample jako argumentem, aby zwrócić listę N różnych wektorów z przykładowego polecenia. Sposób, w jaki aktualnie używam przykładowego polecenia (bez żadnych ograniczeń), wymaga, aby N był bardzo duży, aby uzyskać jak najwięcej możliwych wektorów, które sumują się do dokładnie 7, jak to możliwe. Sądzę, że musi być łatwiejszy sposób na zrobienie tego!

Oto mój kod dla tej strony:

x <- replicate(100000, sample(0:7, 7, replace=T))  

Idealnie chcę 10000 100000 lub wektory w X podsumować do 7, ale musiałby ogromną wartość N, aby to zrobić. Dzięki za pomoc.

Answer 1

aby upewnić się, że pobieranie próbek równomiernie, można po prostu generowania wszystkich permutacji i ograniczeń co do tych, które suma 7:

library(gtools) 
perms <- permutations(8, 7, 0:7, repeats.allowed=T) 
perms7 <- perms[rowSums(perms) == 7,] 

Od nrow(perms7) widzimy istnieją tylko 1716 możliwych permutacji tej sumy do 7. Teraz można równomiernie próbka z permutacji:

set.seed(144) 
my.perms <- perms7[sample(nrow(perms7), 100000, replace=T),] 
head(my.perms) 
#  [,1] [,2] [,3] [,4] [,5] [,6] [,7] 
# [1,] 0 0 0 2 5 0 0 
# [2,] 1 3 0 1 2 0 0 
# [3,] 1 4 1 1 0 0 0 
# [4,] 1 0 0 3 0 3 0 
# [5,] 0 2 0 0 0 5 0 
# [6,] 1 1 2 0 0 2 1 

Zaletą tego podejścia jest to, że łatwo zrozumieć, że jesteśmy próbkowania równomiernie losowo. Ponadto, jest to dość szybkie - tworzenie perms7 trwało 0,3 sekundy na moim komputerze, a zbudowanie 1 miliona wierszy my.perms trwało 0,04 sekundy. Jeśli chcesz narysować wiele wektorów, będzie to trochę szybciej niż podejście rekursywne, ponieważ używasz indeksowania macierzy do perms7 zamiast generowania każdego wektora osobno.

Oto rozkład liczby numerów w próbce:

#  0  1  2  3  4  5  6  7 
# 323347 188162 102812 51344 22811 8629 2472 423 

Answer 2

start ze wszystkimi zerami, dodać jeden do dowolnego elementu, zrobić 7 razy:

sumTo = function(){ 
    v = rep(0,7) 
    for(i in 1:7){ 
     addTo=sample(7)[1] 
     v[addTo]=v[addTo]+1 
    } 
    v 
} 

albo równoważnie, wystarczy wybrać, który z 7 elementów masz zamiar zwiększyć w jednej próbce o długości 7, następnie tabularyzować tych, upewniając się tabularyzować do 7:

sumTo = function(){tabulate(sample(7, 7, replace = TRUE), 7)} 


> sumTo() 
[1] 2 1 0 0 4 0 0 
> sumTo() 
[1] 1 3 1 0 1 0 1 
> sumTo() 
[1] 1 1 0 2 1 0 2 

nie wiem, czy to będzie uzyskania jednorodnej próbki ze wszystkich możliwych kombinacji ...

dystrybucja ind Pozostałe elementy powyżej 100 000 powtórzeń to:

> X = replicate(100000,sumTo()) 
> table(X) 
X 
    0  1  2  3  4  5  6 
237709 277926 138810 38465 6427 627  36 

Nie trafił w tym czasie 0,0,0,0,0,7!

Answer 3

Ten algorytm rekurencyjny wyprowadzi dystrybucję z większym prawdopodobieństwem dla dużych liczb niż inne rozwiązania. Chodzi o to, aby rzucić losową liczbę y w 0:7 w każdym z siedmiu dostępnych szczelinach, a następnie powtórz z liczbą losową w 0:(7-y), etc:

sample.sum <- function(x = 0:7, n = 7L, s = 7L) { 
    if (n == 1) return(s) 
    x <- x[x <= s] 
    y <- sample(x, 1) 
    sample(c(y, Recall(x, n - 1L, s - y))) 
} 

set.seed(123L) 
sample.sum() 
# [1] 0 4 0 2 0 0 1 

Rysunek 100,000 wektorów trwało 11 sekund na moim komputerze i tu jest dystrybucja uzyskać:

#  0  1  2  3  4  5  6  7 
# 441607 98359 50587 33364 25055 20257 16527 14244 

Answer 4

Nie może być łatwiej i/lub bardziej elegancki sposób, ale tutaj jest to metoda brute-force przy użyciu funkcji LSPM:::.nPri. Łącze zawiera definicję algorytmu tylko dla R, dla zainteresowanych.

#install.packages("LSPM", repos="http://r-forge.r-project.org") 
library(LSPM) 
# generate all possible permutations, since there are only ~2.1e6 of them 
# (this takes < 40s on my 2.2Ghz laptop) 
x <- lapply(seq_len(8^7), nPri, n=8, r=7, replace=TRUE) 
# set each permutation that doesn't sum to 7 to NULL 
y <- lapply(x, function(p) if(sum(p-1) != 7) NULL else p-1) 
# subset all non-NULL permutations 
z <- y[which(!sapply(y, is.null))] 

Teraz można spróbować z z i mieć pewność, że dostajesz permutacji że sum 7.

Answer 5

znajdę to pytanie intrygujące i dał mu kilka dodatkowych myśli. Inne (bardziej ogólne) podejście do (przybliżonej) próbki równomiernie od wszystkich możliwych rozwiązań, bez generowania i przechowywania wszystkich permutacji (co jest oczywiście niemożliwe w przypadku znacznie więcej niż 7 liczb), w R przez sample(), może być prostym MCMC implementacja:

S <- c(0, 1, 1, 2, 3, 0, 0) #initial solution 
N <- 100 #number of dependent samples (or burn in period) 
series <- numeric(N) 
for(i in 1:N){ 
    b <- sample(1:length(S), 2, replace=FALSE) #pick 2 elements at random 
    opt <- sum(S[-b]) #sum of complementary elements 
    a <- sample(0:(7-opt), 1) #sample a substistute 
    S[b[1]] <- a #change elements 
    S[b[2]] <- 7 - opt - a 
} 
S #new sample 

To oczywiście bardzo szybko dla kilku próbek. W „dystrybucja”:

#"distribution" N=100.000:  0  1  2  3  4  5  6  7 
#       321729 189647 103206 52129 22287 8038 2532 432 

Oczywiście w tym przypadku, gdy jest to rzeczywiście możliwe, aby znaleźć i zapisać wszystkie kombinacje, a jeśli chcesz ogromny próbki ze wszystkich możliwych wyników, wystarczy użyć partitions::compositions(7, 7), jak również sugerowane przez Josh O'Brien w komentarzach, aby uniknąć obliczania wszystkich permutacji, gdy potrzebna jest tylko niewielka część:

perms7 <- partitions::compositions(7, 7) 

>tabulate(perms7[, sample(ncol(perms7), 100000, TRUE)]+1, 8) 
#"distribution" N=100.000:  0  1  2  3  4  5  6  7 
#       323075 188787 102328 51511 22754 8697 2413 435