Solution (source code)

= Solution

Use $A=a+\sum_ix_i$ and $B=b+n$ as above. The future observation is independent of the old observations conditional on $\theta$, so the <posterior predictive probability> is obtained by integrating its <Poisson distribution> against the <posterior distribution> of $\theta$:
$$
\begin{aligned}
\Pr(X_{n+1}=k\mid x)
&=\int_0^\infty\frac{e^{-\theta}\theta^k}{k!}\frac{B^A}{\Gamma(A)}\theta^{A-1}e^{-B\theta}\,d\theta\\
&=\frac{B^A}{\Gamma(A)k!}\int_0^\infty\theta^{A+k-1}e^{-(B+1)\theta}\,d\theta\\
&=\boxed{\frac{\Gamma(A+k)}{\Gamma(A)k!}\left(\frac B{B+1}\right)^A
\left(\frac1{B+1}\right)^k,\quad k=0,1,\ldots.}
\end{aligned}
$$
The last integral is the normalizing <gamma function> integral, equal to $\Gamma(A+k)/(B+1)^{A+k}$. This is a <negative binomial distribution> in the count-of-failures convention with shape $A$ and success probability $B/(B+1)$; the gamma-function expression remains valid for noninteger $A>0$. As checks, its <mean> is $A/B$, and the <law of total variance> gives <variance> $A/B+A/B^2$. The additional term is the contribution of posterior uncertainty in $\theta$.