# deriving posterior conditionals for gibbs sampling

I'm new to Bayesian inference and Gibbs sampling in general, and I'm struggling trying to derive the conditional posteriors for a particular data generating process I'm trying to model. The model I am considering is as follows:

$k \sim Beta(a,b) \\ \lambda \sim Gamma(c,f)$

for i = 1 to N

$d_i \sim Bernoulli(k) \\ s_i \sim Poisson(\lambda) \; if d_i = 1 \\ s_i = 0 \; if d_i = 0$

The idea is that for each $s_i$, we flip a coin with some probability, and if the coin is heads we produce a value from the Poisson, and if it is tails, we produce a 0. I hope my notation above isn't too incorrect.

Given observed $s_1 \ldots s_n$, I want to infer the parameters which most likely generated the observations. I think Gibbs sampling can help with this.

The way I've written out the joint posterior is as follows:

$P(k, \lambda, d_{1:n}, s_{1:n}) = \left[k^{a-1}(1-k)^{b-1}\right] \cdot \left[\lambda^{c-1} e^{-\lambda/f}\right] \cdot \prod\limits_{i=1}^{n} \left[ k^{d_i}(1-k)^{1-d_i} \right] \cdot \left[(\lambda^{s_i}e^{-\lambda})^{d_i} (0^{s_i})^{1 - d_i} \right]$

The first square-bracket enclosed part corresponds to the beta prior, the second to the Gamma prior, the third to the Bernoulli toss, and the fourth to the "Poisson or 0" component. I think I've defined this correctly, but I had some trouble figuring out how to write out the "conditional if" part of the process notationally.

My question now is, how can I go about deriving the conditional posteriors from here? I think ultimately I need to derive conditional posteriors for $k$ and $\lambda$, but the $d_i$ terms are confusing me -- should I have a conditional posterior for each $d_i$ as well? It seems that to derive and draw from the posterior conditionals of any variables, I need to know what the values of $d_i$ are for each $i$. Sorry if this vaguely worded, but I would really appreciate any help in figuring out the direction here.

