# Problem

Is it possible to simplify/rewrite the following expression, preferably without explicit sums, such that it can be computed without numerical issues when the $n_*$ are in the range of thousands?

$$\left( \prod_{i,j} {n_{ij} \choose n_{ij0}} \right) \sum_{k_{11}}^{n_{111}} {n_{111} \choose k_{11}}(-1)^{k_{11}} \dots \sum_{k_{33}}^{n_{331}} {n_{331} \choose k_{33}}(-1)^{k_{33}} \prod_i \frac{1}{1 + \sum_j n_{ij0} + k_{ij}} \prod_j \frac{1}{1 + \sum_i n_{ij0} + k_{ij}}$$ (There are 9 sums in total for $k_{11},k_{12},k_{13},k_{21}...k_{33}$, $i$ and $j$ always range from 1 to 3.)

For example by rewriting it in terms of digamma functions, like in this simplified problem: $${m+n \choose m} \sum_{i=0}^n {n \choose i}(-1)^i \frac{1}{(1+m+i)^2} = \frac{\psi(m+n+2) - \psi(m+1)}{m+n+1}$$ where $\psi(x)$ is the digamma function. Any characterization or helpful reference would be much appreciated.

# Background

I am working with the following integral $$\int_0^1 \dots \int_0^1 \prod_{i,j \in \{1,2,3\}^2} {n_{ij} \choose n_{ij0}} (p_iq_j)^{n_{ij0}}(1-p_iq_j)^{n_{ij1}}f(p_i;\alpha_i, \beta_i)f(q_j;\alpha_j,\beta_j)dp_idq_j$$ where $f(p_i;\alpha_i,\beta_i)$ is the beta distribution for $p_i$ with parameters $\alpha_i$ and $\beta_i$, and correspondingly for $f(q_j;\alpha_j,\beta_j)$. The $n_*$ variables are integer constants in the range of thousands. The $i,j$ indices always range from $1$ to $3$. The expression itself is the model posterior for a 3x3 matrix of binomials with a certain parameter structure. Specifically, each cell contains two counts, successes and failures, and the success probability for cell $i,j$ is $p_iq_j$. The integral above is the model posterior of this model.

By using the binomial theorem and the definition of the beta function I can solve the integral and get (splitted into two lines): $$\left( \prod_{i,j} {n_{ij} \choose n_{ij0}} \right) \sum_{k_{11}}^{n_{111}} {n_{111} \choose k_{11}}(-1)^{k_{11}} \dots \sum_{k_{33}}^{n_{331}} {n_{331} \choose k_{33}}(-1)^{k_{33}}\cdot$$ $$\prod_i \frac{B(\alpha_i + \sum_j n_{ij0} + k_{ij}, \beta_i)}{B(\alpha_i,\beta_i)} \prod_j \frac{B(\alpha_j + \sum_i n_{ij0} + k_{ij}, \beta_j)}{B(\alpha_j,\beta_j)}$$

where $B(x,y)$ is the Beta function.

For my own purposes it is enough if I can simplify it in the case $\alpha_i = \beta_i = \alpha_j = \beta_j = 1$, that is the expression above: $$\left( \prod_{i,j} {n_{ij} \choose n_{ij0}} \right) \sum_{k_{11}}^{n_{111}} {n_{111} \choose k_{11}}(-1)^{k_{11}} \dots \sum_{k_{33}}^{n_{331}} {n_{331} \choose k_{33}}(-1)^{k_{33}} \prod_i \frac{1}{1 + \sum_j n_{ij0} + k_{ij}} \prod_j \frac{1}{1 + \sum_i n_{ij0} + k_{ij}}$$

A potentially helpful observation is that it can be written as the finite difference $$\left( \prod_{i,j} {n_{ij} \choose n_{ij0}} \right)\Delta_{k_{11}}^{n_{111}} \dots \Delta_{k_{33}}^{n_{331}} \prod_i \frac{B(\alpha_i + \sum_j n_{ij0} + k_{ij}, \beta_i)}{B(\alpha_i,\beta_i)} \prod_j \frac{B(\alpha_j + \sum_i n_{ij0} + k_{ij}, \beta_j)}{B(\alpha_j,\beta_j)}$$ where $\Delta_k^n$ is the finite difference operator for the $n$:th difference with respect to $k$.

# Numerical issues

Even though the integral is solved it is hard to compute numerically due to the large binomial coefficients. The issue is that the terms of the sum have to cancel properly for the expression to be between 0 and 1. For big $n$ this does not happen due to the addition of very small and very large floating point terms. It seems like both the lower and higher order bits are important for the cancellation to happen properly. If I use 2000 bit floats I get a reasonable answer (which I also can verify for small $n_*$), this however takes too long to compute.

I have also tried (also in combination):

• Formulating it as a recursion based on finite differences and thus getting rid of explicit binomial coefficients also failed.
• Adding terms of similar order of magnitude first.
• Various summation algorithms like Kahan summation.
-