Approximation to a compounded Binomial distribution I need to find an approximation, from which I can easily sample, to the following compounded Binomial distribution:
$X \sim \mathrm{Binomial}(e^{-\epsilon}, \ n)$ where $\epsilon \sim \mathrm{Gamma}(\alpha, \beta)$
I know that if $e^{-\epsilon}$ were a Beta then $X$ would be a Beta-Binomial. 
Moreover, I know that $N\ \mathrm{Beta(\alpha, N)} \rightarrow \mathrm{Gamma}(\alpha, 1)$ in distribution. 
However $e^{-\epsilon}$ is not a Gamma and I cannot use this nice approximation. $\epsilon$ has parameters such that its mode is around 0.1 so I thought I could use a Taylor expansion of $e^{-\epsilon}$, but it does not help in making my compound Binomial tractable.
Since I want to sample from this in a particle filter, I don't want to use any importance sampling (or related) because I need something quite fast, so I'm looking for a known distribution which could be a good approximation.
EDIT 1:
I think I've found an answer to my problem.
Notice that $e^{-\epsilon} \in ]0, 1]$. So it's the support of a Beta. Natural idea is to minimize the KL divergence between the distribution of $e^{-\epsilon}$ (we know its density its just an easy change of variable) and a Beta. 
If we denote $P=e^{-\epsilon}$ and $q_P(p)=\frac{\beta^\alpha}{\Gamma(\alpha)}(-\log p)^{\alpha-1}p^{\beta-1}$ its density and if we denote $q(\ \cdot \ ; a, b)$ the density of a $\mathrm{Beta}(a,b)$ we have:
\begin{equation*}
\begin{split}
D_{KL}(a,b) & = \int_{0}^{1}q_P(p)\frac{q_P(p)}{q(p; a, b)}\mathrm{d}p \\
& \propto - \int_{0}^{1}q_P(p)\log(q(p; a, b))\mathrm{d}p \\
& \propto - \int_{0}^{1}q_P(p)\log(\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}p^{a-1}(1-p)^{b-1})\mathrm{d}p \\
& \propto -\log\Gamma(a+b) + \log\Gamma(a) + \log\Gamma(b) \\ & \qquad  + (a-1)\int_{0}^{1}(-\log p)\frac{\beta^\alpha}{\Gamma(\alpha)}(-\log p)^{\alpha-1}p^{\beta-1}\mathrm{d}p \\ & \qquad - (b-1)\int_{0}^{1}\log (1-p)\frac{\beta^\alpha}{\Gamma(\alpha)}(-\log p)^{\alpha-1}p^{\beta-1}\mathrm{d}p
\end{split}
\end{equation*}
Now notice that $-\log p(-\log p)^{\alpha-1}p^{\beta-1}= (-\log p)^{\alpha}p^{\beta-1}$ is proportional to the density of $e^{-X}$ where $X \sim \mathrm{Gamma}(\alpha+1, \beta)$ and if we denote $K=\int_{0}^{1}\log(1-p)\frac{\beta^\alpha}{\Gamma(\alpha)}(-\log p)^{\alpha-1}p^{\beta-1}\mathrm{d}p$ we have:
\begin{equation*}
\begin{split}
D_{KL}(a,b) & = \propto -\log\Gamma(a+b) + \log\Gamma(a) + \log\Gamma(b) \\ 
& \qquad  + (a-1)\frac{\beta^\alpha}{\Gamma(\alpha)}\frac{\Gamma(\alpha+1)}{\beta^{\alpha+1}} - (b-1)K \\
& \propto -\log\Gamma(a+b) + \log\Gamma(a) + \log\Gamma(b) + (a-1)\frac{\alpha}{\beta} - (b-1)K
\end{split}
\end{equation*}
Now we minimize by finding a critical point:
\begin{equation}
\frac{\partial }{\partial a}D_{KL} = -\psi^{(0)}(a+b) + \psi^{(0)}(a) + \frac{\alpha}{\beta} = 0
\end{equation}
\begin{equation}
\frac{\partial }{\partial b}D_{KL} = -\psi^{(0)}(a+b) + \psi^{(0)}(b) - K = 0
\end{equation}
where $\psi^{(0)}$ is the digamma function.
$K$ is easily calculated numerically and so is the system of two equations above
I did a numerical example with  $\alpha = \beta = 10$ and found $a = 6.739281$ and $b=10.736471$. I sampled 10000 times from $e^{-\epsilon}$ and the Beta with coefficient $a$ and $b$ and drew the following QQ-plot which seems to indicate that the approximation is accurate and I can use a Beta-Binomial to sample from instead of from $\mathrm{Binomial}(e^{-\epsilon}, n)$:

 A: The density of $\pi = e^{-\epsilon}$ where $\epsilon \sim \operatorname{Gamma}(a,b)$ is parametrized by rate, is given by $$f_\pi(p) = \frac{b^a (-\log p)^{a-1} p^{b-1}}{\Gamma(a)}, \quad 0 < x < 1.$$  Therefore $$\begin{align*} \Pr[X = x] &= \int_{p=0}^1 \Pr[X = x \mid \pi = p] f_\pi(p) \, dp \\ &= \binom{n}{x} \frac{b^a}{\Gamma(a)} \int_{p=0}^1  p^{x+b-1} (1-p)^{n-x} (-\log p)^{a-1} \, dp \\ &= \binom{n}{x} \frac{b^a}{\Gamma(a)} \int_{v=0}^\infty e^{-(x+b-1)v} (1-e^{-v})^{n-x} v^{a-1} e^{-v} \, dv. \end{align*}$$  At this point we can see that if $n, x$ are nonnegative integers (which of course they are), we could expand the integrand, integrate term by term, and obtain a finite linear combination of gamma functions in terms of $a$ and $b$.  The resulting formula is lengthy but closed form and programmable.
$$\begin{align*} \Pr[X = x] &= \binom{n}{x} \frac{b^a}{\Gamma(a)} \sum_{k=0}^{n-x} \binom{n-x}{k} (-1)^k \int_{v=0}^\infty v^{a-1} e^{-(x+b+k)v} \, dv \\ &= \binom{n}{x} \frac{b^a}{\Gamma(a)} \sum_{k=0}^{n-x} \binom{n-x}{k} (-1)^k \frac{\Gamma(a)}{(x+b+k)^a} \\ &= \binom{n}{x} b^a \sum_{k=0}^{n-x} \binom{n-x}{k} \frac{(-1)^k}{(x+b+k)^a}. \end{align*}$$
To utilize this result, it is probably best to choose your parameters $n, a, b$ in advance and compute the probability mass function for each $x = 0, 1, 2, \ldots, n$, provided that $n$ is not huge.
