4
$\begingroup$

I am sure I am just misunderstanding this but rbinom and sample seem to be generating disparate results and I am not sure why. I am trying to create a categorical factor with a balanced third in each levels (so .33 get 0, .33 get 1, .33 get 2).

If I attempt to sample evenly 1/3 from rbinom(), I get uneven distributions:

table(rbinom(100,2,c(.33,.33,.33)))
 0  1  2 
47 42 11 
 0  1  2 
49 41 10
 0  1  2 
46 40 14 
 0  1  2 
45 43 12 
 0  1  2 
46 38 16 

The third level is not getting many successes and I am not sure why.

If I do the same thing for the sample() function, I get more satisfactory results.

table(sample(c(1,2,3),100,replace = T,prob = c(.33,.33,.33)))
 1  2  3 
32 32 36 
 1  2  3 
39 33 28 
 1  2  3 
30 41 29 
 1  2  3 
33 34 33 
1  2  3 
36 25 39 
1  2  3 
39 31 30 

Here, each level (1, 2, and 3) is getting around 33% of successes which is what I wanted.

How do I get rbinom() to do the same thing, distributing successes evenly across 1, 2, and 3?

$\endgroup$
2
  • 3
    $\begingroup$ You want the the multinomial, not the binomial distribution. Try rmultinom(5, 100, c(.33,.33,.33)). $\endgroup$
    – dipetkov
    Commented Jun 21, 2022 at 6:06
  • 2
    $\begingroup$ for rbinom, the probability of getting a 2 is from a size 2 with success of 0.33 is 0.1089. Hence 11% of the time you will get a 2. and that is exactly what you have. dbinom is not same as discrete uniform. You should use rmultinom as suggested or use sample as you are doing so far $\endgroup$
    – Onyambu
    Commented Jun 21, 2022 at 6:52

2 Answers 2

8
$\begingroup$

As noticed in the comments, your code does not do what you want it to do. First, many of R's functions vectorized the passed parameter values, in case of your code rbinom(100,2,c(.33,.33,.33)) is the same as rbinom(100,2,.33), or rbinom(100,2,c(.33,.33)), or rbinom(100,2,rep(.33, 100)) because all the values provided to prob are the same, if they differed, it'd cycle through them until drawing the intended 100 samples.

The code itself produces draws from a binomial distribution with sample size=2, so you are independently tossing a biased coin with prob=0.33 two times. Your second snippet draws each of the values in c(1,2,3) with an equal probability of 0.33. So the first code would produce the result 2 with probability 0.33 * 0.33 = 0.1089 (two independent coin throws), while the second the result 3 with a probability of 0.33 (you roll a three-sided dice). In the first case, the probability distribution is

$$ \Pr(X=k) = {2 \choose k} \, 0.33^k \,0.66^{2-k} $$

while in the second case, it's just $\Pr(X=k) = 0.33$. They are not the same.

This is what you observe if you run the code.

> n <- 1000000
> table(rbinom(n, 2, .33)) / n

       0        1        2 
0.448398 0.443002 0.108600 
> table(rbinom(n, 2, c(.33,.33,.33))) / n

       0        1        2 
0.449160 0.441886 0.108954 
> table(sample(c(1,2,3), n, replace = T, prob = c(.33,.33,.33))) / n

       1        2        3 
0.333505 0.332922 0.333573

You can achieve the same result as when using rbinom with sample if you use this unnecessarily slow code:

> table(replicate(n, sum(sample(c(0, 1), size=2, replace=T, prob=c(.66,.33))))) / n

       0        1        2 
0.444294 0.444688 0.111018 
$\endgroup$
0
5
$\begingroup$

When calling

rbinom(100,2,c(a,b,c))

the outcome is equivalent to

c(rbinom(2,a),rbinom(2,b),rbinom(2,c),rbinom(2,a),rbinom(2,b),...)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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