gem_dist.py

#!/usr/bin/env python
import sys
import numpy as np
import numpy.random as random
import scipy.stats as stats
import matplotlib.pyplot as plt

if __name__ == '__main__':
    if len(sys.argv) < 2:
        print "usage: $ ./gem_dist.py alpha"
        sys.exit()

    alpha = float(sys.argv[1])
    beta_dist = stats.beta(1, alpha)
    xs = np.linspace(0, 1, num = 100)

    plt.figure(figsize = (10, 8))
    plt.subplot(3, 3, 1)
    plt.plot(xs, beta_dist.pdf(xs))
    plt.title("Beta PDF")

    for i in range(2, 10):
        us = [random.beta(1, alpha) for _ in range(9)]
        betas, rest = reduce(lambda r, x: (r[0] + [r[1] * x], r[1] * (1 - x)), us, ([], 1.0))
        betas += [rest]

        plt.subplot(3, 3, i)
        plt.ylim(0, 1)
        plt.xlim(0, 10)
        plt.bar(range(10), betas)

    plt.show()