这个问题来自贝叶斯统计建模项目。为了使用我的模型进行计算，我需要执行积分，其中被积函数的一部分是“Pólya”或“Dirichlet-Multinomial”分布，

$$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!}$$

在哪里$n_i$和$N = \sum_{i=1}^K n_i$是整数，$n = \left(n_1, n_2, \dots, n_K\right)$， 和$\alpha > 0$. 我想计算的积分，$\int_0^\infty (\text{other terms})p(n|\alpha) d\alpha$, 适用于小型$N$，但我尝试过的求积方法（在 MATLAB 中）分解为$N$变大。我没试过蒙特卡洛；一个准确、快速的求积方法对我的项目来说非常好。

目前，“最好”的方法是$N$大是计算$\log[p(n|\alpha)]$在 alpha 的网格上，归一化和取幂。这是不准确的（我基本上失去了关于分布的所有细节，除了它的峰值），但至少产生了一个数字。

我将不胜感激有关改进此计算的任何建议，或指向不同算法/方法或现有软件的指针。

编辑：我也许应该补充一点，我对$p(n|\alpha)$, 通过计算执行$\log p(n|\alpha)$使用一些精心编写的代码来计算$\log \Gamma(x)$对于大$x$, 似乎没有引起任何问题。

编辑2：此外，“大”值将是$N\sim 10^8$, 最大$n_i\sim 10^5$，以及许多小的值$n_i$. 其他项在数值上表现良好。作为大致适当的尾部行为的简化，您可以采取

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