3
$\begingroup$

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)

$\endgroup$
2
  • $\begingroup$ Can you actually provide $f$ and the numerical method you are using ? Regards $\endgroup$
    – Amzoti
    Commented Jul 25, 2013 at 12:50
  • $\begingroup$ Have added further details to question. Hope that helps :) $\endgroup$
    – tristan
    Commented Jul 25, 2013 at 13:08

2 Answers 2

2
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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

$\endgroup$
3
  • $\begingroup$ Hi Ron, thanks. Sorry I wasn't clear enough in my clarifications! I meant that $f_1(\rho)=a+b\rho$ . Could you provide an updated answer? $\endgroup$
    – tristan
    Commented Jul 25, 2013 at 21:46
  • $\begingroup$ Notice also that there is a beta pdf to include in the integrand... I'd vote up because it's a helpful start but not enough privilege at the moment! $\endgroup$
    – tristan
    Commented Jul 26, 2013 at 6:17
  • $\begingroup$ @tristan: sorry about taking so long. The idea does not change then: the integrand is essentially a product of six binomials raised to various powers. Again, the difficult binomial coefficients are pulled out of the integral and you need not worry; meanwhile, the integral may be directly attacked. $\endgroup$
    – Ron Gordon
    Commented Jul 26, 2013 at 6:21

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .