Bayesian hurdle model and predictions

Question

I fitted a Bayesian hurdle model assuming a Poisson distribution:

$ P(Y_i=y_i) = \begin{cases} 1-\pi_{i} & \mbox{if }y_i=0\\[0.5em] \pi_{i} \frac{\lambda_{i}^{y_i} e^{-\lambda_i}}{y_i!(1-e^{-\lambda_i})} &\mbox{if }y_i > 0 \end{cases} $

Because the log-likelihood is separable with respect to $\pi$ and $\lambda$, I fitted two separate regressions (i.e., logistic and zero-truncated Poisson regression) using a set of predictors. In other words:

$ \lambda = \log(\beta_{1}^T X) \\ \pi = \mbox{logit}(\beta_{2}^T X) $

My issue arises in trying to predict new values. In a Bayesian context, I have posterior distributions for the coefficients from the logistic and zero-truncated Poisson regressions. For a new predicted value, I will have a posterior distribution for $\pi$ and $\lambda$. To get the $E(\pi \cdot \lambda)$ is easy, I can just take the product of the posterior means. But what about getting the 2.5 and 97.5% percentiles for the credible interval around for the new predicted value?

Answer 1

Let's call the posterior predictive distribution $\tilde{Y}$.

$$\begin{align}\tilde{Y}_{\pi} &\sim\text{Bernoulli}(\pi)\\ \tilde{Y}_{\lambda} &\sim\text{Poisson}(\lambda)\\ \tilde{Y} &\sim\tilde{Y}_{\pi}*\tilde{Y}_{\lambda}\end{align}$$

You would need to truncate $\tilde{Y}_{\lambda}$ to be non-zero. Since you would then have a full predictive distribution, it would be easy to calculate intervals, expectations, etc.

Computation wise, this could be implemented in a number of ways. One possibility (in pseudo-code):

y_new = 0
y_pi = bernoulli_rng(pi)
if(y_pi==1)
  while(y_new==0)
    y_new = poisson_rng(lambda)
return y_new