rand(5) #returns 5 random numbers between 0 and 1
sum([x<1./6 for x in rand(10)]) #True=1, False=0
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)
:
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\))
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']);