The pmf of Ising model is considered as $p(\boldsymbol{x})=\frac{1}{Z(\theta)} exp\left\{ \underset{\left(s,t\right)\in E}{\sum\theta_{st}}x_{s}x_{t}\right\},\quad \boldsymbol{x}\in \{-1,1\}^n$, where $E$ is the set that contains the edges of the graph and $Z(\theta)=\underset{x}\sum{ exp\left\{ \underset{\left(s,t\right)\in E}{\sum\theta_{st}}x_{s}x_{t}\right\}}$. I have read that the joint of two random variables $x_i, x_j$ is $p\left(x_i,x_j\right)=\left[1+x_{i}x_{j}E\left[X_{i}X_{j}\right]\right]/4$. in https://arxiv.org/pdf/1604.06749.pdf on page 8. How can we show that?
I know that we can express $E\left[X_{i}X_{j}\right]$ in terms of $\theta$: $E\left[X_{i}X_{j}\right]=\frac{\partial \ln Z(\theta)}{\partial\theta_{i,j}}$ but i can not use it to find the joint pmf.
I have posted the same question on cross-validate.