TVDTemperingLift.py

From Werner KRAUTH

This is the program TVDTemperingLift.py, which allows one to compute the TVD as a function of time for the simulated tempering problem with the flat and the V-shaped probability distribution.

#
# Simulated tempering for the V-shaped stationary distribution --- Lifting
# algorithm
#
import random
import pylab
import numpy as np
import scipy.linalg as la
ReplicaChange = 0.1
n = 100
const = 4.0 / n ** 2
PiStat = {}
Table = []

for x in range(1, n + 1):
    Table.append((x, -1, 0))
    Table.append((x,  1, 0))
    Table.append((x, -1, 1))
    Table.append((x,  1, 1))
#
#   factor of 1/4 because the total must be normalized, with two liftings and
#   two replicas  
#
    PiStat[(x, -1,  0)] = 1.0 / float(n) / 4.0
    PiStat[(x,  1,  0)] = 1.0 / float(n) / 4.0
    PiStat[(x, -1, 1)] = const * abs( (n + 1) / 2 - x) / 4.0
    PiStat[(x,  1, 1)] = const * abs( (n + 1) / 2 - x) / 4.0
PiStat[(0, -1, 0)] = 0.0
PiStat[(0,  1, 0)] = 0.0
PiStat[(0, -1, 1)] = 0.0
PiStat[(0,  1, 1)] = 0.0
PiStat[(n + 1, -1,  0)] = 0.0
PiStat[(n + 1,  1,  0)] = 0.0
PiStat[(n + 1, -1, 1)] = 0.0
PiStat[(n + 1,  1, 1)] = 0.0
Position = (1,  1, 0)
print(Table)
PTrans   = np.zeros((4 * n, 4 * n))
for x in range(1, n + 1):
    for Dir in [-1, 1]:
        for Rep in [0, 1]:
            i = Table.index((x, Dir, Rep))
            k = Table.index((x, -Dir, Rep))
            if PiStat[(x + Dir, Dir, Rep)] > 0.0:
                j = Table.index((x + Dir, Dir, Rep))
                PTrans[i, j] = min(1.0, PiStat[(x + Dir, Dir,  Rep)] / PiStat[(x, Dir, Rep)])
                PTrans[i, k] = 1.0 - PTrans[i, j]
            else:
                PTrans[i, k] = 1.0
print('PTrans')
print(PTrans)

PReplica = np.zeros((4 * n, 4 * n))
for x in range(1, n + 1):
    for Dir in [-1, 1]:
        i = Table.index((x, Dir, 0))
        j = Table.index((x, Dir, 1))
        PReplica[i, j] = ReplicaChange * min(1.0, PiStat[(x, Dir, 1)] / PiStat[(x, Dir, 0)])
        PReplica[i, i] = 1.0 - PReplica[i, j]
        PReplica[j, i] = ReplicaChange * min(1.0, PiStat[(x, Dir,  0)] / PiStat[(x, Dir, 1)])
        PReplica[j, j] = 1.0 - PReplica[j, i]
print('PReplica')
print(PReplica)
P = PTrans @ PReplica
eigvals, eigvecsl, eigvecsr = la.eig(P, left=True)
eigvals = eigvals.real
print('Eigenvalues')
print(eigvals)
gap = 1000.0
for k in range(4 * n):
    if abs(eigvals[k] - 1.0) < 1.e-5: max_val = k
    else:
        gap = min(gap, abs(eigvals[k] - 1.0))
print('gap', gap, ReplicaChange)
print('maxval', max_val)
PiStat = eigvecsl[:, max_val]
PiStat = PiStat.real
PiStatSum = PiStat.sum()
PiStat0 = []
PiStat1 = []
for k in range(4 * n):
    if k % 4 == 0:
        PiStat0.append(4.0 * PiStat[k] / PiStatSum)
    elif k % 4 == 2:
        PiStat1.append(4.0 * PiStat[k] / PiStatSum)
XVec = [k for k in range(1, n + 1)]
print(PiStat0)
pylab.plot(XVec, PiStat0)
pylab.plot(XVec, PiStat1)
pylab.xlabel("$x$ (position) ")
pylab.title("Simulated tempering (Lifted MCMC), Left EV transition matrix")
pylab.show()