number_of_samples.py

import numpy as np 
from matplotlib import pyplot as plt

def sample_gamma(alpha):
    total = 0
    scale = 1.

    if (alpha < 1.0):
        scale *= (1 - np.random.uniform() ** (1.0 / alpha))
        total += 1
        alpha += 1.0

    d = alpha - 1.0 / 3.0
    c = 1.0 / np.sqrt(9.0 * d)

    while True:
        x, y = -1, -1

        while y <= 0:
            x = np.random.randn()
            y = 1.0 + c * x
            total += 1

        v = y * y * y
        u = 1 - np.random.random()
        total += 1
        xx = x * x
        if (u < 1.0 - 0.0331 * xx * xx):
            return total
        if (np.log(u) < 0.5 * xx + d * (1.0 - v + np.log(v))):
            return total

def main():
    alpha = np.linspace(0.01, 10, 100)
    counts = np.zeros(len(alpha))
    maxes = np.zeros(len(alpha))
    for i,a in enumerate(alpha):
        for j in range(2000):
            samples = sample_gamma(a)
            counts[i] += samples
            maxes[i] = max(maxes[i], samples)
    counts /= 2000
    plt.plot(alpha, counts, label='average')
    plt.plot(alpha, maxes, label='max samples (n=2000)')
    plt.legend()
    plt.tight_layout()
    plt.xlabel(r'$\alpha$')
    plt.ylabel('random samples')
    plt.show()
 
main()

And here's the plot of the maximum number of samples used across all 2000 trials: