# TVDMetroLift.py

### From Werner KRAUTH

(Difference between revisions)

Revision as of 10:10, 8 September 2022Werner (Talk | contribs) ← Previous diff |
Current revisionWerner (Talk | contribs) |
||

Line 44: |
Line 44: | ||

PResampling[i, j] = 1.0 / n | PResampling[i, j] = 1.0 / n | ||

PResampling[i, i] = 1.0 - PResampling[i,j] | PResampling[i, i] = 1.0 - PResampling[i,j] | ||

- | + | ||

P = PTrans @ PResampling | P = PTrans @ PResampling | ||

Pit = np.zeros([2 * n]) | Pit = np.zeros([2 * n]) |

## Current revision

# # TVD for the Lifted 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, -1)) Table.append((x, 1)) if model == 'Flat': PiStat[(x, -1)] = 1.0 / float(n) / 2.0 PiStat[(x, +1)] = 1.0 / float(n) / 2.0 elif model == 'VShape': PiStat[(x, -1)] = const * abs( (n + 1) / 2 - x) / 2.0 PiStat[(x, 1)] = const * abs( (n + 1) / 2 - x) / 2.0 PiStat[(0, -1)] = 0.0 PiStat[(0, 1)] = 0.0 PiStat[(n + 1, -1)] = 0.0 PiStat[(n + 1, 1)] = 0.0 PTrans = np.zeros((2 * n, 2 * n)) Pi = np.zeros([2 * n]) for x in range(1, n + 1): for Sigma in [-1, 1]: i = Table.index((x, Sigma)) Pi[i] = PiStat[(x, Sigma)] k = Table.index((x, -Sigma)) if PiStat[(x + Sigma, Sigma)] > 0.0: j = Table.index((x + Sigma, Sigma)) PTrans[i, j] = min(1.0, PiStat[(x + Sigma, Sigma)] / PiStat[(x, Sigma)]) PTrans[i, k] = 1.0 - PTrans[i, j] else: PTrans[i, k] = 1.0 PResampling = np.zeros((2 * n, 2 * n)) for x in range(1, n + 1): for Sigma in [-1, 1]: i = Table.index((x, Sigma)) j = Table.index((x, -Sigma)) PResampling[i, j] = 1.0 / n PResampling[i, i] = 1.0 - PResampling[i,j] P = PTrans @ PResampling Pit = np.zeros([2 * n]) Pit[0] = 1.0 xvalues = [] yvalues = [] iter = 0 while True: iter += 1 Pit = np.array(Pit) 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$ (rescaled time) ") pylab.ylabel("TVD") pylab.title("TVD for the lifted Metropolis algorithm on the path graph of $n$ sites") pylab.show()

~