0
$\begingroup$

Ex:Show the procedures to simulate an random variable that follows a binomial distribution with parameter $p$, using the $\mathscr{U}(0,1)$(Uniform distribution on the interval (0,1)).

I tried to solve this question by using the following theorem:

Theorem: Let $U\sim\mathscr{U}(0,1) $. Let $X$ be a random variable with distribution $F_X(x)$. Therefore the random variable Y=F^{-1}(U) has a distribution function equal to $F_X$, the distribution of the $X$ variable.

According to this theorem I would need to find a the inverse of the binomial c.d.f, define it as a function in python and generate random numbers.

However I have no idea on how to invert the Binomial distribution.

Questions:

1) Is this the simplest method to simulate a Binomial distribution with the Uniform(0,1)? Are there other methods?

2) How do I compute the inverse of the Binomial distribution?

Thanks in advance!

$\endgroup$
  • 5
    $\begingroup$ Simulating $n$ iid Bernoulli variables is easier. $\endgroup$ – GNUSupporter 8964民主女神 地下教會 Apr 1 '18 at 14:40
  • $\begingroup$ @GNUSupporter Indeed, but this requires generating $n$ uniform random variates. The problem doesn't specify whether this is allowed :) $\endgroup$ – Math1000 Apr 1 '18 at 15:18
  • $\begingroup$ @Math1000 When one can simulate one uniform random variable, one can also simulate $n$ copies. If you insist on generating one single variable and use the inversion theorem, you'll have to compute a discrete sum $\sum_{i = 1}^k p^i (1-p)^{n-i}$ in order to find out $P(X \le k)$. I don't think it's a sensible simulating algo. $\endgroup$ – GNUSupporter 8964民主女神 地下教會 Apr 1 '18 at 15:25
1
$\begingroup$

It may be helpful to see how this procedure for random sampling from $\mathsf{Binom}(n=5, p=.4)$ is implemented in R statistical software. First, some R notation: runif (without extra parameters) is a source of pseudorandom values from standard uniform; dbinom, pbinom and qbinom denote binomial PDF, CDF, and quantile funcions (invese CDF) respectively.

So in R we can generate $m = 100,000$ observations from $\mathsf{Binom}(n=5, p=.4)$ in a vector xas follows.

set.seed(4118)
m = 10^5;  u = runif(m)
x = qbinom(u, 5, .4)     # inverse CDF transformation

Then we can tally the results, make a histogram of them, and plot exact PDF values on the histogram for comparison:

table(x)/m
x
      0       1       2       3       4       5 
0.07790 0.25775 0.34608 0.23105 0.07744 0.00978 

hist(x, prob=T, br=(0:6)-.5, col="skyblue2", main="100,000 Realizations of BINOM(5,.5)")
k = 0:5; pdf = dbinom(k, 5, .4)
points(k, pdf, col="red") 

enter image description here

Within the accuracy of the graph, it is not possible to distinguish the simulated results (heights of histogram bars) from the theoretical ones (small red circles).

Finally, by using the same seed for the pseudoranom generator as above, we can access exactly the same values u as above. Thus, we can see that R implements this (inverse CDF) method to generate $m$ observations from $\mathsf{Binom}(n=5, p=.4)$ by using the function rbinom defined in R. The tallied results are exactly the same below as above.

set.seed(4118)
m = 10^5;  x = rbinom(m, 5, .4)
table(x)/m
x
      0       1       2       3       4       5 
0.07790 0.25775 0.34608 0.23105 0.07744 0.00978 

Note: By contrast, when $p >.5,$ R uses a slight modification of the inverse CDF method, so that the two approaches give slightly different results. (I used $m = 10,000$ so that differences would be more obvious.)

# Inverse CDF method 
set.seed(401);  m = 10^4;  u = runif(m);  x1 = qbinom(u, 5, .7)
table(x1)/m
x1
     0      1      2      3      4      5 
0.0022 0.0296 0.1323 0.3076 0.3651 0.1632 

# Built-in function 'rbinom'
set.seed(401);  x2 = rbinom(m, 5, .7)
table(x2)/m
x2
     0      1      2      3      4      5 
0.0023 0.0278 0.1288 0.3126 0.3599 0.1686 

Addendum 1: Graphs of CDF of $X \sim \mathsf{Binom}(5, .4)$ and its inverse function. The latter shows $F_X^{-1}(u) = 0,$ for $u < 0.6^5=0.07776,$ as in a Comment.

enter image description here

Addendum 2: Generating $X \sim \mathsf{Binom}(5, .4)$ as the sum of five independent Bernoulli random variables with $p=.4.$ [This is the method suggested in the Comment by @GNUSupporter.]

First let $U_1, U_2, \dots, U_5$ be a random sample from $\mathsf{Unif}(0,1).$ Then $B_i = 1,$ if $U_i \le .4$ and $0$ otherwise. This is essentially a trivial application of the quantile method to a variable that takes only values $0$ and $1$. Then $X = \sum_{i=1}^5 B_i \sim \mathsf{Binom}(5, .4).$ We generate four such binomial random variables below (results: 1, 2, 2, 4). Notice that five pseudorandom uniform values are required for each binomial.

