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

  • $\begingroup$ Can you give a typical/concrete example integral? For example, giving values for n_i, N, a, and the (other terms)? $\endgroup$ Feb 12, 2012 at 9:39
  • $\begingroup$ I certainly will! $\endgroup$
    – yep
    Feb 12, 2012 at 21:39
  • $\begingroup$ @sydeulissie: did you solve the problem (since you did not post a typical example)? $\endgroup$
    – GertVdE
    Aug 8, 2013 at 14:40
  • $\begingroup$ 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. $\endgroup$
    – Bort
    Oct 27, 2016 at 11:46

1 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.