div - sidebarpanel shiny

R: comando de muestra() sujeto a una restricción

Question

Estoy tratando de muestrear aleatoriamente 7 números del 0 al 7 (con reemplazo), pero sujeto a la restricción de que los números elegidos sumen 7. Entonces, por ejemplo, la salida 0 1 1 2 3 0 0 está bien, pero la salida 1 2 3 4 5 6 7 no lo es. ¿Hay alguna manera de usar el comando de ejemplo con restricciones agregadas?

Tengo la intención de usar la función replicar () con el comando de ejemplo como un argumento, para devolver una lista de N vectores diferentes del comando de ejemplo. En la forma en que actualmente estoy usando el comando de muestra (sin ninguna restricción), necesito que N sea muy grande para obtener tantos vectores posibles que sumen exactamente 7 como sea posible. ¡Me imagino que debe haber una manera más fácil de hacer esto!

Aquí está mi código para esa parte:

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

Idealmente, quiero que 10,000 o 100,000 vectores en x sumen 7, pero necesitaría un enorme valor de N para hacer esto. Gracias por cualquier ayuda.

Answer 1

Comience con todos los ceros, agregue uno a cualquier elemento, haga 7 veces:

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

O de manera equivalente, simplemente elija cuál de los 7 elementos que va a incrementar en una muestra de longitud 7, luego tabúlelos, asegurándose de tabular hasta 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

No sé si esto producirá una muestra uniforme de todas las combinaciones posibles ...

La distribución de elementos individuales a lo largo de 100,000 repeticiones es:

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

No golpeó un 0,0,0,0,0,7 esa vez!

Answer 2

Encuentro esta pregunta intrigante y pensé un poco más. Otro enfoque (más general) para (aproximar) muestrear uniformemente de todas las soluciones factibles, sin generar y almacenar todas las permutaciones (lo que claramente no es posible en el caso con mucho más de 7 números), en R por sample() , podría ser un Implementación sencilla de MCMC:

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 

Esto es, por supuesto, muy rápido para algunas muestras. La distribución":

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

Por supuesto, en este caso, donde es realmente posible encontrar y almacenar todas las combinaciones, y si desea una gran muestra de todos los resultados posibles, simplemente use partitions::compositions(7, 7) , como también lo sugiere Josh O''Brien en Los comentarios, para evitar calcular todas las permutaciones, cuando solo se necesita una pequeña fracción:

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

Answer 3

Este algoritmo recursivo generará una distribución con mayor probabilidad de grandes números que las otras soluciones. La idea es lanzar un número aleatorio y en 0:7 en cualquiera de las siete ranuras disponibles, luego repetir con un número aleatorio en 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

Dibujar 100,000 vectores tomó 11 segundos en mi máquina y aquí está la distribución que obtengo:

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

Answer 4

Para asegurarte de que estás muestreando uniformemente, puedes generar todas las permutaciones y limitarlas a las que suman 7:

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

Desde nrow(perms7) , vemos que solo hay 1716 permutaciones posibles que suman 7. Ahora puede muestrear uniformemente de las permutaciones:

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

Una ventaja de este enfoque es que es fácil ver que estamos muestreando uniformemente al azar. Además, es bastante rápido: la creación de perms7 tomó 0.3 segundos en mi computadora y la creación de 1 millón de filas my.perms tomó 0.04 segundos. Si necesita dibujar muchos vectores, esto será un poco más rápido que un enfoque recursivo porque solo está usando la indexación matricial en perms7 en lugar de generar cada vector por separado.

Aquí hay una distribución de cuentas de números en la muestra:

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

Answer 5

Puede haber una forma más fácil y / o más elegante, pero aquí hay un método de fuerza bruta que usa la función LSPM:::.nPri . El enlace incluye la definición de una versión solo en R del algoritmo, para aquellos interesados.

#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))]

Ahora puede muestrear desde z estar seguro de que está obteniendo una permutación que suma 7.