``````
from pylab import *
import numpy
from decimal import *
from random import *

getcontext().prec = 200

def RandList(n, t):
return [randrange(1, n) for i in range(0, t)]

def DuplicatesP(t):
return len(unique(t)) < len(t)

def TryBirthday(n, t, tr=1000):
j = 0
for i in range(0, tr):
u = RandList(n, t)
if DuplicatesP(u):
j += 1.0
return j / (tr * 1.0)

def GenPlot(n, tr=5000):
x = list(range(0, n + 2))
y = [TryBirthday(n, i, tr) for i in x]
return [x, y]

def ActualEstimate(i, n):
p = 1
for j in range(0, i):
p *= 1.0 - j / float(n)
return 1 - p

def ExpEstimate(i, n):
l = -i * (i + 1) / Decimal(2 * n)
# SEE Heuer for this constant
const = sqrt(-2 * Decimal(log(1 - 1 / 2.0)))
v = 1 - const * (Decimal(l).exp())
v = float(v)
if v > 0:
return v
return 0

def PrintProbs(n, tr=5000):
for i in range(0, n + 2):
print(
(
"{0:6d} -> value={1:04f} Est={2:04f} Exp={3:04f}".format(
i,
TryBirthday(n, i, tr),
ActualEstimate(i, n),
ExpEstimate(i, n),
)
)
)

def PlotProbs(n, tr=5000):
[x, y] = GenPlot(n, tr)
fig = figure()
ax1.plot(x, y)
ax1.axhline(y=0.5, color="r")
ax1.axvline(x=sqrt(n), color="g")
ax1.axvline(x=1.77 * sqrt(n), color="g")
l = list(range(0, n + 2))
ma = [ExpEstimate(i, n) for i in l]
ax1.plot(l, ma, color="y")
mp = [ActualEstimate(i, n) for i in l]
ax1.plot(l, mp, color="m")
return fig

```
```