Chinese Restaurant Process implementation in Python
Question: I have wrote a code in Python for CRP problem. The problem itself can be found
here: <http://cog.brown.edu/~mj/classes/cg168/slides/ChineseRestaurants.pdf>

And to give a short description of it: Suppose we want to assign people
entering to a restaurants to potentially infinite number of tables. If $z_i$
represents the random variable assigned for the $i$'th person entering the
restaurant the following should hold:

With probability $p(z_i=a|z_1,...,z_{i-1})=\frac{n_a}{i-1+\alpha} for $n_a>0$,
$i$'th person will sit in table $a$ and with probability
$p(z_i=a|z_1,...,z_{i-1})=\frac{\alpha}{i-1+\alpha} $i$'th person will sit
around a new table.

I am not quite sure if my code is correct cause I am surprised how small the
final number of tables are.

I would be happy if somebody could say if the implementation is correct and if
so are there any possible improvements.

    
    
    import numpy as np
    def CRP(alpha,N):
        """Chinese Restaurant Process with alpha as concentration parameter and N 
        the number of sample"""
        #Array which will save for each i, the number of people people sitting
        #until table i
        summed=np.ones(1) #first person assigned to the first table
        for i in range(1,N):
            #A loop that assigns the people to tables
    
            #randind represent the random number from the interval [1,i-1+alpha]
            randind=(float(i)+alpha)*np.random.uniform(low=0.0, high=1.0, size=1)
            #update is the index for the table that the person should be placed which
            #if greater than the total number, will be placed in a new table
            update=np.searchsorted(summed,randind,side='left')
            if randind>i:
                summed=np.append(summed,i+1)
            else:
                zerovec=np.zeros(update)
                onevec=np.ones(summed.size-update)
                summed+=np.append(zerovec,onevec)
        #This part converts summed array to tables array which indicates the number
        #of persons assigned to that table
        tables=np.zeros(summed.size)
        tables[0]=summed[0]
        for i in range(1,summed.size):
            tables[i]=summed[i]-summed[i-1]
        return tables
    a=CRP(0.9999,1000)
    print a
    


Answer: Suggestion. Forget about the code you have written. Construct declarative
tests of the code. By taking that approach, you start with examples for which
you know the correct answer. That would have answered Brainiac's question, for
example.

Then write your program. You will likely find that if you start approaching
problems this way, you may create sub-problems first, for which you can also
write tests. Until they all pass, there is no need to rush on to the full
problem.

