maxent_joint.py

bnds = tuple([(0, None)]*(num_sides2*num_sides1))

# remember number of probabilities
lookup_map = {}
for ni, piqj in enumerate(itertools.product(sides1, sides2)):
    lookup_map[ni] = piqj

# common guesses
guess_common = np.ones(len(lookup_map))/len(lookup_map)

# common constrain on sum to 1
ccons0 = ({'type': 'eq', 'fun': lambda p: np.sum(p) - 1.})

# constrains on mean
def cons20_joint(params):
    s = 0
    for ni, (pij) in enumerate(params):
        s += pij*lookup_map[ni][0] # conditions for p
    return s

def cons21_joint(params):
    s = 0
    for ni, (pij) in enumerate(params):
        s += pij*lookup_map[ni][1] # conditions for q
    return s

ccons20 = ({'type': 'eq', 'fun': lambda x: cons20_joint(x)-mean1})
ccons21 = ({'type': 'eq', 'fun': lambda x: cons21_joint(x)-mean2})

sol_common = minimize(sum_form,
                      guess_common,
                      bounds=bnds,
                      constraints=[ccons0,
                                   ccons20,
                                   ccons21],
                      options={'maxiter':10001})