# How do I generate a random number according to the binomial distribution?

 1 I would like to generate a random integer according to the binomial distribution. That is, I would like to generate a Bin(n,p) random value. That's number between 0 and n in which the probability we get the value k is C(n,k)p^k(1-p)^(n-k). Here is an inefficient method (which requires n calls to random() ): def bin_rv(n,p=0.5): """ Generate a binomial random variable with parameters n,p. """ return sum( random() < p for _ in range(n))  Is there a better way to do this (already built into Sage, I hope)? (In Matlab, this can be done via the binord function and in R with the rbinom function.) asked Apr 02 '12 Ed Scheinerman 31 ● 1 ● 2 ● 5

 1 I usually use scipy for this: sage: import scipy.stats sage: scipy.stats.[TAB] Display all 213 possibilities? (y or n) scipy.stats.Tester scipy.stats.dweibull scipy.stats.gmean scipy.stats.linregress scipy.stats.pointbiserialr scipy.stats.t scipy.stats.alpha scipy.stats.entropy scipy.stats.gompertz scipy.stats.loggamma scipy.stats.poisson scipy.stats.test scipy.stats.anderson scipy.stats.erlang scipy.stats.gumbel_l scipy.stats.logistic scipy.stats.powerlaw scipy.stats.threshold scipy.stats.anglit scipy.stats.expon scipy.stats.gumbel_r scipy.stats.loglaplace [.. etc.] sage: B = scipy.stats.binom(10, 0.5) sage: B sage: B.[TAB] B.args B.cdf B.dist B.entropy B.interval B.isf B.kwds B.mean B.median B.moment B.pdf B.pmf B.ppf B.rvs B.sf B.stats B.std B.var sage: B.mean() 5.0 sage: B.rvs(15) array([6, 5, 5, 5, 4, 7, 4, 7, 6, 4, 6, 7, 4, 7, 5])  posted Apr 02 '12 DSM 4882 ● 12 ● 65 ● 105 Thanks, but this does not appear to be faster. In my application, I'm generating lots of Bin(n,0.5) RVs, but with a different n each time. Ed Scheinerman (Apr 02 '12)
 0 Since you like R: sage: L = r.rbinom(100,10,.5) sage: L [1] 5 4 6 4 3 6 9 1 7 5 5 4 6 3 5 6 5 1 5 6 2 4 7 4 7 7 5 3 3 7 6 3 2 2 4 6 4 4 6 4 4 5 6 5 7 3 7 [48] 7 3 5 4 7 7 6 4 3 2 4 5 5 7 2 5 2 3 2 3 4 6 7 8 6 6 3 4 2 6 3 5 3 8 6 8 5 5 2 5 5 5 2 5 7 5 5 [95] 7 3 5 5 6 4 sage: M = L._sage_() sage: type(M), len(M) (, 100)  You do get Python ints, not Sage Integers, though that wouldn't be hard to change with a list comprehension, presumably. posted Apr 02 '12 kcrisman 7427 ● 17 ● 76 ● 166 Or map, which is faster.Eviatar Bach (Apr 03 '12)Unfortunately, the R routine is much slower than the naive one for my application. Ed Scheinerman (Apr 05 '12)

[hide preview]