Flawed random number generator?

Question

I used this weighted random number generator.

import random

def weighted_choice(weights):
    totals = []
    running_total = 0

    for w in weights:
        running_total += w
        totals.append(running_total)

    rnd = random.random() * running_total
    for i, total in enumerate(totals):
        if rnd < total:
            return i

as follows:

# The meaning of this dict is a little confusing, so here's the explanation:
# The keys are numbers and values are weights of its occurence and values - 1
# are weights of its disoccurence. You can imagine it like biased coins
# (except for 2 which is fair coin).
probabilities = { 0 : 1.0, 1 : 1.0, 2 : 0.5, 3 : 0.45, 4 : 0.4, 5 : 0.35,
                    6 : 0.3, 7 : 0.25, 8 : 0.2, 9 : 0.15, 10 : 0.1
                  }
  numberOfDeactivations = []
  for number in probabilities.keys():
    x = weighted_choice([probabilities[number], 1 - probabilities[number]])
    if x == 0:
      numberOfDeactivations.append(number)
  print "chance for ", repr(numberOfDeactivations)

I see quite often 7, 8, 9, 10 in the result.

Is there some proof or guarantee that this is correct to probability theory?

Answer 1

This is mathematically correct. It's an application of inverse transform sampling (although the reason it works in this case should be relatively intuitive).

I don't know Python, so I can't say whether there are any subtleties that make this particualr implementation invalid.

Answer 2

Edit: as a side note: I think your code is equivalent to

import random
probabilities = { 0 : 1.0, 1 : 1.0, 2 : 0.5, 3 : 0.45, 4 : 0.4, 5 : 0.35,
                    6 : 0.3, 7 : 0.25, 8 : 0.2, 9 : 0.15, 10 : 0.1}
numberOfDeactivations=filter(
         lambda kv:random.random()<=probabilities[kv] , probabilities)

Original answer:

The method is correct. Below is a complete example, creating the frequency table and comparing it with the requested weights.

With 100000 iterations there's nothing indicating that you don't get what you requested. The 'expected' is the probability you requested, 'got' is the fraction of times you actually got that value. Ratio should be close to 1 and it is:

  0, expected: 0.2128 got: 0.2107 ratio: 1.0100
  1, expected: 0.2128 got: 0.2145 ratio: 0.9921
  2, expected: 0.1064 got: 0.1083 ratio: 0.9825
  3, expected: 0.0957 got: 0.0949 ratio: 1.0091
  4, expected: 0.0851 got: 0.0860 ratio: 0.9900
  5, expected: 0.0745 got: 0.0753 ratio: 0.9884
  6, expected: 0.0638 got: 0.0635 ratio: 1.0050
  7, expected: 0.0532 got: 0.0518 ratio: 1.0262
  8, expected: 0.0426 got: 0.0418 ratio: 1.0179
  9, expected: 0.0319 got: 0.0323 ratio: 0.9881
 10, expected: 0.0213 got: 0.0209 ratio: 1.0162

 A total of 469633 numbers where generated for this table. 

Here's the code:

import random

def weighted_choice(weights):
    totals = []
    running_total = 0
    for w in weights:
        running_total += w
        totals.append(running_total)
    rnd = random.random() * running_total
    for i, total in enumerate(totals):
        if rnd < total:
            return i


counts={ k:0 for k in range(11)}
probabilities = { 0 : 1.0, 1 : 1.0, 2 : 0.5, 3 : 0.45, 4 : 0.4, 5 : 0.35,
                    6 : 0.3, 7 : 0.25, 8 : 0.2, 9 : 0.15, 10 : 0.1
                  }

for x in range(100000):
  numberOfDeactivations = []
  for number in probabilities.keys():
    x = weighted_choice([probabilities[number], 1 - probabilities[number]])
    if x == 0:
      numberOfDeactivations.append(number)
  for k in numberOfDeactivations:
    counts[k]+=1.0

sums=sum(counts.values())
counts2=[x*1.0/sums for x in counts.values()]

print "ratio expected frequency to requested:":

# make the probabilities real probabilities instead of weights:
psum=sum(probabilities.values())
for k in probabilities:
    probabilities[k]=probabilities[k]/psum

for k in probabilities:
    print "%3d, expected: %6.4f got: %6.4f ratio: %6.4f" %(k,probabilities[k],counts2[k], probabilities[k]/counts2[k])