I've just red the great 2012 blog post of Edwin Chen about Dirichlet Process with companion code in R and Ruby. Then I'm trying to translate the Stick-Breaking Process from R to Python.
I've got this small piece of code. Is-it good?
# A direct translation of the R code import numpy as np def Stick_Breaking(num_weights,alpha): betas = np.random.beta(1,alpha, size=num_weights) remaining_stick_lengths =+list(np.cumprod(1-betas))[0:num_weights-1] return remaining_stick_lengths * betas # A nice correction suggested by Tomáš Tunys def Stick_Breaking(num_weights,alpha): betas = np.random.beta(1,alpha, size=num_weights) betas[1:] *= np.cumprod(1 - betas[:-1]) return betas
In order to show small distribution histograms:
import matplotlib.pyplot as plt for _ in range(5): num_weights = 10 alpha = 1 weights = Stick_Breaking(num_weights,alpha) plt.axis([0, num_weights+1, 0, 1]) plt.bar(range(1,num_weights+1),weights) plt.show() print(sum(weights))