In [1]: 
rand(5)  #returns 5 random numbers between 0 and 1
Out[1]:
array([ 0.26846914,  0.17006191,  0.65676757,  0.58864605,  0.4034436 ])
In [2]: 
sum([x<1./6 for x in rand(10)])  #True=1, False=0
Out[2]:
1

Recall "Bernoulli process": probability \(p\) of success in each trial.
Probability of \(m\) successes in N trials: \(P(m)={N\choose m}p^m(1-p)^{N-m}\).
In numpy \({N\choose m}={N!\over m!(N-m)!}\) is given by the function comb(N,m):

In [3]: 
from scipy.misc import comb
def bernoulli(m,N,p): return comb(N,m)* p**m * (1-p)**(N-m)

Now roll a die, say "success" is roll a 6. Then \(p=1/6\). In 60 rolls, what is the probability of \(m\) successes? (where necessarily \(0\le m\le 60\))

In [4]: 
T=10000  # number of "experiments"
N=60     # number of trials per "experiment"
plot(T*bernoulli(arange(N/2),N,1./6),'go--') #prediction
hist([sum([x<1./6 for x in rand(N)]) for t in xrange(T)],bins=range(N/2),align='left')
xlim(-.5,N/2)
xlabel('# of sixes in {} rolls'.format(N))
ylabel('# times occurred in {} "experiments"'.format(T))
legend(['prediction','data']);
In [14]: