# Coding a simple Stick-Breaking Process in Python

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 =[1]+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))


Correct me if I am wrong, but shouldn't betas sum to 1.0 to make a valid distribution?

I have changed your code snippet just to do that (mind that the last sample from beta distribution is not used):

def Stick_Breaking(num_weights, alpha):
betas = np.random.beta(1, alpha, size=num_weights)
betas[1:] *= np.cumprod(1 - betas[:-1])
return betas


For example, one sample from Stick_Breaking(10, 1.0) give array([0.8698612 , 0.11382917, 0.01005555, 0.00410595, 0.00137072]) and np.sum(Stick_Breaking(10, 1.0)) gives values close to 1.0.

EDIT: I stand corrected that the outcomes is not a distribution for any finite num_weights, but the code should be correct.

• Nice! Indeed, it should sum up to one. I was suspecting something was wrong, but you put your finger on it. I would correct it mention you! Many thanks. – Claude COULOMBE Mar 9 '19 at 19:09
• That said with big alpha parameter (like alpha=10) the spreading effect is obvious but the sum up is under one. It's more complicated than I figure out at the beginning. – Claude COULOMBE Mar 9 '19 at 19:57
• @ClaudeCOULOMBE True, I was wrong assuming the result should be a distribution. It actually is when num_weights tends to $inf$, but it has nothing to do with it. I should remind myself some of this stuff :). If you are interested into DP I highly recommend this talk by Michael Jordan on YT: youtube.com/watch?v=yfLoxwjCGNY. Have fun! ;) – Tomáš Tunys Mar 10 '19 at 11:58
• @ Tomáš Tunys Thanks Tomáš, great talk by a great researcher! Thanks for the sharing. By the way, you help me a lot. I'm fascinated by all those China restaurant, Indian Buffet, Urn-Draw and Stick Breaking. I was wrong, you were imprecise, we're correct now. It's how the science and open source software progress and share. Keep contributing. – Claude COULOMBE Mar 10 '19 at 18:31