Parallel Tempering

General Chain-Structured Models
There is an important probability distribution used in applications (Liu, 2008):
π(x)∝exp⁡{−∑i=1dhi(xi−1,xi)}≡exp⁡(−H(x))\pi(\mathbf x) \propto \exp\Big\{-\sum_{i=1}^dh_i(x_{i-1},x_i)\Big\} \equiv \exp(-H(\mathbf x))π(x)∝exp{−i=1∑d​hi​(xi−1​,xi​)}≡exp(−H(x))
where x=(x0,x1,…,xd)\mathbf x = (x_0,x_1,\ldots,x_d)x=(x0​,x1​,…,xd​)
. This type of model exhibits a so-called "Markovian structure" because
π(xi∣x[−i])∝exp⁡{−h(xi−1,xi)−h(xi,xi+1)} .\pi(x_i\mid \mathbf x_{[-i]}) \propto \exp\{-h(x_{i-1},x_i)-h(x_i,x_{i+1})\}\,.π(xi​∣x[−i]​)∝exp{−h(xi−1​,xi​)−h(xi​,xi+1​)}.
We can simulate from π(x)\pi(\mathbf x)π(x)
 by an exact method. Note that
Z≡∑xexp⁡{−H(x)}=∑xd[⋯[∑x1{∑x0e−h(x0,x1)}e−h2(x1,x2)]⋯] .Z\equiv \sum_{\mathbf x}\exp\{-H(\mathbf x)\} = \sum_{x_d}\Big[\cdots\Big[ \sum_{x_1}\Big\{\sum_{x_0}e^{-h(x_0,x_1)}\Big\}e^{-h_2(x_1,x_2)} \Big] \cdots\Big]\,.Z≡x∑​exp{−H(x)}=xd​∑​[⋯[x1​∑​{x0​∑​e−h(x0​,x1​)}e−h2​(x1​,x2​)]⋯].
The simulating procedure is as follows:
Exact Method for Ising Model
For a one-dimensional Ising Model,
π(x)=Z−1exp⁡{β(x0x1+⋯+xd−1xd)} ,\pi(\mathbf x)=Z^{-1}\exp\{\beta(x_0x_1+\cdots+x_{d-1}x_d)\}\,,π(x)=Z−1exp{β(x0​x1​+⋯+xd−1​xd​)},
where xi∈{−1,1}x_i\in\{-1,1\}xi​∈{−1,1}
. And thus the simulating procedure is much easy. Note that for t=1,…,dt=1,\ldots,dt=1,…,d
,
Vt(x)=(e−β+eβ)tV_t(x) = (e^{-\beta}+e^{\beta})^tVt​(x)=(e−β+eβ)t
and Z=2(e−β+eβ)dZ=2(e^{-\beta}+e^{\beta})^dZ=2(e−β+eβ)d
.
We can use the following Julia program to simulate from π(x)\pi(\mathbf x)π(x)
:
function exactMethod(M = 2000; beta = 10, d = 100)
    V = ones(d)
    V[1] = exp(beta) + exp(-beta)
    for t = 2:d
        V[t] = V[1] * V[t-1]
    end
    Z = 2*V[d]
    x = ones(Int, d)
    x[d] = 2*rand(Bernoulli(1/2))-1
    for t = (d-1):1
#        p1 = V[t] * exp(x[d] * 1)
#        p2 = V[t] * exp(x[d] * (-1))
        p1 = exp(x[d])
        p2 = exp(-x[d])
        x[t] = 2*rand(Bernoulli(p1/(p1+p2)))-1
    end
    return x
end
Parallel Tempering for Ising Model
Liu (2008) also introduced the parallel tempering strategy:
I wrote the following Julia program to implement this procedure and reproduce the simulation example of Liu (2008):
function parallelTemper(M = 2000, T = [0.1, 0.2, 0.3, 0.4]; d = 100, alpha0 = T[1])
    I = length(T)
    # initial 
    x = 2*rand(Bernoulli(1/2), d, I) .- 1
    num1 = zeros(3)
    num2 = zeros(3)
    res = zeros(M, 4)
    for m = 1:M
        if rand() < alpha0 # parallel step
            #for i = 1:I
                # !!! Doesn't work !!!
                # y1 = x[:,i]
                # sampleX!(x[:,i], T[i])
                # y2 = x[:,i]
                # println(sum(abs.(y2.-y1)))
            #end
            sampleX!(x, T)
        else
#            idx = sample(1:I, 2, replace = false) # not neigbor
            idx1 = sample(1:(I-1))
            num1[idx1] += 1
            idx = [idx1, idx1+1]
#            rho = pdfIsing(x[:,idx[2]], T[idx[1]]) * pdfIsing(x[:,idx[1]], T[idx[2]]) / (pdfIsing(x[:,idx[1]], T[idx[1]]) * pdfIsing(x[:,idx[2]], T[idx[2]]))
            rho = logpdfIsing(x[:,idx[2]], T[idx[1]]) + logpdfIsing(x[:,idx[1]], T[idx[2]]) - (logpdfIsing(x[:,idx[1]], T[idx[1]]) + logpdfIsing(x[:,idx[2]], T[idx[2]]))
            if log(rand()) < rho
                # swappin step
                num2[idx1] += 1
                tmp = copy(x[:, idx[1]])
                x[:,idx[1]] .= x[:, idx[2]]
                x[:,idx[2]] .= tmp
            end
        end
        res[m, :] = sum(x, dims=1)
    end
    return res, num2 ./ num1
end
It is important to note that we should use logpdf to avoid too big values, in which the rejection part doesn't work well. Refer to Ising-model.jl for complete source code.
Finally, I reproduce the Figure 10.1 as follows:
and the corresponding autocorrelation plots:
References
​Liu, J. S. (2008). Monte Carlo strategies in scientific computing. Springer Science & Business Media.​
Tempering
Misc
Last modified 4yr ago