title: Chinese Restaurant Process implementation in Python
body: <p>I have wrote a code in Python for CRP problem. The problem itself can be found here:
  <a href="http://cog.brown.edu/~mj/classes/cg168/slides/ChineseRestaurants.pdf" rel="nofollow">http://cog.brown.edu/~mj/classes/cg168/slides/ChineseRestaurants.pdf</a></p>
  
  <p>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:</p>
  
  <p>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.</p>
  
  <p>I am not quite sure if my code is correct cause I am surprised how small the final number of tables are.</p>
  
  <p>I would be happy if somebody could say if the implementation is correct and if so are there any possible improvements.</p>
  
  <pre><code>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&gt;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
  </code></pre>