pbinom - poisson distribution in r

Comportamiento errático de semillas con rbinom(prob=0.5) en R

Question

He encontrado lo que consideraría una conducta errática (pero espero que haya una explicación simple) en el uso de semillas de R junto con rbinom() cuando se usa prob=0.5 . Idea general: para mí, si configuro la semilla, ejecute rbinom() una vez (es decir, lleve a cabo un único proceso aleatorio), a pesar de que el valor prob esté establecido en, la semilla aleatoria debe cambiar en un incremento. Luego, si vuelvo a establecer la semilla en el mismo valor, y ejecuto otro proceso aleatorio (como rbinom() otra vez, pero tal vez con un valor diferente de prob ), la semilla debería volver a cambiar al mismo valor que lo hizo para el proceso individual previo aleatorio.

He encontrado que R hace exactamente esto siempre que use rbinom() con cualquier prob!=0.5 . Aquí hay un ejemplo:

Compara el vector de semillas, .Random.seed , para dos probabilidades distintas de 0.5:

set.seed(234908)
x <- rbinom(n=1,size=60,prob=0.4)
temp1 <- .Random.seed
set.seed(234908)
x <- rbinom(n=1,size=60,prob=0.3)
temp2 <- .Random.seed
any(temp1!=temp2)
> [1] FALSE

Compare el vector de semillas, .Random.seed , para prob = 0.5 vs. prob! = 0.5:

set.seed(234908)
x <- rbinom(n=1,size=60,prob=0.5)
temp1 <- .Random.seed
set.seed(234908)
x <- rbinom(n=1,size=60,prob=0.3)
temp2 <- .Random.seed
any(temp1!=temp2)
> [1] TRUE
temp1==temp2
> [1] TRUE FALSE TRUE TRUE TRUE TRUE TRUE
> [8] TRUE TRUE TRUE TRUE TRUE TRUE TRUE
...

He encontrado esto para todas las comparaciones de prob=0.5 contra todas las otras probabilidades en el conjunto {0.1, 0.2, ..., 0.9}. De manera similar, si comparo cualquier valor de prob de {0.1, 0.2, ..., 0.9} distinto de 0.5, el vector .Random.seed siempre es igual a elemento por elemento. Estos hechos también son válidos para el size impar o par dentro de rbinom() .

Para hacerlo aún más extraño (pido disculpas por el hecho de que es un poco intrincado, es relevante para la forma en que se escribe mi función), cuando uso probabilidades guardadas como elementos en un vector, tengo el mismo problema si 0.5 es primer elemento, pero no segundo. Aquí está el ejemplo para este caso:

Primer caso: 0.5 es la primera probabilidad referenciada en el vector

set.seed(234908)
MNAR <- c(0.5,0.3)
x <- rbinom(n=1,size=60,prob=MNAR[1])
y <- rbinom(n=1,size=50,prob=MNAR[2])
temp1 <- .Random.seed
set.seed(234908)
MNAR <- c(0.1,0.3)
x <- rbinom(n=1,size=60,prob=MNAR[1])
y <- rbinom(n=1,size=50,prob=MNAR[2])
temp2 <- .Random.seed
any(temp1!=temp2)
> [1] TRUE
any(temp1!=temp2)
> [1] TRUE FALSE TRUE TRUE TRUE TRUE TRUE
> [8] TRUE TRUE TRUE TRUE TRUE TRUE TRUE

Segundo caso: 0.5 es la segunda probabilidad referenciada en el vector

set.seed(234908)
MNAR <- c(0.3,0.5)
x <- rbinom(n=1,size=60,prob=MNAR[1])
y <- rbinom(n=1,size=50,prob=MNAR[2])
temp1 <- .Random.seed
set.seed(234908)
MNAR <- c(0.1,0.3)
x <- rbinom(n=1,size=60,prob=MNAR[1])
y <- rbinom(n=1,size=50,prob=MNAR[2])
temp2 <- .Random.seed
any(temp1!=temp2)
> [1] FALSE

