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.

  • 1
    $\begingroup$ 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. $\endgroup$
    – 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$.

  • $\begingroup$ 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)$? $\endgroup$ Dec 28 '17 at 19:55
  • $\begingroup$ 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. $\endgroup$
    – aleshing
    Dec 28 '17 at 20:15
  • $\begingroup$ 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... $\endgroup$ Dec 28 '17 at 20:21
  • $\begingroup$ 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. $\endgroup$
    – aleshing
    Dec 28 '17 at 20:38

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.