1
$\begingroup$

I have stumbled upon R behavior involving set.seed(), sample() and rbinom() which I don't understand.

Background: I am solving problems from Richard McElreath's "Statistical Rethinking". Medium problems 1-3 in Ch.3 basically ask to do the following:

Suppose you throw a coin 15 times and get 8 heads. Given this and a uniform prior for $p_{success}$, approximate posterior distribution using grid approximation, sample from it, and conduct posterior predictive check using these samples (i.e. generate new data from the model and see whether the observed data is likely).

Here is the code I use to compute samples and posterior predictive check:

# calculate posterior
p_grid<- seq( from=0 , to=1 , length.out=1000 )
prior <- rep( 1 , 1000 )
likelihood <- dbinom( 8 , size=15 , prob=p_grid )
posterior <- likelihood * prior
posterior <- posterior / sum(posterior)

# sample from it
set.seed(500)
samples <- sample( p_grid , prob=posterior , size=1e4 , replace=TRUE )

# posterior predictive check
set.seed(500)
ppc <- rbinom(1e4, size = 15, prob = samples)

# plot the results of posterior predictive check
rethinking::simplehist(ppc)

Question: if I use two different values in the two calls to set.seed(), or don't call it at all, I get the expected picture: 8 successes in 15 trials are predicted to be the most frequent outcome.

enter image description here

However, if I use the same value in the two calls to set.seed(), as in the code above, the output is completely unexpected:

enter image description here

Does anybody have any idea what is going on here? I know that normally one would just call set.seed() once in the beginning of the session, so I am guessing the pattern I see is somehow related to using the same random numbers in sample() and rbinom(); but does anybody know in what exact way this leads to this weird bimodal histogram?

$\endgroup$
1
$\begingroup$

Both these functions (sample and rbinom) will use uniform pseudo-random numbers to obtain their random values. By resetting the seed you'll start the sequence of uniforms in the same place for both the call to sample and the call to rbinom. As a result, you will have dependent sequences.

That's most likely what the issue is in your results.

For example, try comparing:

set.seed(500)
x=runif(10000)
y=runif(10000)
hist(x+y,n=50)

with

set.seed(500)
x=runif(10000)
set.seed(500)
y=runif(10000)
hist(x+y,n=50)

Now check out:

summary(x-y)
$\endgroup$
  • $\begingroup$ Thanks! I thought something along those lines... Is there any particular mathematical reason for why the resulting "weird" distribution is bimodal? Also, I've experimented further and I think the issue I describe relates to this. If in the rbinom call I set size = 60, I still get a weirdly looking distribution, but for size = 61 and higher the distribution is much more reasonable and looks like in the first plot. So it looks like it is related to the differences between np>=30 and np<30 branches in the underlying C code. $\endgroup$ – Adahn Oct 17 '17 at 3:57
  • 1
    $\begingroup$ A detailed answer would require examining exactly how sample and rbinom use the numbers -- studying the code in enough detail to figure out what happened to each and every value generated from your seed. All manner of weird things may happen (e.g. 1. for some parts of the domain the uniforms may end up flipped or sliced up in various ways and so have complex relationships with whatever happens to the other variate; e.g. 2. if rbinom uses accept-reject anywhere, the sequence may be shifted after each rejection, so they may start of perfectly dependent but end up shifted ... ctd $\endgroup$ – Glen_b Oct 17 '17 at 4:12
  • 1
    $\begingroup$ ctd... along the sequence of uniforms by such a rejection, and hence at some point may start to act like they're independent, even though initially they were perfectly lined up). This is way more of a task that I care to undertake -- it could well end up in 'takes a small book to answer" territory. $\endgroup$ – Glen_b Oct 17 '17 at 4:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.