# TVDTemperingLift.py

Revision as of 21:56, 16 March 2022; view current revision

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()
```