# Parallel Tempering

## General Chain-Structured Models

There is an important probability distribution used in applications (Liu, 2008):

$$
\pi(\mathbf x) \propto \exp\Big{-\sum\_{i=1}^dh\_i(x\_{i-1},x\_i)\Big} \equiv \exp(-H(\mathbf x))
$$

where $$\mathbf x = (x\_0,x\_1,\ldots,x\_d)$$. This type of model exhibits a so-called "Markovian structure" because

$$
\pi(x\_i\mid \mathbf x\_{\[-i]}) \propto \exp{-h(x\_{i-1},x\_i)-h(x\_i,x\_{i+1})},.
$$

We can simulate from $$\pi(\mathbf x)$$ by an exact method. Note that

$$
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],.
$$

The simulating procedure is as follows:

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LUZSTjNo_5AXQwol1GC%2F-LUZSUR9iRbYqcNtEF2r%2Fexact.png?generation=1545726196447794\&alt=media)

## Exact Method for Ising Model

For a one-dimensional Ising Model,

$$
\pi(\mathbf x)=Z^{-1}\exp{\beta(x\_0x\_1+\cdots+x\_{d-1}x\_d)},,
$$

where $$x\_i\in{-1,1}$$. And thus the simulating procedure is much easy. Note that for $$t=1,\ldots,d$$,

$$
V\_t(x) = (e^{-\beta}+e^{\beta})^t
$$

and $$Z=2(e^{-\beta}+e^{\beta})^d$$.

We can use the following Julia program to simulate from $$\pi(\mathbf x)$$:

```julia
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:

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LUZRIDv_JiJ40YqHscJ%2F-LUZRIlC-0oFSD2tadFH%2Fpt.png?generation=1545725886591462\&alt=media)

I wrote the following Julia program to implement this procedure and reproduce the simulation example of Liu (2008):

```julia
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](https://github.com/szcf-weiya/MonteCarlo/blob/master/Tempering/Ising-Model/Ising-model.jl) for complete source code.

Finally, I reproduce the Figure 10.1 as follows:

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LUZRIDv_JiJ40YqHscJ%2F-LUZRIlEarA0Rue7Ja5h%2FIsingHist-rep9.png?generation=1545725886567951\&alt=media)

and the corresponding autocorrelation plots:

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LUZRIDv_JiJ40YqHscJ%2F-LUZRIlG2zTatDHRhbFT%2FIsingACF-rep9.png?generation=1545725887131585\&alt=media)

## References

[Liu, J. S. (2008). Monte Carlo strategies in scientific computing. Springer Science & Business Media.](https://github.com/szcf-weiya/MonteCarlo/tree/ea95daee44829f9736e1cef968ba1095b6c21c60/References/Monte-Carlo-Strategies-in-Scientific-Computing.pdf)