scipy.integrate.quad precision on big numbers
<p>I try to compute such an integral (actually cdf of exponential distribution with its pdf) via <code>scipy.integrate.quad()</code>:</p>

<pre><code>import numpy as np
from scipy.integrate import quad

def g(x):
    return .5 * np.exp(-.5 * x)

print quad(g, a=0., b=np.inf)
print quad(g, a=0., b=10**6)
print quad(g, a=0., b=10**5)
print quad(g, a=0., b=10**4)
</code></pre>

<p>And the result is as follows:</p>

<pre><code>(1.0, 3.5807346295637055e-11)
(0.0, 0.0)
(3.881683817604194e-22, 7.717972744764185e-22)
(1.0, 1.6059202674761255e-14)
</code></pre>

<p>All the attempts to use a big upper integration limit yield an incorrect answer though the usage of <code>np.inf</code> solves the problem.</p>

<p>Similiar case is discussed in <a href="https://github.com/scipy/scipy/issues/5428" rel="nofollow">scipy issue #5428 at GitHub</a>.</p>

<p>What should I do to avoid such an error in integrating other density functions?</p>

<p>I believe the issue is due to <code>np.exp(-x)</code> quickly becoming very small as <code>x</code> increases, which results in evaluating as zero due to limited numerical precision. For example, even for <code>x</code> as small as <code>x=10**2*</code>, <code>np.exp(-x)</code> evaluates to <code>3.72007597602e-44</code>, whereas <code>x</code> values of order <code>10**3</code> or above result in <code>0</code>.</p>

<p>I do not know the implementation specifics of <code>quad</code>, but it probably performs some kind of sampling of the function to be integrated over the given integration range. For a large upper integration limit, most of the samples of <code>np.exp(-x)</code> evaluate to zero, hence the integral value is underestimated. (Note that in these cases the provided absolute error by <code>quad</code> is of the same order as the integral value which is an indicator that the latter is unreliable.)</p>

<p>One approach to avoid this issue is to restrict the integration upper bound to a value above which the numerical function becomes very small (and, hence, contributes marginally to the integral value). From your code snipet, the value of <code>10**4</code> appears to be a good choice, however, a value of <code>10**2</code> also results in an accurate evaluation of the integral.</p>

<p>Another approach to avoid numerical precision issues is to use a module that performs computation in <em>arbitrary</em> precision arithmetic, such as <code>mpmath</code>. For example, for <code>x=10**5</code>, <code>mpmath</code> evaluates <code>exp(-x)</code> as follows (using the native <code>mpmath</code> exponential function) </p>

<pre><code>import mpmath as mp
print(mp.exp(-10**5))
</code></pre>

<blockquote>
  <p><code>3.56294956530937e-43430</code></p>
</blockquote>

<p>Note how small this value is. With the standard hardware numerical precision (used by <code>numpy</code>) this value becomes <code>0</code>. </p>

<p><code>mpmath</code> offers an integration function (<code>mp.quad</code>), which can provide an accurate estimate of the integral for arbitrary values of the upper integral bound.</p>

<pre><code>import mpmath as mp

print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, mp.inf]))
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**13]))
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**8]))
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**5]))
</code></pre>

<blockquote>
<pre><code>1.0
0.999999650469474
0.999999999996516
0.999999999999997
</code></pre>
</blockquote>

<p>We can also obtain even more accurate estimates by increasing the precision to, say, <code>50</code> decimal points (from <code>15</code> which is the standard precision)</p>

<pre><code>mp.mp.dps = 50; 

print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, mp.inf]))
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**13]))
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**8]))
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**5]))
</code></pre>

<blockquote>
<pre><code>1.0
0.99999999999999999999999999999999999999999829880262
0.99999999999999999999999999999999999999999999997463
0.99999999999999999999999999999999999999999999999998
</code></pre>
</blockquote>

<p>In general, the cost for obtaining this accuracy is an increased computation time. </p>

<p>P.S.: It goes without saying that if you are able to evaluate your integral analytically in the first place (e.g., with the help of <code>Sympy</code>) you can forget all the above. </p>

python numpy scipy calculus