In [16]: from numpy.random import binomial

In [17]: nHeads = binomial(500_001, 0.5, size=int(1e6))

In [18]: sum(nHeads == 250_000) / len(nHeads)
Out[18]: 0.001171

In [19]: 1/800
Out[19]: 0.00125