import scipy.stats as st
rands = np.random.geometric(1./50,10000)
obs = bincount(rands)
rv = st.geom(1./50)
exp = rv.pmf(range(len(obs))) * obs.sum()
plot(range(1,len(obs)),obs[1:],'r',range(1,len(exp)),exp[1:],'b')
legend(["observed","expected"]);
print "chisq stat: %.2f, p-value: %f" % st.chisquare(obs[1:],exp[1:])
chisq stat: 500.61, p-value: 0.249091
st.chisqprob(500.61, len(obs)-2)
0.2490955165306051