[SciPy-User] [SciPy-user] Speeding up scipy.integrate.quad
Tue Oct 25 09:54:51 CDT 2011
On Tue, Oct 25, 2011 at 4:48 AM, tomrichardson <
> I'm currently trying to implement a Metropolis Hastings algorithm for which
> numerical integration has to be called to obtain the acceptance probability
> of any given trial in the iteration. I'm currently using the
> scipy.integrate.quad function but it's too slow.Only part of the integrand
> has an analytical form - in fact it depends on an interpolation function
> that uses data in two numpy arrays. I've plotted the integrand and it
> diverges near the lower limit of integration, there is a term
> 1/(r^2-R^2) where r is the the dummy variable over which I am integrating
> and R is the lower limit of integration [I've already added a small
> increment to R to avoid a zero division error]
The integral of 1/(r^2 - R^2) from r=R to r=a>R is divergent (i.e. the area
under that curve is infinite). Does some other term in your integrand
cancel the singularity at r=R to give a convergent integral? If so, can you
rewrite your integrand so that these terms are handled separately (ideally
analytically), and only use quad for the well-behaved part of the integrand?
> Things I've tried:
> -I've implemented all of the other QUADPACK integration methods in scipy
> which are all slower
> -reduced the epsabs and epsrel limits
> -split the interval into subranges that weight the divergence
> Are there any special purpose integrators that anyone can recommend that
> deal with this integrand and compete favourably with quad for speed?
> I'm not familiar with the C language but I have Cython and Weave installed
> as part of the Enthought Python pack. Would implementing the numerical
> integration in the C compiler increase the speed dramatically? If so how do
> I introduce the numpy arrays in the cdef terminology for Cython / are there
> any C++ recipes that use a sampled data input that I can implement
> (relatively pain free as a C++ noob).
> View this message in context:
> Sent from the Scipy-User mailing list archive at Nabble.com.
> SciPy-User mailing list
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the SciPy-User