# Suggestions for numerical integral over Pólya Distribution

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)$

• Can you give a typical/concrete example integral? For example, giving values for n_i, N, a, and the (other terms)? Feb 12 '12 at 9:39
• I certainly will!
– yep
Feb 12 '12 at 21:39
• @sydeulissie: did you solve the problem (since you did not post a typical example)? Aug 8 '13 at 14:40
• This is inaccurate (I lose essentially all detail about the distribution except its peaks), but at least produces a number.  I don't understand why this should be a problem. The result in a bayesian approach is always dominated by the peaks (think of occam's razor). Local features will give you a neglectible contribution to the final probabilities.
– Bort
Oct 27 '16 at 11:46