Python code to sample from a probability distribution described by the logarithms of the probabilities (log probabilities). Includes code for log sum (log_add)

log_add.py

def log_add(x, y):
    maximum = max(x,y)
    minimum = min(x,y)
    if(abs(maximum - minimum) > 30):
        # the difference is too small, return the just the maximum
        return maximum
    return maximum + log(1 + pow(2, minimum - maximum) ,2)

sample_from_log_prob.py

import numpy as np
import random
from math import log
import log_add

# Takes as input a log probability vector. The base of the logarithm is 2.
# Returns a sampled position according to the corresponding log probabilities
# Uses the formula from
# http://blog.smola.org/post/987977550/log-probabilities-semirings-and-floating-point-numbers
def sample_from_log_prob(A):
    max_log_prob = max(A)
    C = [A[0]]
    for a in A[1:]:
        C.append(log_add(C[-1], a))
    C_pos = [ -c for c in reversed(C)]
    r = log(random.random(),2)
    pos = np.searchsorted(C_pos,-r)
    return len(C) - pos