TVDTemperingRev.py

From Werner KRAUTH

---

Revision as of 10:23, 8 September 2022

2. Simulated tempering for the V-shaped stationary distribution --- Metropolis
3. algorithm

import random import pylab import numpy as np for n in [10, 20, 40, 80, 160, 320]:

   ReplicaChange = 0.1
   const = 4.0 / n ** 2
   PiStat = {}
   Table = []

   for x in range(1, n + 1):
       Table.append((x, 0))
       Table.append((x, 1))

2. factor of 1/2 because the total must be normalized

       PiStat[(x, 0)] = 1.0 / float(n) / 2.0
       PiStat[(x, 1)] = const * abs( (n + 1) / 2 - x) / 2.0
   PiStat[(0, 0)] = 0.0
   PiStat[(0, 1)] = 0.0
   PiStat[(n + 1, 0)] = 0.0
   PiStat[(n + 1, 1)] = 0.0
   PTrans   = np.eye(2 * n)
   Pi = np.zeros([2 * n])
   for x in range(1, n + 1):
       for Rep in [0, 1]:
           i = Table.index((x, Rep))
           Pi[i] = PiStat[(x, Rep)]
           for Dir in [-1, 1]:
               if PiStat[(x + Dir, Rep)] > 0.0:
                   j = Table.index((x + Dir, Rep))
                   PTrans[i, j] = min(1.0, PiStat[(x + Dir, Rep)] / PiStat[(x, Rep)]) / 2.0
                   PTrans[i, i] -= PTrans[i, j]

   PReplica = np.zeros((2 * n,2 * n))
   for x in range(1, n + 1):
       i = Table.index((x,0))
       j = Table.index((x,1))
       PReplica[i, j] = ReplicaChange * min(1.0, PiStat[(x, 1)] / PiStat[(x, 0)])
       PReplica[i, i] = 1.0 - PReplica[i, j]
       PReplica[j, i] = ReplicaChange * min(1.0, PiStat[(x, 0)] / PiStat[(x, 1)])
       PReplica[j, j] = 1.0 - PReplica[j, i]
   P = PTrans @ PReplica
   Pit = np.zeros([2 * n])
   Pit[0] = 1.0
   xvalues = []
   yvalues = []
   iter = 0
   while True:
       iter += 1
       Pit = Pit @ P
       TVD = sum(np.absolute(Pi - Pit) / 2.0)
       xvalues.append(iter / float(n ** 2))
       yvalues.append(TVD)
       if TVD < 0.1: break
   pylab.plot(xvalues,yvalues, label='$n =$ '+str(n))

pylab.legend(loc='upper right') pylab.xlabel("$t/ n^2$ (rescaled time) ") pylab.ylabel("TVD") pylab.title("TVD rev tempering on the path graph of $n$ sites") pylab.show()