title: How to extract unsupervised clusters from a Dirichlet Process in PyMC3?
body: <p>I just finished the <a href="https://rads.stackoverflow.com/amzn/click/com/1785883801" rel="noreferrer" rel="nofollow noreferrer">Bayesian Analysis in Python</a> book by <a href="https://github.com/aloctavodia" rel="noreferrer">Osvaldo Martin</a> (great book to understand bayesian concepts and some fancy numpy indexing).</p>
  
  <p>I really want to extend my understanding to bayesian mixture models for unsupervised clustering of samples. All of my google searches have led me to <a href="https://pymc-devs.github.io/pymc3/notebooks/dp_mix.html" rel="noreferrer">Austin Rochford's tutorial</a> which is really informative.  I understand what is happening but <strong>I am unclear in how this can be adapted to clustering</strong> (especially using multiple attributes for the cluster assignments but that is a different topic).  </p>
  
  <p>I understand how to assign the priors for the <code>Dirichlet distribution</code> but I can't figure out how to get the clusters in <code>PyMC3</code>.  It looks like the majority of the <code>mus</code> converge to the centroids (i.e. the means of the distributions I sampled from) but they are still separate <code>components</code>.  I thought about making a cutoff for the <code>weights</code> (<code>w</code> in the model) but that doesn't seem to work the way I imagined since multiple <code>components</code> have slightly different mean parameters <code>mus</code> that are converging. </p>
  
  <p><strong>How can I extract the clusters (centroids) from this <code>PyMC3</code> model?</strong> I gave it a maximum of <code>15</code> components that I want to converge to <code>3</code>. The <code>mus</code> seem to be at the right location but the weights are messed up b/c they are being distributed between the other clusters so I can't use a weight threshold (unless I merge them but I don't think that's the way it is normally done). </p>
  
  <pre><code>import pymc3 as pm
  import numpy as np
  import matplotlib.pyplot as plt
  import multiprocessing
  import seaborn as sns
  import pandas as pd
  import theano.tensor as tt
  %matplotlib inline
  
  # Clip at 15 components
  K = 15
  
  # Create mixture population
  centroids = [0, 10, 50]
  weights = [(2/5),(2/5),(1/5)]
  
  mix_3 = np.concatenate([np.random.normal(loc=centroids[0], size=int(150*weights[0])), # 60 samples
                          np.random.normal(loc=centroids[1], size=int(150*weights[1])), # 60 samples
                          np.random.normal(loc=centroids[2], size=int(150*weights[2]))])# 30 samples
  n = mix_3.size
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/sg1vZ.png" rel="noreferrer"><img src="https://i.stack.imgur.com/sg1vZ.png" alt="enter image description here"></a></p>
  
  <pre><code># Create and fit model
  with pm.Model() as Mod_dir:
      alpha = pm.Gamma('alpha', 1., 1.)
  
      beta = pm.Beta('beta', 1., alpha, shape=K)
  
      w = pm.Deterministic('w', beta * tt.concatenate([[1], tt.extra_ops.cumprod(1 - beta)[:-1]]))
  
      component = pm.Categorical('component', w, shape=n)
  
      tau = pm.Gamma("tau", 1.0, 1.0, shape=K)
  
      mu = pm.Normal('mu', 0, tau=tau, shape=K)
  
      obs = pm.Normal('obs',
                      mu[component], 
                      tau=tau[component],
                      observed=mix_3)
  
      step1 = pm.Metropolis(vars=[alpha, beta, w, tau, mu, obs])
  #     step2 = pm.CategoricalGibbsMetropolis(vars=[component])
      step2 = pm.ElemwiseCategorical([component], np.arange(K)) # Much, much faster than the above
  
      tr = pm.sample(1e4, [step1, step2], njobs=multiprocessing.cpu_count())
  
  #burn-in = 1000, thin by grabbing every 5th idx
  pm.traceplot(tr[1e3::5])
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/1h1Ok.png" rel="noreferrer"><img src="https://i.stack.imgur.com/1h1Ok.png" alt="enter image description here"></a></p>
  
  <p><em>Similar questions below</em></p>
  
  <p><a href="https://stats.stackexchange.com/questions/120209/pymc3-dirichlet-distribution">https://stats.stackexchange.com/questions/120209/pymc3-dirichlet-distribution</a> for regression and not clustering</p>
  
  <p><a href="https://stats.stackexchange.com/questions/108251/image-clustering-and-dirichlet-process">https://stats.stackexchange.com/questions/108251/image-clustering-and-dirichlet-process</a> theory on the DP process</p>
  
  <p><a href="https://stats.stackexchange.com/questions/116311/draw-a-multinomial-distribution-from-a-dirichlet-distribution">https://stats.stackexchange.com/questions/116311/draw-a-multinomial-distribution-from-a-dirichlet-distribution</a> explains DP</p>
  
  <p><a href="https://stackoverflow.com/questions/25685104/dirichlet-process-in-pymc-3">Dirichlet process in PyMC 3</a> directs me to Austin Rochford's tutorial above</p>