# Gibbs sampler for Dirichlet Process concentration parameter

I am trying to implement a Gibbs sampler for Hierarchical Dirichlet process, but I cannot seem to correctly estimate the concentration parameters. I therefore started testing just this part of a sampler, and found an issue that I don't believe was mentioned here before. Then again, I do suspect I'm just making some silly mistake - all inputs highly appreciated.

Dirichlet process:

I first simulate a Dirichlet process, particularly the number of components $$k$$ that is generated in the process, which depends on the number of observations $$n$$ and the concentration parameter $$\alpha$$.

A classic explanation of this is the Chinese restaurant process, i.e. each observation can be though of as a customer entering a restaurant. The customer joins a table with probability proportional to the number of customers already sitting at this table, or sits at an empty (new) table with probability proportional to the concentration parameter $$\alpha$$. So for the $$i$$-th observation the probability for being assigned to a particular component can be expressed as: $$p(t_i = t) = \begin{cases} \frac{n_t}{\sum_t (n_t) + \alpha}, & \mbox{for } t \leq t_{\text{max}} \\ \frac{\alpha}{\sum_t (n_t) + \alpha}, & \mbox{for } t = t_{\text{max}}+1 \end{cases}$$ where $$n_t$$ is the number of observations already assigned to component $$t$$ and $$t_{\text{max}}$$ is the current number of components. Important for later:

$$\mathbb{E}(k) \approx \alpha \ln\Big(\frac{\alpha+n}{\alpha}\Big)$$

or more generally - the number of components $$k$$ increases with both $$\alpha$$ and $$n$$.

Now, having simulated the Dirichlet process I have the number of observations $$n$$, the number of components $$k$$ and know the true concentration parameter $$\alpha$$. I would like to test if my sampler is able to identify the true $$\alpha$$. I am using a method proposed by Escobar and West (1995) - eq. (13-14), particularly the form used by Teh et al. (2005) - appendix A, and Heinrich (2008) - eq. (18-19). While the form differs a bit and to be honest I don't see how they are equivalent, I tested both distributions and they are indeed the same.

Assuming a gamma prior, $$\alpha \sim \text{Gamma}(a, b)$$, the sampling procedure first samples two auxiliary variables based on the old value of $$\alpha$$, and then uses them to sample the new value for $$\alpha$$.

$$u \sim \text{Bern}\Big(\frac{n}{n+\alpha}\Big)$$ $$v \sim \text{Beta}(\alpha + 1, n)$$ $$\alpha \sim \text{Gamma}(a + K - 1 + u, b - \ln v)$$

Problem:

This sampler provides absolutely wrong estimates for $$\alpha$$. Example: I simulate a DP with $$n=1000$$, $$\alpha=1$$, resulting with $$k=7$$ - a number very close to the expected value. I then run the sampler, with prior $$\text{Gamma}(1,1)$$ (doesn't really have much impact with such number of observations) for 1000 iterations, starting from the true value for $$\alpha$$. The sampler then quickly achieves stability, with average $$\hat{\alpha}\approx 35$$. Utterly wrong.

I feel like I'm missing some crucial part of the method. Especially since it seems to me like the formula for the sampler makes no sense - here's my rationale:

For a fixed value of $$k$$ my estimate of $$\alpha$$ should decrease in $$n$$ (if I have a huge sample with very few components, that tells me my $$\alpha$$ is tiny). Yet if $$n \to \infty$$, then $$u \sim \text{Bern}(\frac{n}{n+\alpha})$$ is almost surely $$1$$, while $$v \sim \text{Beta}(\alpha + 1, n)$$ is almost surely $$0$$ since $$\mathbb{E}(v) = \frac{\alpha + 1}{\alpha + 1 + n}$$.

Then the expected value of my new $$\alpha$$ is: $$\mathbb{E}(\alpha) = (a+k-1+u) (b - \ln v) \to \infty$$ because of "$$\ln 0$$". So for huge number of observations I will have huge estimate of $$\alpha$$. Running the Gibbs sampler with fixed $$k$$ and initial $$\alpha$$, while increasing $$n$$ confirms that suspicion. Now, this doesn't really make any sense, does it? Or am I making a mistake somewhere?

So just as I thought, it was a silly mistake. While the authors repeatedly refer to the second Gamma distribution parameter as scale, it's actually the rate parameter (so to write in accordance with e.g. NumPy standards, $$\alpha \sim \text{Gamma}(a+k-1+u,\frac{1}{b-\ln v})$$. I'm a bit ashamed of how long it took me to realise that, but I'll leave the question on the very off chance someone has the same problem.