• Unless you want to marginalize out the $d_i$, then yes you need to sample the conditional posterior for the $d_i$ as well (they're latent variables). It might help to draw out your model as a graph (something like this en.wikipedia.org/wiki/Plate_notation) so it's clear where the dependencies lie. – aleshing Dec 28 '17 at 18:28

Alright so I decided to actually write this all out for practice on my own part. Let's go through the conditional posterior for each parameter/ latent variable (let me know if anything is unclear). I first start off by writing the model down in a way that makes it easier for me to see all the dependencies: $$s_i|d_i=0 \sim \delta_0(\cdot), \ s_i|d_i=1,\lambda \sim \textsf{Pois}(\lambda),\ d_i|k\sim\textsf{Bern}(k),$$ $$\lambda\sim\textsf{Gamma}(c,f),\ k\sim\textsf{Beta}(a,b),$$ where $\delta_x(\cdot)$ is a point mass at $x$. Thus we see that the posterior conditional of $k$ will only depend on the $d_i$, the posterior conditional of the $d_i$ will depend on the corresponding $s_i$, $k$, and $\lambda$, and the posterior conditional of $\lambda$ will only depend on the $d_i$ and the $s_i$.

We can write the posterior conditional of $k$ as: $$p(k|\{d_i\})\propto p(\{d_i\}|k)p(k)=\textsf{Beta}(k|a+\#\{i|d_i=1\}, b+\#\{i|d_i=0\})$$ due to beta Bernoulli conjugacy (here $\textsf{Beta}(\cdot|\alpha,\beta)$ is the beta pdf with parameters $\alpha$ and $\beta$ and $\#$ just denotes the cardinality of a set).

We can write the posterior conditional of each $d_i$ as: $$p(d_i=1|s_i,k,\lambda)\propto p(s_i|d_i=1,\lambda)p(d_i=1|k)=\textsf{Pois}(s_i|\lambda)k,$$ $$p(d_i=0|s_i,k,\lambda)\propto p(s_i|d_i=0)p(d_i=0|k)=\delta_{0}(s_i)(1-k).$$ We thus see that after normalizing that $$p(d_i|s_i,k,\lambda)=\textsf{Bern}\left(d_i\bigg|\frac{\textsf{Pois}(s_i|\lambda)k}{\textsf{Pois}(s_i|\lambda)k+\delta_{0}(s_i)(1-k)}\right).$$

We can write the posterior conditional of $\lambda$ as: $$p(\lambda|\{d_i\}, \{s_i\})\propto p(\{s_i\}|\{d_i\},\lambda)p(\lambda)=\prod_{i=1}^n(\textsf{Pois}(s_i|\lambda))^{\delta_1(d_i)}\textsf{Gamma}(\lambda|c,f)=\prod_{i|d_i=1}\textsf{Pois}(s_i|\lambda)\textsf{Gamma}(\lambda|c,f)=\textsf{Gamma}(\lambda|\hat{n}\hat{s}+c,\hat{n}+f)$$ where $\hat{n}=\#\{i|d_i=1\}$ and $\hat{s}=\frac{1}{\hat{n}}\sum_{i|d_i=1}s_i$, due to the gamma Poisson conjugacy.

A good reference (imo) for starting out with Gibbs samplers is http://www.cs.columbia.edu/~blei/fogm/2015F/notes/mixtures-and-gibbs.pdf. The plate notation for the Gaussian mixture model there will be nearly the same as your model since your model is also a mixture (instead with a point mass and a Poisson instead of with different Gaussians).

When we're done running our Gibbs sampler according to the posterior conditionals above, we'll have samples from the joint posterior $p(\lambda, k, \{d_i\}|\{s_i\})$. We can throw out the $\{d_i\}$, which would be equivalent to marginalizing them out and having samples from $p(\lambda, k|\{s_i\})$, which I'll call $\{\lambda_j\}$, $\{k_j\}$ for $j=1,\cdots,T$, where $T$ is the number of samples. We then want to compute the posterior predictive distribution $$p(s^*|\{s_i\})=\int p(s^*|k,\lambda)p(\lambda, k|\{s_i\})dkd\lambda$$ for a new data point $s^*$. We can write $$p(s^*|k,\lambda)=\sum_{i=0}^1p(s^*|d^*=i,\lambda)p(d^*=i|k),$$ so that the integral we want to calculate is $$p(s^*|\{s_i\})=\int \left[\sum_{i=0}^1p(s^*|d^*=i,\lambda)p(d^*=i|k)\right]p(\lambda, k|\{s_i\})dkd\lambda.$$ Using the sample, we can approximate this integral using Monte Carlo integration as $$p(s^*|\{s_i\})\approx\frac{1}{T}\sum_{j=1}^T \left[\sum_{i=0}^1p(s^*|d^*=i,\lambda_j)p(d^*=i|k_j)\right],$$ where $p(s^*|d^*=i,\lambda_j)$ will be a Poisson pdf with parameter $\lambda_j$ evaluated at $s^*$ if $i=1$, or the point mass $\delta_0(\cdot)$ evaluated at $s^*$ if $i=0$, and $p(d^*=i|k_j)$ is the Bernoulli pmf with parameter $k_j$ evaluated at $i$.

• Thanks for your detailed answer! I'm in the process of grokking through it, but on a high-level it makes sense. I had one more question: The main reason I wanted to build this model was to evaluate the likelihood of observing possible future $s_i$. For those, I wouldn't know $d_i$, but could compute likelihood using posterior means for k and lambda, right? If that's my use-case, would I want to marginalize out $d_i$ for parameter estimation? If I don't marginalize out $d_i$, would I essentially have to compute an expected likelihood over $p(d_i | s_i, k, \lambda)$? – voidobscura Dec 28 '17 at 19:55
• What you're talking about is the posterior predictive distribution, which are talked about on pages 3 and 4 of cs.columbia.edu/~blei/fogm/2015F/notes/mixtures-and-gibbs.pdf. I'll update my answer with how to compute this since it's too long for the comments. – aleshing Dec 28 '17 at 20:15
• Thanks so much. I'm perusing the document as well, but I appreciate your walkthrough as it's refreshingly clear given my comparative lack of stats background... – voidobscura Dec 28 '17 at 20:21
• No problem! If you have the time, it'd probably be worthwhile, if you haven't already, to go through an introductory probability book (like Introduction to Probability Theory by Hoel et al or A First Course in Probability by Ross) along with an introductory bayesian statistics book (like Statistical Rethinking by McElreath or A First Course in Bayesian Statistical Methods by Hoff). It'll help build intuition on how to work with these types of probabilistic models. – aleshing Dec 28 '17 at 20:38