Подтвердить что ты не робот

Ошибочное поведение семян с rbinom (prob = 0,5) в R

Я нашел то, что я считаю неустойчивым поведением (но для которого я надеюсь, что есть простое объяснение) в R использовании семян в сочетании с rbinom(), когда используется prob=0.5. Общая идея. Для меня, если я установил семя, запустите rbinom() один раз (т.е. Выполните один случайный процесс), несмотря на то, что установлено значение prob, случайный семя должно меняться на один шаг. Затем, если я снова установил семя в одно и то же значение и запустил еще один случайный процесс (например, rbinom() снова, но, возможно, с другим значением prob), семя снова должно измениться до того же значения, что и для предыдущего одиночного случайного процесса.

Я нашел, что R делает именно это, пока я использую rbinom() с любым prob!=0.5. Вот пример:

Сравните вектор-семестр, .Random.seed, для двух вероятностей, отличных от 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

Сравнить вектор-семпл, .Random.seed, для 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
...

Я нашел это для всех сравнений prob=0.5 против всех других вероятностей в множестве {0,1, 0,2,..., 0,9}. Аналогично, если я сравню любые значения prob с {0,1, 0,2,..., 0,9}, кроме 0,5, вектор .Random.seed всегда является поэтапным. Эти факты справедливы и для нечетного или даже size в пределах rbinom().

Чтобы сделать его еще более странным (я извиняюсь, что это немного запутанно - это имеет отношение к тому, как моя функция написана), когда я использую вероятности, сохраненные как элементы в векторе, у меня такая же проблема, если 0.5 является первым элементом, но не второй. Вот пример для этого случая:

Первый случай: 0.5 - первая вероятность, на которую ссылается вектор

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

Второй случай: 0.5 - вторая вероятность, указанная в векторе

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

Опять же, я считаю, что, несмотря на значения, используемые для prob и size, этот шаблон сохраняется. Может ли кто-нибудь объяснить эту тайну мне? Это вызывает проблему, потому что результаты, которые должны быть одинаковыми, различаются, потому что семя по какой-либо причине используется/вычисляется по-разному, когда prob=0.5, но ни в каком другом экземпляре.

4b9b3361

ОТВЕТЫ

Ответ 1

Итак, давайте ответом на наши комментарии. Спасибо Ben Bolker за то, что вы поместили нас на правильный путь со ссылкой на код: https://svn.r-project.org/R/trunk/src/nmath/rbinom.c и предложение отследить, где вызывается unif_rand().

Быстрое сканирование, и кажется, что код разбит на два раздела, ограниченные комментариями:

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

и

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

Внутри каждого из них количество вызовов unif_rand не совпадает (два против одного).

Итак, для данного size (n) ваше случайное семя может оказаться в другом состоянии в зависимости от значения prob (p): size * prob >= 30 или нет.

С учетом этого все результаты, которые вы получили с вашими примерами, теперь должны иметь смысл:

# 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

Ответ 2

Я собираюсь принять противоположную позицию по этому вопросу и утверждать, что ожидания не подходят и не поддерживаются документацией. В документации нет никаких претензий в отношении того, какие побочные эффекты (в частности, на .Random.seed) можно ожидать, вызвав rbinom или как эти побочные эффекты могут или не могут быть одинаковыми в различных случаях.

rbinom имеет три параметра: n, size и prob. Ваше ожидание состоит в том, что для случайного набора семян перед вызовом rbinom, .Random.seed будет таким же после вызова rbinom для заданных n и любых значений size и prob (или, может быть, любые конечные значения size и prob). Вы, конечно, понимаете, что для разных значений n это будет отличаться. rbinom не гарантирует этого или подразумевает это.

Не зная внутренних функций, это невозможно узнать; как показал другой ответ, алгоритм отличается от продукта size и prob. И внутренности могут измениться, чтобы эти конкретные детали могли измениться.

По крайней мере, в этом случае результат .Random.seed будет таким же после каждого вызова rbinom, который имеет те же n, size и prob. Я могу построить патологическую функцию, для которой это не так:

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

В принципе, эта функция смотрит, является ли десятая часть второго времени нечетной или даже решить, нарисовать ли случайное число или нет. Запустите эту функцию и .Random.seed может или не может измениться:

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

Надеемся, что лучше всего вы можете (следует?) надеяться, что данный набор кодов со всеми входами/параметрами одинаковый (включая семя) всегда будет давать одинаковые результаты. Ожидание идентичных результатов, когда только большинство (или только некоторые) параметров одинаковы, нереалистично, если все функции, называемые этими гарантиями, не обеспечивают эти гарантии.