Integration of Poisson-Exponential Product I have an unnormalised p.d.f I would like to normalise by taking the integral, conditioned on a given value, $k$,
My distribution is the product of a Poisson distribution and an Exponential distribution;
\begin{equation}
f_{X}(x,k) = \text{Pois}(k, \lambda(x))\text{Exp}(\lambda_{0})
\end{equation}
Where
\begin{equation}
\lambda(x) = (1 - e^{-\lambda_{1} x})
\end{equation}
I would therefore like to determine the total probability as a function of $k$
\begin{equation}
P(k) = \int^{\infty}_{0} \frac{\lambda(x)^{k}e^{-\lambda(x)}}{k!} \lambda_{0}e^{-\lambda_{0}x}.dx
\end{equation}
Attempt at a Solution
Let us first attempt the simplest case, $k=0$,
\begin{equation}
P(0) = \lambda_{0} \int^{\infty}_{0} e^{-\lambda(x)} e^{-\lambda_{0}x}.dx
\end{equation}
Substituting in $\lambda(x)$:
\begin{equation}
P(0) = \lambda_{0} \int^{\infty}_{0} e^{-(1 - e^{-\lambda_{1}x})} e^{-\lambda_{0}x}.dx
\end{equation}
This integral looks nasty, so I try to substitute the entire equation and integrate using $\lambda(x)$ as the integration variable, because I know $\lambda(x) \in (0, 1)$
\begin{equation}
x = -\frac{1}{\lambda_{1}}\log(1 - \lambda(x))
\end{equation}
\begin{equation}
dx= \frac{{d\lambda(x)}}{\lambda_{1}(1 - \lambda(x))}
\end{equation}
We can substitute $x$ into  to get:
\begin{equation}
e^{-\lambda_{0}x} = e^{-\lambda_{0}(-\frac{1}{\lambda_{1}}\log(1 - \lambda(x)))} \rightarrow (1 - \lambda(x))^{\lambda_{0}/\lambda_{1}}
\end{equation}
Substituting in:
\begin{equation}
P(0) = \frac{\lambda_{0}}{\lambda_{1}}\int^{1}_{0} e^{-\lambda(x)}(1 - \lambda(x))^{\lambda_{0}/\lambda_{1}-1}.d\lambda(x)
\end{equation}
Wolfram Alpha tells me this bizarre integral is equal to:
\begin{equation}
P(0) = \left[ \frac{\lambda_{0}}{\lambda_{1}} \frac{(1-\lambda(x))^{\lambda_{0}/ \lambda_{1}}(\lambda(x)-1)^{-\lambda_{0} / \lambda_{1}} \Gamma\left(\frac{\lambda_{0}}{\lambda_{1}}, \lambda(x)-1\right)}{e} \right]^{1}_{0}
\end{equation}
However, the same integral, this time with $P(k)$, yields something that Wolfram Alpha cannot integrate:
\begin{equation}
P(k) = \frac{\lambda_{0}}{\lambda_{1}k!}\int^{1}_{0} \lambda(x)^{k} e^{-\lambda(x)}(1 - \lambda(x))^{\lambda_{0}/\lambda_{1}-1}.d\lambda(x)
\end{equation}
Of course, I do not assume this integral is unsolvable because Wolfram Alpha says so. However, I do not know how to proceed.
Can anyone offer some commentary on how to proceed with this integral, assuming my method is correct thus far?
 A: An idea is to use the integral $I(a)=\int_0^{1}e^{-au}(1-u)^{\gamma-1}du$ to generate the others by taking derivatives. It is easy to see that with $\gamma=\lambda_0/\lambda_1$
$$P(k)=\gamma \frac{(-1)^k}{k!}I^{(k)}(1)$$
Then, calculating the generating function of the $P(k)$'s is trivial:
$$\sum_{k=0}^\infty P(k)x^k=\gamma I(1-x)$$
After some manipulations the function can be rewritten in the form
$$I(x)=\frac{e^{-x}}{x^\gamma}\int_0^{x} e^t t^{\gamma-1}dt :=e^{-x}\Delta(x):=\frac{e^{-x}}{x^\gamma}K(x)$$
This allows probabilities of arbitrary order to be computed, noting that
$$I^{(n)}(1)=\frac{1}{e}\sum_{m=0}^n(-1)^m{n\choose m} \Delta^{(m)}(1)$$
$$\Delta^{(n)}(1)=K(1)C_n(-\gamma)+e\sum_{m=0}^{n-1}{n\choose m+1}\sum_{l=0}^m{m \choose l}C_l(\gamma-1)C_{n-m-1}(-\gamma)$$
Here we have also defined $(x^{a})^{(n)}=C_{n}(a)x^{a-n}$ for brevity. This expression is relatively cumbersome to parse, but clear; all higher order probabilities $P(k)$ can be expressed in terms of $P(0)$ plus some other $\gamma$-dependent
polynomials of degree $k-1$. Also, the advantage of this expression is that only a finite number of terms of order $\mathcal{O}(k^3)$ need to be computed for each $k$.
There is at least one other interesting infinite series representation for the probabilities, obtained by directly expanding $e^{-\lambda(x)}$ in the last equation given in the OP, which allows an expression in terms of sums of beta functions. The final form is
$$P(k)=
\frac{\Gamma(\gamma+1)}{k!}\sum_{m=0}^\infty\frac{(-1)^m}{m!}\frac{\Gamma(k+m+1)}{\Gamma(\gamma+k+m+1)}$$
which after some manipulations can be expressed in terms of a low order hypergeometric function
$$P(k)=\frac{1}{(\gamma)_k} {}_1F_1(k;\gamma+k, -1)$$
where as usual $(\gamma)_k=(\gamma+1)...(\gamma+k)$.
