title: Dealing with a precision error when generating probability distribution in python
body: <p>In the context of a model I'm generating, I at one point need to generate a probability distribution from a an array of real numbers. I'll leave out the relevant details, but essentially have a function (we'll just call it "f" for now), that generates an array of n floats:</p>
  
  <pre><code>arr = [value_1, value_2, ..., value_n]
  </code></pre>
  
  <p>Now, these values are <strong>proportional</strong> to probabilities I next need to use in a multinomial sampling procedure, so the obvious approach is simply this:</p>
  
  <pre><code>result = np.random.multinomial(number_of_samples,arr/arr.sum())
  </code></pre>
  
  <p>But this (sometimes) doesn't work! Basically the sum of arr/arr.sum() ends up being greater than 1. In principle this should be mathematically impossible, but I'm assuming this boils down to a floating-point precision issue. Here's a trivial example of how this can happen:</p>
  
  <pre><code>In [58]: arr = np.array([1/20.]*20)
  
  In [59]: arr/arr.sum()
  Out[59]:
  array([ 0.05,  0.05,  0.05,  0.05,  0.05,  0.05,  0.05,  0.05,  0.05,
      0.05,  0.05,  0.05,  0.05,  0.05,  0.05,  0.05,  0.05,  0.05,
      0.05,  0.05])
  
  In [60]: (arr/arr.sum()).sum()
  Out[60]: 1.0000000000000002
  </code></pre>
  
  <p>So long story short, my question is how best to deal with this. I can cheat by simply adding a very small number to the sum, i.e.:</p>
  
  <pre><code>probs = arr / (arr.sum()+0.000001)
  </code></pre>
  
  <p>But this is really hackish, and I fear it may introduce further unwanted precision issues. Is there a better solution?</p>