TVDMetroPath.py

From Werner KRAUTH

#
# TVD for the Metropolis algorithm on the path graph.
#
import random
import pylab
import numpy as np
model = 'Flat'
model = 'VShape'
for n in [10, 20, 40, 80, 160, 320]:
    const = 4.0 / n ** 2
    PiStat = {}
    Table = []
    for x in range(1, n + 1):
        Table.append(x)
#
#   flat stationary probability distribution.
#
        if model == 'Flat':
            PiStat[x] = 1.0 / float(n)
        elif model == 'VShape':
            PiStat[x] = const * abs( (n + 1) / 2 - x)
    PiStat[0] = 0.0
    PiStat[n + 1] = 0.0
    PTrans   = np.eye(n)
    Pi = np.zeros([n])
    for x in range(1, n + 1):
        i = Table.index(x)
        Pi[i] = PiStat[x]
        for Dir in [-1, 1]:
            if PiStat[x + Dir] > 0:
                j = Table.index(x + Dir)
                PTrans[i, j] = min(1.0, PiStat[x + Dir] / PiStat[x]) / 2.0
                PTrans[i, i] -= PTrans[i, j]
    Pit = np.zeros([n])
    Pit[0] = 1.0
    xvalues = []
    yvalues = []
    iter = 0
    while True:
        iter += 1
        Pit = np.array(Pit)
        Pit = Pit@PTrans
        TVD = sum(np.absolute(Pi - Pit) / 2.0)
        xvalues.append(iter / float(n ** 2 * np.log(n)))
        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 for the Metropolis algorithm on the path graph of $n$ sites")
pylab.show()