Comportement erratique des graines avec rbinom(prob=0.5) dans R

Question

J'ai trouvé ce que je considère comme un comportement erratique (mais pour lequel j'espère qu'il y a une explication simple) en R L'utilisation des semences par l'UE en conjonction avec rbinom() quand prob=0.5 est utilisé. Idée générale : Pour moi, si je mets la graine, exécute rbinom() une fois (c'est-à-dire mener un seul processus aléatoire), malgré la valeur prob est fixé à, la valeur aléatoire devrait changer d'un incrément. Ensuite, si je fixe à nouveau la graine à la même valeur, et que j'exécute un autre processus aléatoire (tel que rbinom() à nouveau, mais peut-être avec une valeur différente de prob ), la graine doit à nouveau prendre la même valeur que pour le processus aléatoire unique précédent.

J'ai trouvé R fait exactement cela tant que j'utilise rbinom() avec n'importe quel prob!=0.5 . Voici un exemple :

Comparer le vecteur de semences, .Random.seed pour deux probabilités différentes 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

Comparer le vecteur de semences, .Random.seed pour prob=0,5 et 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
...

J'ai trouvé cela pour toutes les comparaisons de prob=0.5 contre toutes les autres probabilités dans l'ensemble {0,1, 0,2, ..., 0,9}. De même, si je compare n'importe quelle valeur de prob de {0,1, 0,2, ..., 0,9} autre que 0,5, l'indice .Random.seed est toujours égale élément par élément. Ces faits sont également valables pour les vecteurs pairs ou impairs. size sur rbinom() .

Pour rendre les choses encore plus étranges (je m'excuse si c'est un peu alambiqué - c'est lié à la façon dont ma fonction est écrite), lorsque j'utilise des probabilités enregistrées comme éléments d'un vecteur, j'ai le même problème si 0,5 est le premier élément, mais pas le second. Voici l'exemple pour ce cas :

Premier cas : 0.5 est la première probabilité référencée dans le vecteur

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

Deuxième cas : 0,5 est la deuxième probabilité référencée dans le vecteur

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

Encore une fois, je trouve que malgré les valeurs utilisées pour prob y size cette tendance se maintient. Quelqu'un peut-il m'expliquer ce mystère ? Cela pose un problème car les résultats qui devraient être identiques sont différents parce que la graine est, pour une raison ou une autre, utilisée/calculée différemment lorsque l'on utilise la méthode de la graine. prob=0.5 mais dans aucun autre cas.

Answer 1

Transformons donc nos commentaires en une réponse. Merci à Ben Bolker de nous avoir mis sur la bonne voie avec un lien vers le code : https://svn.r-project.org/R/trunk/src/nmath/rbinom.c et la suggestion de trouver où unif_rand() s'appelle.

Un examen rapide et il semble que le code soit divisé en deux sections, délimitées par les commentaires :

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

et

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

À l'intérieur de chacun d'entre eux, le nombre d'appels à unif_rand n'est pas le même (deux contre un.)

Ainsi, pour un size ( n ), votre graine aléatoire peut se retrouver dans un état différent selon la valeur de prob ( p ) : si size * prob >= 30 ou pas.

Avec cela en tête, tous les résultats que vous avez obtenus avec vos exemples devraient maintenant avoir un sens :

# 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

Je vais adopter une position contraire sur cette question et prétendre que les attentes ne sont pas appropriées et ne sont pas soutenues par la documentation. La documentation ne fait aucune déclaration sur les effets secondaires (en particulier sur la santé). .Random.seed ) peut être attendu en appelant rbinom ou comment ces effets secondaires peuvent ou non être les mêmes dans différents cas.

rbinom a trois paramètres : n , size y prob . Vous vous attendez à ce que, pour une graine aléatoire définie avant d'appeler rbinom , .Random.seed sera le même après avoir appelé rbinom pour un n y tout valeurs de size y prob (ou peut-être toute valeur finie de size y prob ). Vous réalisez certainement que ce serait différent pour différentes valeurs de n . rbinom ne le garantit pas ou ne l'implique pas.

Sans connaître les internes de la fonction, on ne peut pas le savoir ; comme le autre réponse a montré, l'algorithme est différent selon le produit de size y prob . Et les internes peuvent changer, donc ces détails spécifiques peuvent changer.

Au moins, dans ce cas, le résultat .Random.seed sera le même après chaque appel de rbinom qui a le même n , size et prob . Je peux construire une fonction pathologique pour laquelle ce n'est même pas vrai :

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

En gros, cette fonction cherche à savoir si les dixièmes de seconde du temps sont impairs ou pairs pour décider de tirer ou non un nombre aléatoire. Exécutez cette fonction et .Random.seed peut ou non changer :

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

Le mieux que vous puissiez (devriez ?) espérer est qu'un ensemble donné de code avec todos les entrées/paramètres identiques (y compris la graine) donnera toujours des résultats identiques. S'attendre à des résultats identiques lorsque seulement le plus (ou seulement certains) des paramètres sont les mêmes n'est pas réaliste, à moins que toutes les fonctions appelées ne donnent ces garanties.