Searching for proof - bayesian inference for exponential distribution According to Wikipedia (https://en.wikipedia.org/wiki/Conjugate_prior) the gamma distribution is a conjugate prior for the exponential distribution (with unknown rate-parameter, $\lambda$, and hyperparameters $\alpha$ and $\beta$). Moreover the posterior predictive is the Lomax (a.k.a. Pareto type II) distribution.
While I have no doubt that these results are correct I have not been able to find any proof leading to the Lomax distribution (the part concerning the gamma distribution is easy to find). I would appreciate if someone would share a reference.
 A: It's easy to show if you have an integer value for $\alpha$ in the prior.
We have that the posterior is $\textrm{Gamma}(\alpha^* = \alpha+n,\beta^* = \beta+\sum_{i=1}^n x_i)$, which means that the posterior density is given by:
\begin{equation}
p[\lambda | \alpha^*, \beta^*] = \frac{(\beta^*)^{\alpha^*}}{\Gamma(\alpha^*)}\lambda^{\alpha^*-1}\exp[-\beta^*\lambda]
\end{equation}
The likelihood is an exponential density:
\begin{equation}
p[x_{new} | \lambda] = \lambda\exp[-\lambda x_{new}]
\end{equation}
Multiplying these two together and moving things around gives us:
\begin{equation}
\frac{(\beta^*)^{\alpha^*}}{\Gamma(\alpha^*)}\lambda^{\alpha^*}\exp[-(\beta^*+x_{neq})\lambda]
\end{equation}
Integrating out $\lambda$, we get:
\begin{equation}
\begin{split}
& \hspace{6mm }\frac{(\beta^*)^{\alpha^*}}{\Gamma(\alpha^*)}\int_{\lambda \in (0,\infty)}\lambda^{\alpha^*}\exp[-(\beta^*+x_{neq})\lambda]\\ & = \frac{(\beta^*)^{\alpha^*}}{\Gamma(\alpha^*)}\int_{\lambda \in (0,\infty)}\frac{(\beta^*+x_{neq})}{(\beta^*+x_{neq})}\lambda^{\alpha^*}\exp[-(\beta^*+x_{neq})\lambda]\\& =  \frac{(\beta^*)^{\alpha^*}}{\Gamma(\alpha^*)}\frac{1}{(\beta^*+x_{neq})}\int_{\lambda \in (0,\infty)}(\beta^*+x_{neq})\lambda^{\alpha^*}\exp[-(\beta^*+x_{neq})\lambda]
\end{split}
\end{equation}
Which is the $\alpha^*$ moment of an $\textrm{Exp}(\beta^*+x_{neq})$ distribution. This is equivalent to $\frac{\alpha^*!}{(\beta^*+x_{neq})^{\alpha^*}}$; giving us that the previous expression is:
\begin{equation}
\begin{split}
\frac{(\beta^*)^{\alpha^*}}{\Gamma(\alpha^*)}\frac{1}{(\beta^*+x_{neq})}\int_{\lambda \in (0,\infty)}(\beta^*+x_{neq})\lambda^{\alpha^*}\exp[-(\beta^*+x_{neq})\lambda] & = \frac{(\beta^*)^{\alpha^*}}{\Gamma(\alpha^*)}\frac{\alpha^*!}{(\beta^*+x_{neq})^{\alpha^*+1}}
\end{split}
\end{equation}
The $\frac{\alpha^*!}{\Gamma(\alpha^*)}$ term reduces to $\alpha^*$, giving us that the expression reduces to:
\begin{equation}
\begin{split}
\frac{(\beta^*)^{\alpha^*}}{\Gamma(\alpha^*)}\frac{\alpha^*!}{(\beta^*+x_{new})^{\alpha^*+1}} & =\frac{\beta^{\alpha^*}\alpha^*}{(\beta^*+x_{new})^{\alpha^*+1}}
\end{split}
\end{equation}
Which matches the density given here with $\lambda = \beta^*$ and $\alpha = \alpha^*$.

