How to numerically integrate expression using its log transform

Question

I am trying to use numerical integration to compute the integral $\int_0^1 f(\rho) d\rho$ where $f(\rho)$ is calculated as $\ln f(\rho)$ for improved numerical precision and stability.

At the moment I am exponentiating and then integrating, i.e., numerically solving $\int_0^1 \exp(\ln f(\rho)) d\rho$

Is there another method that can be used so I don't have to exponentiate until after the integral is calculated?

I know this is not a programming board but for reference I am using scipy.integrate.quad in Python at the moment.

Many thanks

EDIT: For further information, $f(\rho)$ is the following:

$f(\rho) = Bin(x_1;N_1,f_1(\rho)) \times Bin(x_2;N_2,f_2(\rho)) \times Bin(x_3;N_3,f_3(\rho))$

CORRECTION:

$f(\rho) = Bin(x_1;N_1,f_1(\rho)) \times Bin(x_2;N_2,f_2(\rho)) \times Bin(x_3;N_3,f_3(\rho)) \times Beta(\rho;\alpha,\beta)$

where $Bin(x;N,f) = \binom{N}{x} f^x (1-f)^{N-x}$

The functions $f_1,f_2,f_3$ are linear with respect to $\rho$, e.g., $f_1(\rho) = \rho(1-\pi)$ where $\pi$ held constant.

The method of integration I believe is QAGS from Fortran QUADPACK (global adaptive quadrature based on 21-point Gauss–Kronrod quadrature within each subinterval, with acceleration by Peter Wynn's epsilon algorithm)

Answer 1

One quite simple way (in python) to perform numerical integration, when using the log of the function, is to adapt the trapezium rule. An example would be:

def logtrapz(lnf, x):
  from scipy.misc import logsumexp
  deltas = np.log(np.diff(x))
  return -np.log(2.) + logsumexp([logsumexp(lnf[:-1]+deltas), logsumexp(lnf[1:]+deltas)])

where lnf is the log of the function to be integrated evaluated at the points x. This is obviously not particularly sophisticated and provides no adaptation or estimate of the error on the result, although you could try different x spacings and see how the result is affected.

For something a bit more sophisticated, here I've adapted the GSL gsl_integration_qag function to work when supplied with the log of the function requiring integration.

Answer 2

It seems to me that you may not even need to do this numerically in this fashion. Let's write out your integral as I understand it.

$$\binom{N_1}{x_1} \binom{N_2}{x_2} \binom{N_3}{x_3} \int_0^1 d\rho \, (a_1 \rho)^{x_1} (1-a_1 \rho)^{N_1-x_1} \, (a_2 \rho)^{x_2} (1-a_2 \rho)^{N_2-x_2} \, (a_3 \rho)^{x_3} (1-a_3 \rho)^{N_3-x_3}$$

which may be rewritten as

$$\binom{N_1}{x_1} \binom{N_2}{x_2} \binom{N_3}{x_3} a_1^{x_1} a_2^{x_2} a_3^{x_3} \int_0^1 d\rho \, \rho^{x_1+x_2+x_3} (1-a_1 \rho)^{N_1-x_1} (1-a_2 \rho)^{N_2-x_2} (1-a_3 \rho)^{N_3-x_3}$$

Now this integral may be evaluated numerically with little fuss, depending on the magnitude of the exponents (if they get too large, then you may have some sampling issues to deal with, but nothing that cannot be predicted). Further, if one of the $a$'s (say, $a_1$)are equal to $1$, then there is an analytical expression for the integral (see Gradsteyn & Ryzhik, 3.211).