[FIXED] How to avoid roundoff errors in numpy.random.choice?

Issue

Say x_1, x_2, ..., x_n are n objects and one wants to pick one of them so that the probability of choosing x_i is proportional to some number u_i.

x, u = np.array([x_1, x_2, ..., x_n]), np.array([u_1, ..., u_n])
np.random.choice(x, p = u/np.sum(u))

There are round-off errors due to finite precision arithmetic. Should one be worried about getting an error because the sum of probabilities is not exactly 1? Is there a standard solution to this problem?

Solution

After reading the answer https://stackoverflow.com/a/60386427/6087087 to the question pointed by @Pychopath, I have found the following solution, inspired by the documentation of numpy.random.multinomial https://docs.scipy.org/doc/numpy-1.15.0/reference/generated/numpy.random.multinomial.html

Say p is the array of probabilities which may not be exactly 1 due to roundoff errors, even if we normalized it with p = p/np.sum(p). This is not rare, see the comment by @pd shah at the answer https://stackoverflow.com/a/46539921/6087087.

Just do

p[-1] = 1 - np.sum(p[0:-1])
np.random.choice(x, p = p)

And the problem is solved! The roundoff errors due to subtraction will be much smaller than roundoff errors due to normalization. Moreover, one need not worry about the changes in p, they are of the order of roundoff errors.

Answered By - Fırat Kıyak