Walker.py

From Werner KRAUTH

Context

This page is part of my 2024 Beg Rohu Lectures on "The second Markov chain revolution" at the Summer School "Concepts and Methods of Statistical Physics" (3 - 15 June 2024).

Python program

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)]   # see Comment below

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])

Comment: In the line above that contains the comment, the normalized probabilities are entered into the list-of-lists pi. Each element of the list pi is thus a list of two elements, the probability pi(i) and the index i. At the beginning, this is a bit pointless. Later on in the process, some cutting-up and reaarranging takes place, and it becomes convenient to know the 'alias' of each stick, that is the original element i that each stick originally came from.