Walker.py

From Werner KRAUTH

Revision as of 15:07, 6 June 2024 Werner (Talk | contribs) ← Previous diff		Current revision Werner (Talk | contribs) (→Python program)
Line 34:		Line 34:
	pi = [random.uniform(0.1, 10.0) for i in range(N)]		pi = [random.uniform(0.1, 10.0) for i in range(N)]
	pisum = sum(a for a in pi)		pisum = sum(a for a in pi)
-	pi = [[pi[k] / pisum, k] for k in range(N)]	+	pi = [[pi[k] / pisum, k] for k in range(N)] # see Comment below
	N_walker, walker_mean, walker_table = walker_setup (pi)		N_walker, walker_mean, walker_table = walker_setup (pi)
Line 50:		Line 50:
	for i in range(N):		for i in range(N):
	print(pi[i], data[i])		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.

---

Current revision

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.