Cacher is the code snippet organizer for pro developers

# Plot the model evidence of the Bernoulli distribution with conjugate prior.

import numpy as np
import matplotlib.pyplot as plt
import scipy

def evidence(N1, N0, alpha, beta):
    B = scipy.special.beta
    return c * B(alpha+N1, beta+N0) / B(alpha, beta)

# N1: The num of heads.
# N0: The num of tails.
N1, N0 = 10, 5

# Hypter parameters for conjugate prior
alphas = np.arange(1, 101)
betas = np.arange(1, 101)

aa, bb = np.meshgrid(alphas, betas)
models = np.c_[aa.ravel(), bb.ravel()]

evidences = np.empty(models.shape[0])
for i, m in enumerate(models):
    evidences[i] = evidence(N1, N0, m[0], m[1])

fig = plt.figure(figsize=(8,4))
cnt = plt.contourf(aa, bb, evidences.reshape(aa.shape))
plt.colorbar(cnt)
plt.title('(N1, N0) = (%d, %d)' % (N1, N0))
plt.xlabel('$\alpha$')
plt.ylabel('$\beta$')
plt.savefig('model_evidence_of_.png')
plt.show()