# Sampling from conditional posterior - continuous and discrete terms

I'm working through Hierarchically Supervised LDA by Perotte et all (2011). The conditional posterior I'm supposed to sample values $z_i$ from, however, is zero almost everywhere.

To see why, lets have a look at the conditional posterior I use in my Gibbs sampler: $$\label{cond_post} p(z_{i}=k| \mathbf{z}_{-i}, \mathbf{a}, \mathbf{w}, \eta)\propto \underbrace{\frac{n_{-i, (k)}^{w_i}+\beta}{n_{-i, (k)} + V\beta} \cdot \frac{n_{-i, (d)}^k + \alpha}{n_{-i, (d)} + K\alpha}}_{\text{ "Term1" }=p(z_{i}=k| \mathbf{z}_{-i}, \mathbf{w})} \cdot \underbrace{\prod_{l=1}^L \exp \left\{ -\frac{1}{2}(\bar{\mathbf{z}}_d^T\eta_l - a_{l,d})^2 \right\}}_{ \text{"Term2" }= p(\mathbf{a} | \mathbf{z}, \eta)}$$

where $k \in \{1, 2, ..., K\}$. This expression may seem overwhelming at first, but I'll lead you through the relevant aspects of it further down below. For a more detailed background of HSLDA (such as the generative model and the graphical model) and the derivation of the conditional posterior, please refer to this post.

Sampling values for $z_i$ from this conditional posterior has been causing me trouble. To see why, lets focus on the two terms:

• Term 1: Is a discrete density with positive probability for all $K$ values. By enumeration, this term can easily be transformed to a proper density. I don't think that any more knowledge about this term is required at this point.
• Term 2: Is a Gaussian kernel, where the parameter for the mean is a function of $z_i$ (and $\mathbf{z}_i$ and $\eta$). It would be easy enough to sample values for $a_{l}$ from this term only, because it is essentially a joint normal distribution with $L$ independent Gaussian variables.

Note that the conditional posterior above is equivalent to the joint conditional posterior of $\mathbf{a}$ and $z_i$: \begin{align*} p(z_{i}=k| \mathbf{z}_{-i}, \mathbf{w}) \cdot p(\mathbf{a} | \mathbf{z}, \eta) &= p(\mathbf{a} | z_{i}, \mathbf{z}_{-i}, \mathbf{w}, \eta)\cdot p(z_{i}=k | \mathbf{z}_{-i}, \mathbf{w}, \eta)\\ &= p(\mathbf{a}, z_{i}=k | \mathbf{z}_{-i}, \mathbf{w}, \eta)\end{align*}

The problem: Looking at the expression for conditional posterior of $z_i$, Term2 will always produce a value close to zero, because it is the product of $L$ very small values. In comparison, Term1 always produces value somewhere around $1/K$, that are much larger than Term2. This poses a problem, because normalising (enumeration) gives me an Error: can't divide by zero, which means that $\sum_{k=1}^K p(z_i = k | \mathbf{z}_{-i}, \mathbf{a}, \mathbf{w}, \eta) = 0$. This must be caused by multiplying a value around $1/K$ ( = $p(z_i|\mathbf{z}_{-i}, \mathbf{w})$) with $L$ values close to zero ( = $\prod_{l=1}^L p(\mathbf{a}|\mathbf{z}, \eta)$).

What I tried so far:

• I thought about magnifying Term2 by some constant to avoid getting only values close to zero, but then I'd artificially increase the influence of Term2 compared to Term1.
• Multiplying both terms by some constant wouldn't solve the problem either, because the influence of Term2 would be nil.
• There's one more thing I've been thinking about: $a_{l}$ is a running variable in a probit regression. Its value is only of secondary interest, because the actual important thing is whether $a_{l}$ is larger or smaller than zero. So I thought that perhaps it'd be more sensible to evaluate $\prod_{l=1}^L p(a_{l} > 0 | \mathbf{z}, \eta)$ instead of $\prod_{l=1}^L p(a_{l} | \mathbf{z}, \eta)$ (= Term2). The major problem I have with this is that the author did not do this either. Before I start doubting the author, I reckon it is more likely I made a mistake.

My question: How can I use the conditional posterior for Gibbs sampling values for $z_i$? I am looking for a theoretically sound approach that allows for a practical implementation of this Gibbs sampler.

• How big is $L$? – jbowman Nov 22 '17 at 1:25
• L is usually about 10 – KenHBS Nov 22 '17 at 6:04
• Your problem seems somehow wrong to me. Even if each term in the product had a value of 0.001, the product would only be around $10^{-30}$, which is easily handled by a computer. Since your Gaussian kernel doesn't have a scale parameter, I'm assuming it's 1 for all $L$ components, and the density at the extreme (for a standard Gaussian) value of 3 is only about 0.004, so you really shouldn't be seeing this happen. Are you sure the code isn't substituting "0" for the Gaussian probability calculation somehow? – jbowman Nov 22 '17 at 16:13
• Yes, the scale is one for all components. $L$ is a variable, that's at least 4 and at most 22. I will double-check the small values. My suspicion is that the values for $\bar{\mathbf{z}}_d^T\eta_l$ vary a lot and very often contain quite extreme values, especially in the beginning of the sampling. The problem is that due to these extreme values, there is no way for the sampling process to stabilize to more probable values. – KenHBS Nov 22 '17 at 16:40
• It takes a while before I can check the final results, but first inspections seem very promising. It turns out that the initial values for $a_{l,d}$ had been programmed to be generated from a distribution different than $\mathcal{N}(\bar{\mathbf{z}}_d^T\eta_l)$, so therefore I indeed got values that were 40 standard deviations away from the mean, causing the Gaussian Kernel to result in 0. Thanks a million for pointing this out! – KenHBS Nov 23 '17 at 10:43