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$.
Inference about the concentration parameter:
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?