TVDTemperingRev.py

From Werner KRAUTH

(Difference between revisions)
Jump to: navigation, search
Revision as of 10:23, 8 September 2022
Werner (Talk | contribs)

← Previous diff		 Current revision
Werner (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()
Retrieved from "http://www.lps.ens.fr/~krauth/index.php/TVDTemperingRev.py"