4
$\begingroup$

I am trying to simulate a coin toss using R. For that I used both sample function and rbinom function.

But I am getting different results.

> set.seed(1)
> rbinom(10,1,0.5)
 [1] 0 0 1 1 0 1 1 1 1 0

> set.seed(1)
> sample(c(0,1), 10, replace = TRUE)
 [1] 0 1 0 0 1 0 0 0 1 1

What may be the reason for that?

I thought for the same seed value, both sample function and rbinom should give the same result.

$\endgroup$
5
  • 3
    $\begingroup$ You have used the same seed, but the two procedures rbinom and sample use pseudorandom numbers in different ways. $\endgroup$ – BruceET Jan 30 at 6:12
  • $\begingroup$ Related: stackoverflow.com/questions/64965398/… $\endgroup$ – Peter O. Jan 30 at 13:01
  • 2
    $\begingroup$ You should mention your version of R because the default algorithm for sample has changed over time (IIRC most recently at 3.6.0). $\endgroup$ – Chris Haug Jan 30 at 15:56
  • 2
    $\begingroup$ A specific seed only reliably gives the same results when you call the exact same random functions in the exact same order (basically when you rerun the same code multiple times). Just because 2 functions output the same distribution of values for some sets of inputs doesn't mean what's happening behind the scenes is the same. This could even apply to functions that use the exact same high-level algorithm, since they could either use different lower-level random functions or the order they get their random values in might just not be the same. $\endgroup$ – Bernhard Barker Jan 30 at 19:03
  • $\begingroup$ Would this answer your question? R Difference between rbinom and sample $\endgroup$ – Xi'an Feb 9 at 7:41
8
$\begingroup$

Somewhat related example: One way to generate 10 tosses of a coin with probability $0.4$ of heads is to use rbinom:

set.seed(123); rbinom(10, 1, .4)
[1] 0 1 0 1 1 0 0 1 0 0

Another way is to use the binomial inverse CDF (quantile) function) qbinom to transform uniform random numbers from runif get the desired Bernoulli distribution.

set.seed(123); qbinom(runif(10), 1, .4)
[1] 0 1 0 1 1 0 0 1 0 0

This suggests that R uses qbinom with runif to get rbinom--in this instance.

However, for success probabilities greater than $0.5,$ it seems that R uses a variant of this method, and results starting with the same seed differ.

set.seed(123);  rbinom(10, 1, .6)
 [1] 1 0 1 0 0 1 1 0 1 1
set.seed(123);  qbinom(runif(10), 1, .6)
 [1] 0 1 1 1 1 0 1 1 1 1

If you use the same seed and you access pseudorandom numbers (as with runif) in exactly the same order to do exactly the same procedure, you will get the same result.

Another example: Means of 100 samples of size 10 from the population $\mathsf{Norm}(\mu = 50, \sigma = 7).$

set.seed(1234)
a = replicate(100, mean(rnorm(10, 50, 7)))
summary(a)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  44.64   48.41   49.69   49.81   51.37   55.41 

set.seed(1234)
x = rnorm(1000, 50, 7)
MAT = matrix(x, byrow = T, nrow=100)
a = rowMeans(MAT)
summary(a)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  44.64   48.41   49.69   49.81   51.37   55.41 

However, if we omit byrow=T in making the matrix, R will use its default, which is to fill the matrix by columns, and we will get a different summary (except, of course, for the mean of means).

set.seed(1234)
x = rnorm(1000, 50, 7)
MAT = matrix(x, nrow=100)
a = rowMeans(MAT)
summary(a)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  42.80   48.47   49.80   49.81   51.20   55.27 
$\endgroup$
1
  • 3
    $\begingroup$ For p > 0.5, the random numbers are drawn according to size - rbinom(n, size, 1-p) in case anyone wonders. $\endgroup$ – AlexR Jan 30 at 23:06
7
$\begingroup$

An historical note about sample is that it got recently modified for being biased in some extreme situations (as commented by Chris Haug). In earlier versions of R such as 3.4.4, still running on my ownlaptop, the outcome of the above would be the same as a cdf inversion (and as the complement of the standard Uniform draw):

> set.seed(1)
> rbinom(10,1,0.5)
 [1] 0 0 1 1 0 1 1 1 1 0

> set.seed(1)
> sample(c(0,1), 10, replace = TRUE)
 [1] 0 0 1 1 0 1 1 1 1 0

> set.seed(1)
> qbinom(runif(10),1,0.5) #inverse cdf
 [1] 0 0 1 1 0 1 1 1 1 0

> set.seed(1)
> 1*(runif(10)>0.5) #complement!
 [1] 0 0 1 1 0 1 1 1 1 0

When checking the C code behind the R function rbinom, the adopted approach relies on a single Uniform when $np<30$:

/* inverse cdf logic for mean less than 30 */
repeat {
 ix = 0;
 f = qn;
 u = unif_rand();
 repeat {
 if (u < f)
     goto finis;
 if (ix > 110)
     break;
 u -= f;
 ix++;
 f *= (g / ix - r);

and a much more involved resolution otherwise ($np\ge 30$), resolution including an accept-reject step,

/*------- np = n*p >= 30 : ----- */
repeat {
  u = unif_rand() * p4;
  v = unif_rand();

hence a random number of Uniforms. On the other hand, sample (in a version of 2017!) uses a cascade of C functions:

if (replace) {
    int i, nc = 0;
    for (i = 0; i < n; i++) if(n * p[i] > 0.1) nc++;
    if (nc > 200)
    walker_ProbSampleReplace(n, p, INTEGER(x), k, INTEGER(y));
    else
    ProbSampleReplace(n, p, INTEGER(x), k, INTEGER(y));
} else
    ProbSampleNoReplace(n, p, INTEGER(x), k, INTEGER(y));

For instance, ProbSampleReplace is based on a single Uniform call:

/* compute the sample */
for (i = 0; i < nans; i++) {
  rU = unif_rand();
  for (j = 0; j < nm1; j++) {
      if (rU <= p[j])

and the other ones as well, which is not to say that they return the same outcome than rbinom!

$\endgroup$
2
$\begingroup$

I tried getting to the source code for both, but couldn't find it.

I did, however, find the references for the building of the two algorithms. They do not use the same references, so it is reasonable that they do not generate them in a similar manner. Indeed, sample() function briefly says it uses an easier way to handle random numbers, presumably through different generation.

$\endgroup$

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