<s>[INST] Mixed Poisson distribution with Lindley mixing I want to calculate the unconditional distribution of $x$
$$p(x|p) 
= \int_0^\infty p(x,\theta|p) d\theta = \int_0^\infty p(x|\theta)p(\theta|p) d\theta$$
where $x|\theta \sim Poisson(\lambda)$ and $\lambda$ is a random variable that follows the Lindley distribution $$p(\theta|p) = \frac{p^2}{p+1}(\theta+1)e^{(-\theta p)}$$
The proven result is $$p(x|p) = \frac{p^{2}(p+2+x)}{(p+1)^{(x+3)}}$$
Doing so I have :
\begin{align*}
p(x|\lambda) 
&= \int_0^\infty p(x,\theta|\lambda) d\theta\\
&= \int_0^\infty p(x|\theta)p(\theta|\lambda) d\theta\\
&= \int_0^\infty \frac{e^{-\theta}\theta^x}{x!} \frac{p^2}{p+1}(\theta+1)e^{(-\theta p)}d\theta\\
&= \frac{p^2}{(p+1)x!}\int_0^\infty e^{-\theta(1+p)}\theta^{x}(\theta+1)  d\theta\\
\end{align*}
So far so good. But I split the integral (am I allowed to do that?) into :
$$
\frac{p^2}{(p+1)x!}\int_0^\infty e^{-\theta(1+p)}\theta^{x}(\theta+1)  d\theta
 = \frac{p^2}{(p+1)x!}\int_0^\infty e^{-\theta(1+p)}\theta^{x+1}  d\theta +\int_0^\infty e^{-\theta(1+p)}\theta^{x} d\theta $$

Edit:
Given that :
\begin{align*}
\int_0^\infty \theta^x e^{-\theta (p+1)} d\theta
&=\frac{\Gamma(x+1)}{(p + 1)^{(x+1)}} \int_0^\infty \frac{(p + 1)^{(x+1)}}{\Gamma(x+1)}\theta^x e^{-\theta (p+1)} d\theta\\
&=\frac{\Gamma(x+1)}{(p + 1)^{(x+1)}} \int_0^\infty \text{Gamma}(\theta; x+1, p+1) d\theta\\
&=\frac{\Gamma(x+1)}{(p + 1)^{(x+1)}}
\end{align*}
then becomes : $$\frac{p^2}{(p+1)\Gamma(x+1)} \left(\frac{\Gamma(x+1)}{(p+1)^{x+1}} + \frac{\Gamma(x)}{(p+1)^x}  \right)$$
Am I right ?
 [/INST] Assuming your work up to the last line (before edit) is correct we have with a little algebraic manipulation
\begin{align}
f_X(x|p)
&=\frac{p^2}{(p+1)x!}\int_0^\infty e^{-\theta(1+p)}\theta^{x}(\theta+1)\,\mathrm d\theta\\
&=\frac{p^2}{(p+1)^{x+2}}\int_0^\infty (1+\theta)\frac{(p+1)^{x+1}}{\Gamma(x+1)}\theta^{(x+1)-1}e^{-\theta(p+1)}\,\mathrm d\theta,
\end{align}
which can be written as the expected value
$$
f_X(x|p)=\frac{p^2}{(p+1)^{x+2}}\mathsf E(1+X),\quad X\sim\operatorname{Gamma}(x+1,p+1).
$$
By linearity of the expected value
$$
\mathsf E(1+X)=1+\mathsf EX=1+\frac{x+1}{p+1};
$$
hence,
$$
f_X(x|p)=\frac{p^2}{(p+1)^{x+2}}\left(1+\frac{x+1}{p+1}\right).
$$

Note:
If you need to show the steps in evaluating $\mathsf EX$ then write
$$
\mathsf EX=\int_0^\infty \theta\frac{(p+1)^{x+1}}{\Gamma(x+1)}\theta^{(x+1)-1}e^{-\theta(p+1)}\,\mathrm d\theta.
$$
Substituting $u=(p+1)\theta$ and making use of the integral definition of the gamma function will get you the final result.

Edit:
As requested by the OP:
\begin{align}
f_X(x|p)
&=\frac{p^2}{(p+1)x!}\int_0^\infty e^{-\theta(1+p)}\theta^{x}(\theta+1)\,\mathrm d\theta\\
&=\frac{p^2}{(p+1)}\int_0^\infty(1+\theta) \frac{1}{\Gamma(x+1)}e^{-\theta(1+p)}\theta^{x}\,\mathrm d\theta\\
&=\frac{p^2}{(p+1)}\int_0^\infty(1+\theta) \frac{1}{\Gamma(x+1)}\theta^{(x+1)-1}e^{-(p+1)\theta}\,\mathrm d\theta\\
&=\frac{p^2(p+1)^{-(x+1)}}{(p+1)}\int_0^\infty(1+\theta) \frac{(p+1)^{x+1}}{\Gamma(x+1)}\theta^{(x+1)-1}e^{-(p+1)\theta}\,\mathrm d\theta\\
&=\frac{p^2}{(p+1)^{x+2}}\int_0^\infty(1+\theta) \frac{(p+1)^{x+1}}{\Gamma(x+1)}\theta^{(x+1)-1}e^{-(p+1)\theta}\,\mathrm d\theta\\
&=\frac{p^2}{(p+1)^{x+2}}\int_0^\infty(1+\theta) \operatorname{Gamma}(\theta|x+1,p+1)\,\mathrm d\theta\\
\end{align}
</s>