# TVDTemperingRev.py

(Difference between revisions)
Jump to: navigation, search
 Revision as of 10:23, 8 September 2022Werner (Talk | contribs)← Previous diff Current revisionWerner (Talk | contribs) Line 1: Line 1: - # + # - # Simulated tempering for the V-shaped stationary distribution --- Metropolis + # Simulated tempering for the V-shaped stationary distribution --- Metropolis - # algorithm + # algorithm - # + # - import random + import random - import pylab + import pylab - import numpy as np + import numpy as np - for n in [10, 20, 40, 80, 160, 320]: + for n in [10, 20, 40, 80, 160, 320]: - ReplicaChange = 0.1 + ReplicaChange = 0.1 - const = 4.0 / n ** 2 + const = 4.0 / n ** 2 - PiStat = {} + PiStat = {} - Table = [] + Table = [] - for x in range(1, n + 1): + for x in range(1, n + 1): - Table.append((x, 0)) + Table.append((x, 0)) - Table.append((x, 1)) + Table.append((x, 1)) - # + # - # factor of 1/2 because the total must be normalized + # factor of 1/2 because the total must be normalized - # + # - PiStat[(x, 0)] = 1.0 / float(n) / 2.0 + PiStat[(x, 0)] = 1.0 / float(n) / 2.0 - PiStat[(x, 1)] = const * abs( (n + 1) / 2 - x) / 2.0 + PiStat[(x, 1)] = const * abs( (n + 1) / 2 - x) / 2.0 - PiStat[(0, 0)] = 0.0 + PiStat[(0, 0)] = 0.0 - PiStat[(0, 1)] = 0.0 + PiStat[(0, 1)] = 0.0 - PiStat[(n + 1, 0)] = 0.0 + PiStat[(n + 1, 0)] = 0.0 - PiStat[(n + 1, 1)] = 0.0 + PiStat[(n + 1, 1)] = 0.0 - PTrans = np.eye(2 * n) + PTrans = np.eye(2 * n) - Pi = np.zeros([2 * n]) + Pi = np.zeros([2 * n]) - for x in range(1, n + 1): + for x in range(1, n + 1): - for Rep in [0, 1]: + for Rep in [0, 1]: - i = Table.index((x, Rep)) + i = Table.index((x, Rep)) - Pi[i] = PiStat[(x, Rep)] + Pi[i] = PiStat[(x, Rep)] - for Dir in [-1, 1]: + for Dir in [-1, 1]: - if PiStat[(x + Dir, Rep)] > 0.0: + if PiStat[(x + Dir, Rep)] > 0.0: - j = Table.index((x + Dir, Rep)) + j = Table.index((x + Dir, Rep)) - PTrans[i, j] = min(1.0, PiStat[(x + Dir, Rep)] / PiStat[(x, Rep)]) / 2.0 + PTrans[i, j] = min(1.0, PiStat[(x + Dir, Rep)] / PiStat[(x, Rep)]) / 2.0 - PTrans[i, i] -= PTrans[i, j] + PTrans[i, i] -= PTrans[i, j] - PReplica = np.zeros((2 * n,2 * n)) + PReplica = np.zeros((2 * n,2 * n)) - for x in range(1, n + 1): + for x in range(1, n + 1): - i = Table.index((x,0)) + i = Table.index((x,0)) - j = Table.index((x,1)) + j = Table.index((x,1)) - PReplica[i, j] = ReplicaChange * min(1.0, PiStat[(x, 1)] / PiStat[(x, 0)]) + PReplica[i, j] = ReplicaChange * min(1.0, PiStat[(x, 1)] / PiStat[(x, 0)]) - PReplica[i, i] = 1.0 - PReplica[i, j] + PReplica[i, i] = 1.0 - PReplica[i, j] - PReplica[j, i] = ReplicaChange * min(1.0, PiStat[(x, 0)] / PiStat[(x, 1)]) + PReplica[j, i] = ReplicaChange * min(1.0, PiStat[(x, 0)] / PiStat[(x, 1)]) - PReplica[j, j] = 1.0 - PReplica[j, i] + PReplica[j, j] = 1.0 - PReplica[j, i] - P = PTrans @ PReplica + P = PTrans @ PReplica - Pit = np.zeros([2 * n]) + Pit = np.zeros([2 * n]) - Pit[0] = 1.0 + Pit[0] = 1.0 - xvalues = [] + xvalues = [] - yvalues = [] + yvalues = [] - iter = 0 + iter = 0 - while True: + while True: - iter += 1 + iter += 1 - Pit = Pit @ P + Pit = Pit @ P - TVD = sum(np.absolute(Pi - Pit) / 2.0) + TVD = sum(np.absolute(Pi - Pit) / 2.0) - xvalues.append(iter / float(n ** 2)) + xvalues.append(iter / float(n ** 2)) - yvalues.append(TVD) + yvalues.append(TVD) - if TVD < 0.1: break + if TVD < 0.1: break - pylab.plot(xvalues,yvalues, label='\$n =\$ '+str(n)) + pylab.plot(xvalues,yvalues, label='\$n =\$ '+str(n)) - pylab.legend(loc='upper right') + pylab.legend(loc='upper right') - pylab.xlabel("\$t/ n^2\$ (rescaled time) ") + pylab.xlabel("\$t/ n^2\$ (rescaled time) ") - pylab.ylabel("TVD") + pylab.ylabel("TVD") - pylab.title("TVD rev tempering on the path graph of \$n\$ sites") + pylab.title("TVD rev tempering on the path graph of \$n\$ sites") - pylab.show() + pylab.show()

## Current revision

```#
# Simulated tempering for the V-shaped stationary distribution --- Metropolis
# 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))
#
#   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()
```