Walker.py

From Werner KRAUTH

import random
def walker_setup (pi):
    """compute the lookup table for Walker's algorithm"""
    N_walker = len(pi)
    walker_mean = sum(a[0] for a in pi) / float(N_walker)
    long_s = []
    short_s = []
    for p in pi:
        if p[0] > walker_mean:
            long_s.append(p[:])
        else:
            short_s.append(p[:])
    walker_table = []
    for k in range(N_walker - 1):
        e_plus = long_s.pop()
        e_minus = short_s.pop()
        walker_table.append((e_minus[0], e_minus[1], e_plus[1]))
        e_plus[0] = e_plus[0] - (walker_mean - e_minus[0])
        if e_plus[0] < walker_mean:
            short_s.append(e_plus)
        else:
            long_s.append(e_plus)
    if long_s != []:
        walker_table.append((long_s[0][0], long_s[0][1], long_s[0][1]))
    else:
        walker_table.append((short_s[0][0], short_s[0][1], short_s[0][1]))
    return N_walker, walker_mean, walker_table

N = 12
pi = [random.uniform(0.1, 10.0) for i in range(N)]
pisum = sum(a for a in pi)
pi = [[pi[k] / pisum, k] for k in range(N)]

N_walker, walker_mean, walker_table = walker_setup (pi)

data = [0] * N
NSample = 100000 * N
for iter in range(NSample):
    i = random.randint (0, N - 1)
    Upsilon = random.uniform (0.0, walker_mean)
    if Upsilon < walker_table[i][0]:
        sample = walker_table[i][1]
    else:
        sample = walker_table[i][2]
    data[sample] += 1.0 / NSample
for i in range(N):
    print(pi[i], data[i])