import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import poisson

np.random.seed(seed=1)

left = 0
right = 100
step = 0.1
alpha = 1
beta = 1
times = 20

def bayes_inf(x, step, alpha, beta, mu, times):
    data = poisson.rvs(mu=mu, size=times)
    integral = np.sum(x ** (sum(data) + alpha) * np.exp((-np.shape(data)[0] - beta) * x)) * step
    y = x ** (sum(data) + alpha) * np.exp((-np.shape(data)[0] - beta) * x) / integral
    return y, data

Sampling mu from gamma-like distribution

x = np.arange(left, right, step)
probs = x ** alpha * np.exp(-beta * x) / np.sum(x ** alpha * np.exp(-beta * x) * step)
plt.rcParams['figure.figsize'] = 20, 20
for i in range(1, 10):
    plt.subplot(3, 3, i)
    mu = np.random.choice(x, p=probs/np.sum(probs))
    y, data = bayes_inf(x, step, alpha=alpha, beta=beta, mu=mu, times=times)
    plt.plot(x, probs, color='green', label='Aprior mu distribution')
    plt.axvline(mu, color='red', label="Sampled mu")
    plt.plot(x, y, label="Aposterior distribution of mu")    
    plt.xlim(0, 10)
    plt.ylim(0, 1.75)
    plt.legend()