foo.py

import numpy as np

def redbrain(n=100, b=0.3, sz=100):
    return [(x, np.random.binomial(n - x, b)) for x in np.random.binomial(n, 0.5, sz)]

def foo(n, sz):
    d = {}
    for i in range(1, 20):
        p = i * 0.05
        xs = redbrain(n=n, b=p, sz=sz)
        for an, qn in xs:
            dd = d[an+qn] if an+qn in d else dict()
            zs = dd[p] if p in dd else []
            zs.append((an, qn))
            dd[p] = zs
            d[an+qn] = dd
    return d