# 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:

![](/files/-LUZSUR9iRbYqcNtEF2r)

## 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:

![](/files/-LUZRIlC-0oFSD2tadFH)

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:

![](/files/-LUZRIlEarA0Rue7Ja5h)

and the corresponding autocorrelation plots:

![](/files/-LUZRIlG2zTatDHRhbFT)

## 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)


---

# Agent Instructions: Querying This Documentation

If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question.

Perform an HTTP GET request on the current page URL with the `ask` query parameter:

```
GET https://mc.hohoweiya.xyz/tempering/ising-model.md?ask=<question>
```

The question should be specific, self-contained, and written in natural language.
The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.

Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.