Calculate the posterior of binomial likelihood

conjugate_prior.py

import numpy as np
import scipy.stats as stats

success_prob = 0.3
data = np.random.binomial(n=1, p=success_prob, size=1000) # sucess is 1, failure is 0.

# Domain θ
theta_range = np.linspace(0, 1, 1000)

# Prior P(θ)
a = 2
b = 8
theta_range_e = theta_range + 0.0001 
prior = stats.beta.cdf(x = theta_range_e, a=a, b=b) - stats.beta.cdf(x = theta_range, a=a, b=b) 

# The sampling dist. aka Likelihood P(X|θ)
likelihood = stats.binom.pmf(k = np.sum(data), n = len(data), p = theta_range)

# Posterior
posterior = likelihood * prior
normalized_posterior = posterior / np.sum(posterior)