Integrate Gamma pdf with respect to shape parameter alpha Is there a trick to integrating a Gamma pdf with respect to the shape parameter alpha? I have yet to come up with a way to do it.
 A: It is unusual to encounter any situation where you want to integrate a gamma density with respect to the shape parameter.  (The only instance I can think of is in a Bayesian model with an improper infinite uniform prior on the shape, which is a weird model.)  In any case, the (definite) integral of a gamma density is a scaled version of the following:
$$I(y) \equiv \int \limits_0^\infty \frac{y^{\alpha-1}}{\Gamma(\alpha)} d\alpha.$$
There is no closed form for this integral so you would need to use numerical integration to obtain its value.  The integrand rapidly becomes small for large values of $\alpha$, so we can approximate the integral using a finite upper bound to replace the infinite upper bound.
Numerical approximation of the integral can be done using any numerical technique for numerical integration, but here is an approximating function using Simpson's second rule (which uses cubic interpolation of the integrand).  We select some large upper bound $D$ for the integral (where the integrand has become small) and some large number $N$ (with $N \text{ mod } 3 = 1$) for the number of intervals in the numerical approximation.  We then define the approximating function:
$$\hat{I}(y, D, N) \equiv \frac{3}{8} \cdot \frac{D}{N}  
\begin{bmatrix}
\frac{y^{\alpha_0-1}}{\Gamma(\alpha_0)} 
+ 3 \frac{y^{\alpha_1-1}}{\Gamma(\alpha_1)} + 3\frac{y^{\alpha_2-1}}{\Gamma(\alpha_2)} + 2 \frac{y^{\alpha_3-1}}{\Gamma(\alpha_3)} \quad \quad \quad \quad \quad \quad \\[6pt]
\quad + \text{ } 3 \frac{y^{\alpha_4-1}}{\Gamma(\alpha_4)} + 3\frac{y^{\alpha_5-1}}{\Gamma(\alpha_5)} + 2 \frac{y^{\alpha_6-1}}{\Gamma(\alpha_6)} \quad \quad \quad \quad \quad \\[6pt]
\quad  + \cdots \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \text{ } \text{ } \\[6pt]
\quad \quad + 3 \frac{y^{\alpha_{N-3}-1}}{\Gamma(\alpha_{N-3})} + 3\frac{y^{\alpha_{N-2}-1}}{\Gamma(\alpha_{N-2})} + 2 \frac{y^{\alpha_{N-1}-1}}{\Gamma(\alpha_{N-1})} + \frac{y^{\alpha_N-1}}{\Gamma(\alpha_N)}
\end{bmatrix},$$
where $\alpha_i \equiv i D /N$.  This approximation approaches the true integral as $D \rightarrow \infty$ and $N \rightarrow \infty$, but it should suffice to take $D$ sufficiently large so that the integrand is very small.  Convergence can be tested by looking at how much the value changes as you change $D$ and $N$.
