Numpy: How to avoid rounding errors

Question

I think I have a problem with Numpy and small numbers.

Can you help me finding a solution for the following problem:

import numpy as np

def gaussian(xx, mu=0, sigma=1):
    return 1./(np.sqrt(2*np.pi)*sigma)*np.exp(-(mu-xx)**2/(2*sigma**2))

factors = (1., 0.1, 0.01, 0.001, 0.0001)
for factor in factors:
    xx = np.linspace(5000, 5200, 1000)
    yy = 1.-(factor*xx*(1.+gaussian(xx, 5100)))
    step = xx[1] - xx[0]

    print np.sum((1.-yy/(1.-factor*xx))*step)

This code should evaluate to -1 for all of the different factors. But the output is :

1.0 -1.00019611689
0.1 -1.00196463662
0.01 -1.0200000008
0.001 -1.24390245353
0.0001 1.04081641153

So the problem is, that the smaller the factor is, the more I get into trouble because of what I think has to do with the precision of Numpy/Python.

What to do to evaluate the equation even for the small factors?

Thanks a lot for your help in advance.

I suspect that at least part of the issue is in the mathematics (or, more probably, in the translation of the mathematics to Python + NumPy code). Why do you think this should evaluate to -1? What's this based on? — Mark Dickinson
– Mark Dickinson, Commented Feb 7, 2014 at 23:08
Rewrite the yy-part to 1-ax(1+g(x)) and the then part inside the sum to 1-yy/(1-ax) = 1-((1-ax(1+g(x))/(1-ax)) = 1-(1+g(x)) = -g(x). Since g(x) is a gaussian which has its peak at 5100 and a stddev of 1 and we are integrating over a range of 200 (from 5000 to 5200) the g(x) will be 1 and thus the result will be -1. I think the mathematics is not the problem but the way Numpy/Python deals with the small numbers. — dotcs
– dotcs, Commented Feb 7, 2014 at 23:17
Nope: there's an error in your algebra. Your expression simplifies to axg(x) / (1 - ax), not to -g(x). — Mark Dickinson
– Mark Dickinson, Commented Feb 7, 2014 at 23:19
Perhaps you meant to define yy as (1 - factor*xx)*(1+gaussian(xx, 5100)) instead? — Mark Dickinson
– Mark Dickinson, Commented Feb 7, 2014 at 23:23
Thanks for the hint. When doing that the problem is solved and everything works fine. I tried to simplify a problem I had getting the code above. Obviously something went wrong. I need to check where I made the mistake ... — dotcs
– dotcs, Commented Feb 7, 2014 at 23:32

Mark Dickinson · Accepted Answer · 2014-02-07 23:40:19Z

3

This isn't a NumPy issue. Your expectation that the result should be -1 is incorrect.

You're effectively computing the definite integral of a function f(x) on a large interval around 5100. The function you're integrating simplifies to factor * x * gaussian(x) / (1 - factor * x). We can easily give a back-of-the-envelope estimate for the value of the integral for the factors you're using: the quantity factor * x / (1 - factor * x) will vary fairly slowly across the range of interest, which is about 5095 to 5105; for anything outside this range the gaussian will render the contribution negligible. So to a first approximation we can treat that quantity as constant. Then we're left with that constant times the integral of gaussian(x), which will be near enough 1. Thus the expected output should be somewhere in the region of factor * 5100 / (1 - factor * 5100). (Though that's not going to work so well when factor is close to 1 / 5100.)

So for example in the case that factor is 0.0001, factor * 5100 / (1 - factor * 5100) has value around 1.0408163265306123. This is close enough to the answer you're seeing to make it plausible that NumPy is doing more-or-less the right thing.

answered Feb 7, 2014 at 23:40

Mark Dickinson

31.2k9 gold badges92 silver badges130 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

dotcs Over a year ago

Yeah, you are perfectly correct. I wanted to simplify a problem I had and made the mistake you described above. I am really sorry about that and although I have checked my simplification I messed it up. Thanks for your effort. Shame on me and upvotes for your correction.

Collectives™ on Stack Overflow

Numpy: How to avoid rounding errors

1 Answer 1

1 Comment

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

1 Comment

Your Answer

Sign up or log in

Post as a guest

Linked

Related