1
$\begingroup$

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?

$\endgroup$
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.