Multiple integral over two dependent beta distributions

Question

I want to evaluate this multiple integral: $$ \iiint\limits_{ \sum_{i=1}^4 x_i=1,\ \ \ x_1,\, x_2,\,x_3,\, x_4\, \ge \,0 } (x_1+x_2)^{N1} (x_3+x_4)^{N2} (x_1+x_3)^{N3} (x_2+x_4)^{N4} \, dX$$

For context, this problem comes up when computing the marginal likelihood of data consisting of two dependent binary variables $(a,b)$. The conditional probability distribution for $B$ given $A$ is: \begin{align} P(B=0\mid A=0)=x_1+x_2; & & P(B=1\mid A=0)=x_3+x_4 \\ P(B=0\mid A=1)=x_1+x_3; & & P(B=0\mid A=1)=x_2+x_4 \end{align}

The marginal likelihood turns into this integral, where $N_i$ are the counts of $B\mid A$ in the data.

This looks like a dirichlet integral of type 1, but the problem is that I am unable to separate out the integral into integrals over fewer variables.

I have tried:

1. Use dirac delta function to remove the simplex and then use laplace transform to find a closed form solution. However, it does not work here because the integrals are not separable.
2. Using the variable substitution for beta variables as in this answer.
3. Interpret the equation as product of two Beta functions: $B(N_1 +1, N_2 +1), B(N_3 +1, N_4 +1)$.
   However, I am stuck because $x_i$ are shared across the terms.

Any ideas?

Answer 1

First observe this simplification: \begin{align} & \int_0^1 \int_0^{1-x_1} \int_0^{1-x_1-x_2} (x_1+x_2)^{N_1}(1-x_1-x_2)^{N_2} (x_1+x_3)^{N_3} (1-x_1-x_3)^{N_4} \,dx_3 \,dx_2 \,dx_1 \\[12pt] = {} & \int_0^1 \int_0^{1-x_1}\left( (x_1+x_2)^{N_1}(1-x_1-x_2)^{N_2} \int_0^{1-x_1-x_2} (x_1+x_3)^{N_3} (1-x_1-x_3)^{N_4} \,dx_3\right) \,dx_2 \,dx_1 \end{align} In the innermost integral, $x_1$ and $x_2$ remain constant as $x_3$ goes from $0$ to $1-x_1-x_2$; hence the step above can be done.

Now look at the inner integral: $$ \int_0^{1-x_1-x_2} (x_1+x_3)^{N_3} (1-x_1-x_3)^{N_4} \,dx_3 = \int_{x_1}^{1-x_2} u^{N_3} (1-u)^{N_4} \, du. $$ This is a doubly incomplete Beta function. By expanding via the binomial theorem and integrating term by term, you can write it as a sum, but I'm not optimistic about a closed form.

This last integral is of course a polynomial function of $x_1$ and $x_2$.