Rewriting numpy.random.binomial using a For loop

Question

Currently I have the following code

def approx_binomial(n, p, size=None): gaussian = np.random.normal(n*p, n*p*(1-p), size=size) # Add the continuity correction to sample at the midpoint of each integral bin. gaussian += 0.5 if size is not None: binomial = gaussian.astype(np.int64) else: # scalar binomial = int(gaussian) return binomial

However, it is not very accurate as it uses the random function. Is there any other way in which I can rewrite the function using a for loop?

Another query I have is how to display a graph of the probability mass function against the number of successes?

def approx_binomial(n, p, size=None): gaussian = np.random.normal(n*p, n*p*(1-p), size=size) # Add the continuity correction to sample at the midpoint of each integral bin. gaussian += 0.5 if size is not None: binomial = gaussian.astype(np.int64) else: # scalar binomial = int(gaussian) return binomial plt.plot(n,p) plt.show()

Thanks!

1. What do you mean by the function is not very accurate? The function seems supposed to be returning random samples. What has randomness to do with accuracy, here? Please be more specific in your question about what you expect vs. what you get. 2. The second question seems rather tangential to the first. Please consider to separate them in two posts for clarity. — MB-F
– MB-F, Commented Feb 2, 2016 at 7:52
@kazemakase What I mean is that I need a more accurate recoding of the numpy.random.binomial. The current code I have has values that differ a lot from what you would get when you use the aforementioned function. — user3393463
– user3393463, Commented Feb 2, 2016 at 7:56
What's wrong with numpy's binomial? You could simply use that. — MB-F
– MB-F, Commented Feb 2, 2016 at 7:58
I see. That is important information, you know :) This place is not intended to solve your homework, but I will try to point you in the right direction. Gimme a moment. — MB-F
– MB-F, Commented Feb 2, 2016 at 8:04

Jaime · Accepted Answer · 2016-02-02 11:13:37Z

The obvious way of getting a random sample from a binomial distribution using a for loop is:

import random def binomial_random_sample(n, p): ret = 0 for j in range(n): if random.random() < p: ret += 1 return ret

Or if you prefer a more terse approach:

def binomial_random_sample(n, p): return sum(random.random() < p for j in range(n))

MB-F · Accepted Answer · 2016-02-02 08:51:04Z

The way to solve your task can be found on Wikipedia, but it is a very small section:

One way to generate random samples from a binomial distribution is to use an inversion algorithm. To do so, one must calculate the probability that P(X=k) for all values k from 0 through n. (These probabilities should sum to a value close to one, in order to encompass the entire sample space.) Then by using a linear congruential generator to generate samples uniform between 0 and 1, one can transform the calculated samples U[0,1] into discrete numbers by using the probabilities calculated in step one.

The above is a bit verbose. In a nutshell, you will have to do the following to generate one sample from the binomial distribution: Use the probability mass function p=pmf(n, k) to compute the binomial probability for each possible k, and randomly select one of the ks based on their probability

Compute the cumulative density function of the binomial distribution cdf(n, k) for every k = [0..n].
Generate a uniformly distributed random number r between 0 and 1
Find the value of k for which cdf(n, k) <= r.

Collectives™ on Stack Overflow

Rewriting numpy.random.binomial using a For loop

2 Answers 2

1 Comment

Comments

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

1 Comment

Comments

Related