deriving posterior conditionals for gibbs sampling

Question

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.

Answer 1

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$.