Avoiding inaccuracy when using floating points numbers?

Question

I'm trying to use Monte Carlo integration to approximate the region under a given graph to calculate its area. For this to work, I need the calculated y_min and y_max to be accurate. So as an example I'll use the graph of sin(x) from 0 to pi. To find y_min and y_max I have the following function:

def y_range(f, x_min, x_max, n=100):
    # Step size
    h = float((x_max - x_min)) / n

    # Calculate y for n points between x_min and x_max
    y = [f(x * h) for x in range(0, n + 1)]

    # Get minimum and maximum y
    y_max = max(y)
    y_min = min(y)

    return y_min, y_max

Printing y_min and y_max gives:

y_max = 1.0
y_min = -3.21624529935e-16

I know y_min should equal 0.0, so how do I rectify this inaccuracy?

That's *really close* to zero... – jonrsharpe Jan 18 '16 at 22:09 — jonrsharpe, Jan 18 '16 at 22:09

score 1 · Answer 1 · answered Jan 20 '16 at 16:22

The underlying issue is that max cannot be derived as the result of min + 100*h, for any h. There are a couple of potential causes, but the most straightforward one is simply that the number of steps between them isn't divisible by 100.

How exactly to do it better depends on exactly how careful you want to be. The big two things you need are to interpolate between the two values (rather than based on a beginning point and a step), and to perform the interpolation itself in an accurate fashion.

The following code will produce dependable results:

def interp_at_step(a, b, i, n):
    # separate calculation of alpha and beta to avoid catastrophic cancellation
    alpha = (n-i)/n
    beta = i/n
    return a*alpha + b*beta

def y_range(f, x_min, x_max, n=100):
    # Calculate y for n points between x_min and x_max
    y = [f(interp_at_step(x_min, x_max, x, n)) for x in range(0, n + 1)]

    # Get minimum and maximum y
    y_max = max(y)
    y_min = min(y)

    return y_min, y_max

Of course, as Denys mentioned, pi can't be exactly represented. So this will alleviate errors from the interpolation, but not necessarily from the operands themselves.

score 0 · Answer 2 · answered Jan 18 '16 at 23:09

pi is irrational, so you have to accept that it can only be represented with certain precision, no matter which kind of representation you are using. If you are using float, then the precision will be in order of 10**(-16) as mentioned here. An easiest way forward then is to round your result to 15 decimal digits with round(value,15). If you need better precision than 15 decimal digits, you may indeed look into other types for representing variables, such as decimal

score -1 · Accepted Answer · edited May 23 '17 at 11:44

-1

I wouldn't use float when precision is so important. I suggest you look at decimal.Decimal. Assuming that you change x_min and x_max to Decimal, you only need to remove that float on line 3.

This post might help explain where the inaccuracy is coming from:

edited May 23 '17 at 11:44

Community

1
1

answered Jan 18 '16 at 22:22

Nabil Khoury

22
4

`decimal.Decimal` isn't going to help you exactly represent pi. – Sneftel Jan 20 '16 at 16:25

Avoiding inaccuracy when using floating points numbers?

3 Answers3