set.seed(1234)
u = runif(5); b = (u < .4);  x = sum(b);  x # sum of logical vec b is nr of its TRUEs
[1] 1
u = runif(5); b = (u < .4);  x = sum(b);  x 
[1] 2
u = runif(5); b = (u < .4);  x = sum(b);  x 
[1] 2
u = runif(5); b = (u < .4);  x = sum(b);  x 
[1] 4

Next we use a for loop to iterate this procedure $m = 100,000$ times. Because we start with the same seed as above, the first four iterations repeat the realizations of $\mathsf{Binom}(5, .4)$ shown just above. A tally of all $m$ results shows results similar to those in the initial simulation of this Answer, closely matching the target distribution.

set.seed(1234);  m = 10^5;  x = numeric(m)
for (i in 1:m) {
  u = runif(5); b = (u < .4);  x[i] = sum(b)  }
mean(x);  x[1:4]
[1] 1.99806  # aprx E(X) = 5(.4) = 2
[1] 1 2 2 4  # same first four realizations of X as above
table(x)/m 
x
      0       1       2       3       4       5 
0.07766 0.25971 0.34683 0.22878 0.07675 0.01027 
$\endgroup$
  • $\begingroup$ Thanks for your answer! I think what I am asked to do is to generate or simulate a Binomial distribution from a Uniform distribution. I can easily using the theorem provided in my post to simulate an exponential function by inverting by hand its expression, then generating random numbers ,defining the inverse exponential as function on Python. However for distributions that are difficult to reverse such as this. I do not know how to proceed since I am not allowed to use commands that already relate to the distribution(on this case Binomial). $\endgroup$ – Pedro Gomes Apr 1 '18 at 18:26
  • $\begingroup$ Is there a way to generate it given the requirements? $\endgroup$ – Pedro Gomes Apr 1 '18 at 18:27
  • $\begingroup$ I guess you'll have to use a sequence of if-then statements as implied by @NCh's Answer, which 'inverts the binomial CDF'. Explicitly, for X~BINOM(5, .4), you have $P(X = 0) = P(X \le 0) = .6^5 = 0.07776.$ So take $X = 0$ if $U \le 0.07776,$ And so on. Hope this helps. (You're not working a 'toy' exercise, but learning about something that is actually used. I was not trying to add another answer for your question, but to discuss implementation of the idea in software.) $\endgroup$ – BruceET Apr 1 '18 at 18:46
  • $\begingroup$ Directly to your two questions: (1) Inverse CDF is a simple method. Context may determine what is 'simplest'. Typically software will already have programmed a quantile function, so using it is trivial. Starting from scratch @GNU may be right that it's easier to program $n$ Bernoulli RVs and add. (2) As mentioned in my previous comment, NCh has shown you how to invert the binomial CDF (it seems, without explicitly using those words). // I'm off to brunch, will check back later. $\endgroup$ – BruceET Apr 1 '18 at 18:58
  • $\begingroup$ See two addenda to my Answer. Hope they are helpful. $\endgroup$ – BruceET Apr 2 '18 at 1:26
1
$\begingroup$

For the second question: If $ U\leq (1-p)^n$, set $Y=0$.

If $(1-p)^n< U\leq (1-p)^n+np(1-p)^{n-1}$, set $Y=1$.

If $(1-p)^n+np(1-p)^{n-1}< U\leq(1-p)^n+np(1-p)^{n-1}+\binom{n}{2}p^2(1-p)^{n-2}$, set $Y=2$.

And so on.

$$Y=k \iff \sum_{i=0}^{k-1} \binom{n}{i}p^i(1-p)^{n-i} < U \leq \sum_{i=0}^{k} \binom{n}{i}p^i(1-p)^{n-i}.$$

On your first question: as GNU Supporter noted in comments, this is not the simplest way to generate binomial distribution.

$\endgroup$
  • 1
    $\begingroup$ What do you think about question (1)? Perhaps you won't code this in a real simulation. $\endgroup$ – GNUSupporter 8964民主女神 地下教會 Apr 1 '18 at 15:54
  • $\begingroup$ The answer to this question is obvious and is given in your comment above. Thank you, I'll add this to the answer. $\endgroup$ – NCh Apr 1 '18 at 15:59
  • $\begingroup$ Thanks for reply. I just want your viewpoint due Math1000's comment. $\endgroup$ – GNUSupporter 8964民主女神 地下教會 Apr 1 '18 at 16:01
  • 1
    $\begingroup$ But simulating a Bernoulli r.v. is also "inverting" its CDF. Anyways, thanks for your opinion. $\endgroup$ – GNUSupporter 8964民主女神 地下教會 Apr 1 '18 at 16:11
  • 1
    $\begingroup$ There exist probability spaces on which a binomial-distributed random variable exists but cannot be written as the sum of independent $\{0,1\}$-valued random variables. See here math.stackexchange.com/questions/1920276/… $\endgroup$ – Math1000 Apr 1 '18 at 17:33

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.