0
$\begingroup$

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))

enter image description here

$\endgroup$
1
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Claude COULOMBE Mar 9 '19 at 19:09
  • $\begingroup$ 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. $\endgroup$ – Claude COULOMBE Mar 9 '19 at 19:57
  • 1
    $\begingroup$ @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! ;) $\endgroup$ – Tomáš Tunys Mar 10 '19 at 11:58
  • $\begingroup$ @ 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. $\endgroup$ – Claude COULOMBE Mar 10 '19 at 18:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.