I am working on implementing the Bayesian hierarchical clustering algorithm found here from scratch as a way to practice and learn the algorithm. However, I have hit a snag in calculating the quantity in formula 1 of section 2:

\begin{align} p(D_k|H_1^k) &= \int p(D_k|\theta) p(\theta |\beta) d\theta\\ &= \int \left[ \prod_{x_i \in D_k} p(x_i |\theta) \right] p(\theta | \beta) d\theta \end{align}

In the subsequent paragraph, the authors state,

This is a natural model-based criterion for measuring how well the data fit into one cluster. If we choose models with conjugate priors (e.g. Normal-Inverse-Wishart priors for Normal continuous data or Dirichlet priors for Multinomial discrete data) this integral is tractable

I understand that this is the probability of cluster $k$ being generated from some distribution with parameter(s) $\mathbf{\theta}$ and that, on the first iteration of the algorithm, I should be calculating this quantity for every pair of data points (i.e. vectors) in the set.

However, I am a bit confused regarding the quote above. Does this mean that I should replace the contents of the integral with an appropriate posterior distribution and then perform the integration with respect to $\mathbf{\theta}$?

What would this mean in the Gaussian example mentioned by the authors:

For example, in the case of Gaussians, $p(D_k|H_1^k)$ is a function of the sample mean and covariance of the data in $D_k$

The paper seems to allude to this being a straight-forward calculation, but I have not seen conjugate priors used in this way. Any help would be appreciated!


Hope this isn't too late to help! In short, yes, you have what I believe is the right idea.

I've been messing around with this a little bit myself, and what's being referenced here is the fact that the marginal likelihood (the integral expression you provide) has a quick-to-evaluate closed form that only makes use of the sample statistics of $D_{k}$ when the conjugate prior is employed, no numerical integration required.

Rather than do the mathematical legwork, I'll point you to a resource for the multivariate Gaussian/Normal-Inverse-Wishart case. See Section 9, and specifically, the derivation of the marginal likelihood in 9.5. It can be a little hard to read because the notation is sprinkled throughout the reference, but ultimately, you just need to compute expression 266.

One more tangentially-related thing. Practically speaking, you'll actually want to be dealing with the marginal log-likelihood, as the gamma function risks overflow. I heartily recommend you compute the (pseudo?)log-odds to make BHC merge decisions: $$\log(P(D_{k}|H_{1}^{k})P(H_{1}^{k})) - \log(P(D_{k}|H_{2}^{k})P(H_{2}^{k}))$$ Then, you can just plug the marginal log-likelihood right in.

Forgive any flubs, this is my first StackExchange-style response.

| cite | improve this answer | |
  • 1
    $\begingroup$ Welcome to Stats.SE and thank you for your answer. If you want, you may take the tour. $\endgroup$ – Ertxiem - reinstate Monica Jun 1 '19 at 21:50

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.