016_04.py

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

tau = 5
beta = 1/tau
alpha = 5
t = 15
lam = beta*t

d = sps.gamma(a=alpha, scale=1/beta)
x = np.linspace(0, d.ppf(.999), 1000)
y = d.pdf(x)
p = d.sf(t)

D = sps.poisson(mu=lam)
X = np.arange(0, D.ppf(.999), 1)
Y = D.pmf(X)
P = D.cdf(alpha-1)

fig, ax = plt.subplots(1, 2, figsize=(12, 5), facecolor="w")
ax[0].plot(x, y)
ax[0].fill_between(
    x[x>=t],
    0, y[x>=t],
    color="r", alpha=.25
)
ax[0].plot([0, 0], [0, y[0]], "r")
ax[0].plot([t, t], [0, y[x<=t][-1]], "r")
ax[0].axhline(0, color="k", alpha=.25)
ax[0].set_xlabel("Wating time $t$")
ax[0].set_ylabel("Probability Density")
ax[0].text(16, .005, f"$P(T > 15)={p*100:.3f}\;\%$")
ax[0].set_title(fr"$T \sim \mathrm{{Gamma}}(\alpha={alpha}, \beta={beta})$")
ax[0].set_xticks(np.arange(0, max(x)+5, 5))

ax[1].stem(X, Y)
ax[1].bar(X[X<=alpha-1], Y[X<=alpha-1], color="r", alpha=.25)
ax[1].text(5, .12, f"$P(A < {alpha})={P*100:.3f}\;\%$")
ax[1].set_title(fr"$A \sim \mathrm{{Poisson}}(\lambda={lam:.0f}$)")
ax[1].set_xlabel(r"Completed processes $\alpha$")
ax[1].set_ylabel("Probability Mass")
ax[1].set_xticks(X)

plt.show()