CODEWITHSUNDEEP
pythonnumpyfloating-pointnumbersscipy
abcd
abcd

Reputation: 10771

Tricking numpy/python into representing very large and very small numbers

I need to compute the integral of the following function within ranges that start as low as -150:

import numpy as np
from scipy.special import ndtr

def my_func(x):
    return np.exp(x ** 2) * 2 * ndtr(x * np.sqrt(2))

The problem is that this part of the function

np.exp(x ** 2)

tends toward infinity -- I get inf for values of x less than approximately -26.

And this part of the function

2 * ndtr(x * np.sqrt(2))

which is equivalent to

from scipy.special import erf

1 + erf(x)

tends toward 0.

So, a very, very large number times a very, very small number should give me a reasonably sized number -- but, instead of that, python is giving me nan.

What can I do to circumvent this problem?

Upvotes: 9

Views: 2392

Answers (3)

Simon Byrne
Simon Byrne

Reputation: 7874

There already is such a function: erfcx. I think erfcx(-x) should give you the integrand you want (note that 1+erf(x)=erfc(-x)).

Upvotes: 4

ali_m
ali_m

Reputation: 74222

I think a combination of @askewchan's solution and scipy.special.log_ndtr will do the trick:

from scipy.special import log_ndtr

_log2 = np.log(2)
_sqrt2 = np.sqrt(2)

def my_func(x):
    return np.exp(x ** 2) * 2 * ndtr(x * np.sqrt(2))

def my_func2(x):
    return np.exp(x * x + _log2 + log_ndtr(x * _sqrt2))

print(my_func(-150))
# nan

print(my_func2(-150)
# 0.0037611803122451198

For x <= -20, log_ndtr(x) uses a Taylor series expansion of the error function to iteratively compute the log CDF directly, which is much more numerically stable than simply taking log(ndtr(x)).

Update

As you mentioned in the comments, the exp can also overflow if x is sufficiently large. Whilst you could work around this using mpmath.exp, a simpler and faster method is to cast up to a np.longdouble which, on my machine, can represent values up to 1.189731495357231765e+4932:

import mpmath

def my_func3(x):
    return mpmath.exp(x * x + _log2 + log_ndtr(x * _sqrt2))

def my_func4(x):
    return np.exp(np.float128(x * x + _log2 + log_ndtr(x * _sqrt2)))

print(my_func2(50))
# inf

print(my_func3(50))
# mpf('1.0895188633566085e+1086')

print(my_func4(50))
# 1.0895188633566084842e+1086

%timeit my_func3(50)
# The slowest run took 8.01 times longer than the fastest. This could mean that
# an intermediate result is being cached  100000 loops, best of 3: 15.5 µs per
# loop

%timeit my_func4(50)
# The slowest run took 11.11 times longer than the fastest. This could mean
# that an intermediate result is being cached  100000 loops, best of 3: 2.9 µs
# per loop

Upvotes: 5

Salvador Dali
Salvador Dali

Reputation: 222761

Not sure how helpful will this be, but here are a couple of thoughts that are too long for a comment.

You need to calculate the integral of 2 \cdot e^{x^2} \cdot f(\sqrt{2}x), which you correctly identified would be e^{x^2}*(1 + erf(x)). Opening the brackets you can integrate both parts of the summation.

enter image description here

Scipy has this imaginary error function implemented

The second part is harder:

enter image description here

This is a generalized hypergeometric function. Sadly it looks like scipy does not have an implementation of it, but this package claims it does.

Here I used indefinite integrals without constants, knowing the from to values it is clear how to use definite ones.

Upvotes: 2

Related Questions