Suggestions for numerical integral over Pólya Distribution

Question

This problem arises from a Bayesian statistical modeling project. In order to compute with my model, I need to perform an integration in which part of the integrand is the "Pólya" or "Dirichlet-Multinomial" Distribution,

$$p(n\mid \alpha) = \frac{(N!) \Gamma(K\alpha)}{\Gamma(\alpha)^K \Gamma ( N + K\alpha)} \prod_{i=1}^K \frac{\Gamma(n_i + \alpha)}{ n_i!}$$

where $n_i$ and $N = \sum_{i=1}^K n_i$ are integers, $n = \left(n_1, n_2, \dots, n_K\right)$, and $\alpha > 0$. The integral I wish to compute, $\int_0^\infty (\text{other terms})p(n|\alpha) d\alpha$, works well for small $N$, but the quadrature methods I've attempted (in MATLAB) break down as $N$ becomes large. I haven't tried Monte Carlo; an accurate, fast quadrature method would be very nice for my project.

Currently, the "best" method when $N$ is large is to compute $\log[p(n|\alpha)]$ over a grid in alpha, normalize, and exponentiate. This is inaccurate (I lose essentially all detail about the distribution except its peaks), but at least produces a number.

I would appreciate any advice on improving this computation, or pointers to different algorithms/methods or existing software.

EDIT: I should maybe add that that my evaluation of $p(n|\alpha)$, performed by computing $\log p(n|\alpha)$ using some carefully-written code to compute $\log \Gamma(x)$ for large $x$, does not appear to be causing any problems.

EDIT 2: Additionally, "large" values would be on the order of $N\sim 10^8$, with the largest $n_i\sim 10^5$, along with many small values of $n_i$. The other terms are numerically well-behaved. As a simplification with roughly the appropriate tail behavior, you could take

$(\text{other terms}) = \exp(-\alpha)$

Answer 1

For integration use a gaussian quadrature rule when dealing with exponentials as a rule of thumb. By way of example, integrating a cubic function over a small area with simpsons rule or trapezoidal rule doesn't need a large N, but using them for an exponential function leads to huge N required to achieve convergence. So always go for an exponential based interpolant (i.e: gaussian quadrature rule). If this fails, then you need to plot the function to see how quickly it is oscillating and how it dies/rises with the area of integration.