Create a probability list which sums up to one

Question

I use the numpy.random.choice method to draw a sample from a list.

The entries in the list are short dicts and look like

{'name': 'SomeName', 'strength': 0.75},
{'name': 'SomeName', 'strength': 0.25},
...

I want to use the 'strength' value (between 0 and 1) as a probability indicator.

To use drawResults = choice(data, drawCount, p=resultProb)

I have to create the resultProb list which has to sum up to 1

So I came up with a function:

def makeProbabilities(data) -> list:
    probResult = []

    totalScore = sum(item['strength'] for item in data) # calculate in sum of all 'strength' values

    for item in data:
        if item['strength'] > 0:
            value = (100.0/totalScore)*item['strength'] #how much is this strength in relation to the total
            value2 = (1.0/100)*value                    #how much is the result above in relation to 1

            probResult.append(value2)
        else:
            probResult.append(0.0)

    return probResult

That seems to work, it just have a very small rounding error on the result (the sum is like 1.0000000001) but the numby.choice method does accept it.

But I have the strong impression that this solution is a bit awkward, un-pythonique and perhaps not very efficient on large data sets.

I am just discovering python so I am a bit lost with all the information concerning this language.

Any feedback on that is very welcome.

Answer 1

Binary floating point numbers cannot represent most numbers exactly, a small error cannot be avoided when computing the sum of probabilities:

>>> probs = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
>>> sum(probs)
0.9999999999999999

For more information, see

• Is floating point math broken? on Stack Overflow,
• What Every Computer Scientist Should Know About Floating-Point Arithmetic.

Actually one can see in the source code of mtrand.pyx that this rounding error is cared for in the choice() function:

        atol = np.sqrt(np.finfo(np.float64).eps)
        # ...
        p = <ndarray>PyArray_ContiguousFromObject(p, NPY_DOUBLE, 1, 1)
        pix = <double*>PyArray_DATA(p)
        # ...
        if abs(kahan_sum(pix, d) - 1.) > atol:
            raise ValueError("probabilities do not sum to 1")

An error is only raised if the sum differs from 1.0 by more than some “machine epsilon.”

We also see that NumPy uses the Kahan summation to add the probabilities, an algorithm which significantly reduces the numerical error. I don't know if NumPy publicly exposes the Kahan summation (I could not find it). But it is not difficult to implement, here is a version from Rosetta Code:

def kahansum(input):
    summ = c = 0
    for num in input:
        y = num - c
        t = summ + y
        c = (t - summ) - y
        summ = t
    return summ

which would then be used as

totalScore = kahansum(item['strength'] for item in data)

Example:

>>> probs = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
>>> sum(probs)
0.9999999999999999
>>> kahansum(probs)
1.0

The computation of the probResult array can be simplified. First note that

value = (100.0/totalScore)*item['strength']
value2 = (1.0/100)*value

is equivalent to

value2 = item['strength'] / totalScore

and that would work for the case item['strength'] == 0 as well:

for item in data:
    value2 = item['strength'] / totalScore
    probResult.append(value2)

And now you can replace the loop with a list comprehension:

def makeProbabilities(data) -> list:
    totalScore = kahansum(item['strength'] for item in data)
    return [ item['strength'] / totalScore for item in data ]

Finally note that according to the Python naming conventions, function and variable names should be “snake case”: make_probabilities, total_score.

Answer 2

When using numpy, you should try to avoid for loops and instead favor vectorized operations.

This means that once you have extracted the desired data from your external data-structure, you store it in an np.ndarray and work forward from there. In your specific case, this is as simple as dividing the entire array with its sum; which is easily written:

def make_probabilities(data) -> np.ndarray:
    strengths = np.array([item['strength'] for item in data])
    return strengths / strengths.sum()

Note that:

• the division will be performed element-wise (this is know as broadcasting in numpy);
• division by 0 will return an array of np.nan and np.random.choice will then only return the first element; you may want to account for that depending on your use-case.