Una vez más, me parece que a pesar de los valores utilizados para el prob y el size , este patrón se cumple. ¿Alguien puede explicarme este misterio? Está causando un gran problema porque los resultados que deberían ser los mismos son diferentes porque la semilla se usa / calcula de manera diferente cuando prob=0.5 pero en ninguna otra instancia.

Answer 1

Entonces, convirtamos nuestros comentarios en una respuesta. Gracias a Ben Bolker por ponernos en el camino correcto con un enlace al código: svn.r-project.org/R/trunk/src/nmath/rbinom.c y la sugerencia de encontrar dónde unif_rand() se llama.

Un análisis rápido y parece que el código está dividido en dos secciones, delimitadas por los comentarios:

/*-------------------------- np = n*p >= 30 : ------------------- */

y

/*---------------------- np = n*p < 30 : ------------------------- */

Dentro de cada uno de estos, el número de llamadas a unif_rand no es el mismo (dos contra uno).

Entonces para un size dado ( n ), su semilla aleatoria puede terminar en un estado diferente dependiendo del valor de prob ( p ): si el size * prob >= 30 o no.

Con eso en mente, todos los resultados que obtuvo con sus ejemplos ahora deberían tener sentido:

# these end up in the same state
rbinom(n=1,size=60,prob=0.4) # => np < 30
rbinom(n=1,size=60,prob=0.3) # => np < 30
# these don''t
rbinom(n=1,size=60,prob=0.5) # => np >= 30
rbinom(n=1,size=60,prob=0.3) # => np < 30
# these don''t
{rbinom(n=1,size=60,prob=0.5) # np >= 30
rbinom(n=1,size=50,prob=0.3)} # np < 30
{rbinom(n=1,size=60,prob=0.1) # np < 30
rbinom(n=1,size=50,prob=0.3)} # np < 30
# these do
{rbinom(n=1,size=60,prob=0.3) # np < 30
rbinom(n=1,size=50,prob=0.5)} # np < 30
{rbinom(n=1,size=60,prob=0.1) # np < 30
rbinom(n=1,size=50,prob=0.3)} # np < 30

Answer 2

Voy a tomar una posición contraria con respecto a esta pregunta y afirmar que las expectativas no son apropiadas y no están respaldadas por la documentación. La documentación no establece ningún reclamo sobre qué efectos secundarios (específicamente en .Random.seed ) se pueden esperar al llamar a rbinom , o cómo esos efectos secundarios pueden ser o no iguales en varios casos.

rbinom tiene tres parámetros: n , size y prob . Su expectativa es que, para un conjunto de semillas al azar antes de llamar a rbinom , .Random.seed será el mismo después de llamar a rbinom para un n dado y cualquier valor de size y prob (o tal vez cualquier valor finito de size y prob ). Ciertamente te das cuenta de que sería diferente para diferentes valores de n . rbinom no garantiza eso ni implica eso.

Sin conocer los aspectos internos de la función, esto no puede ser conocido; como mostró la otra respuesta , el algoritmo es diferente en función del producto de size y prob . Y las partes internas pueden cambiar, por lo que estos detalles específicos pueden cambiar.

Al menos, en este caso, la .Random.seed resultante será la misma después de cada llamada de rbinom que tenga el mismo n , size y prob . Puedo construir una función patológica para la cual esto ni siquiera es verdad:

seedtweak <- function() {
if(floor(as.POSIXlt(Sys.time())$sec * 10) %% 2) {
runif(1)
}
invisible(NULL)
}

Básicamente, esta función busca si las décimas del segundo momento son par o impar para decidir si se dibuja un número aleatorio o no. Ejecute esta función y .Random.seed puede o no cambiar:

rs <- replicate(10, {
set.seed(123)
seedtweak()
.Random.seed
})
all(apply(rs, 1, function(x) Reduce(`==`, x)))

Lo mejor que puede (¿debería?) Esperar es que un conjunto dado de códigos con todas las entradas / parámetros de la misma (incluida la semilla) siempre dé resultados idénticos. Esperar resultados idénticos cuando la mayoría (o solo algunos) de los parámetros son iguales no es realista a menos que todas las funciones llamadas hagan esas garantías.