Hammersley-Clifford theorem The paper spatial interaction and the statistical analysis of lattice systems by Besag (1974) presents an alternative proof for the Hammersley-Clifford theorem.
In order to prove the HM theorem, Besag says that that under the following assumptions

*

*The number of possible values at each site/vertex is finite.

*$0$ is a possible value at every site.

the following expansion holds for every cdf $\mathbb{P}(\mathbf{x})$ and is unique on the sample space $\Omega$,
$$
\begin{align}
Q(\textbf{x}) &=  \mathbb{P}(\mathbf{x})/\mathbb{P}(\mathbf{0})
\\  
&= \sum_{1 \leq i \leq n} x_iG_i(x_i)
 +\sum_{1 \leq i \leq j \leq n}x_ix_jG_{i,j}(x_i,x_j)\\ &+ \sum_{1 \leq i \leq j \leq k \leq n} x_ix_jx_kG_{i,j,k}(x_i,x_j,x_k)+\ldots + x_1x_2 \ldots x_nG_{1,2,\ldots,n}(x_1,x_2,\ldots,x_n)
\end{align}
$$
This equality looks slightly related to the inclusion-exclusion principle and in fact Besag says that we can use it to prove a HM for some restricted classes of lattices. How did Besag arrived at this equality? I'm really clueless.
My question is related to this  question that remains unanswered.
 A: $\newcommand{\x}{\mathbf x}\newcommand{\tx}{\tilde{\mathbf x}}\newcommand{\0}{\mathbf0}\newcommand{\R}{\mathbb R}$You are citing Besag's paper (formulas (3.1) and (3.3) there) incorrectly. What the paper actually has is this:
\begin{equation*}
    \begin{aligned}
Q(\x)&:=\ln\frac{P(\x)}{P(\0)}
\\  
&= \sum_{1\le i\le n} x_iG_i(x_i)
 +\sum_{1\le i<j\le n}x_ix_jG_{i,j}(x_i,x_j)\\ 
 &+ \sum_{1\le i<j<k\le n} x_ix_jx_kG_{i,j,k}(x_i,x_j,x_k) +\cdots \\ 
 &+x_1\cdots x_nG_{1,\dots,n}(x_1,\ldots,x_n). 
\end{aligned}
\tag{1}\label{1}
\end{equation*}
This is indeed a version of the inclusion-exclusion formula. To show this, identify the vector $\x=(x_1,\dots,x_n)\in\R^n$ with the function $[n]\ni i\mapsto x_i$, which will also be denoted by $\x$; here, as usual, $[n]:=\{1,\dots,n\}$. For each $J\subseteq[n]$, let $\tx_J$ be the vector/function on $[n]$ such that $(\tx_J)(i):=x_i$ for $i\in J$ and $(\tx_J)(i):=0$ for $i\notin J$. Let $\x_J:=\x|_J=(x_i\colon i\in J)$, the restriction of the function $\x$ to $J$, so that $\tx_J$ is the extension by zeroes of the the function $\x_J$ from $J$ to $[n]$. Let $|J|$ denote the cardinality of $J$.
For any $K\subseteq[n]$, let
\begin{equation*}
    H_K(\x_K):=\sum_{J\subseteq K}(-1)^{|K|-|J|}Q(\tx_J);   
\end{equation*}
note that the last expression depends on $\x$ only through $\x_K$, so that $H_K(\x_K)$ is well defined.
Then
\begin{equation*}
    \begin{aligned}
    &\sum_{K\subseteq[n]}H_K(\x_K) \\ 
    &=\sum_{K\subseteq[n]}\sum_{J\subseteq K}(-1)^{|K|-|J|}Q(\tx_J) \\ 
    &=\sum_{J\subseteq[n]}Q(\tx_J)
    \sum_{K: J\subseteq K\subseteq[n]}(-1)^{|K|-|J|} \\ 
    &=\sum_{J\subseteq[n]}Q(\tx_J)
    \sum_{k=|J|}^n(-1)^{k-|J|}\binom{n-|J|}{k-|J|} \\ 
    &=\sum_{J\subseteq[n]}Q(\tx_J)
    \sum_{m=0}^{n-|J|}(-1)^m\binom{n-|J|}m \\ 
    &=\sum_{J\subseteq[n]}Q(\tx_J)
    \,1(|J|=n) \\ 
    &=Q(\tx_{[n]})=Q(\x).  
\end{aligned}
\tag{2}\label{2}
\end{equation*}
Note also that $H_K(\x_K)=0$ if $x_j=0$ for some $j\in K$; this follows in view of the natural bijection between the (set of all subsets $J$ of the -- necessarily nonempty -- set $K$ such that $J\ni j$) and (the set of all subsets $J$ of $K$ such that $J\not\ni j$).
So, $H_K(\x_K)=0$ if $p_K(\x)=0$, where
\begin{equation*}
    p_K(\x):=\prod_{j\in K}x_j. 
\end{equation*}
So, letting $G_K(\x_K):=H_K(\x_K)/p_K(\x)$ if $p_K(\x)\ne0$ and $G_K(\x_K):=0$ if $p_K(\x)=0$, we get $H_K(\x_K)=p_K(\x)G_K(\x_K)$ for all $\x$ and all $K$. Thus, by \eqref{2},
\begin{equation*}
    Q(\x)=\sum_{K\subseteq[n]}p_K(\x)G_K(\x_K), \tag{3}\label{3}
\end{equation*}
which is just another, more compact way of writing \eqref{1}.
Of course, one may refer to $G_K(\x_K)$ as a partial divided difference (of the function $Q$) of order $|K|$ with respect to $\x_K=(x_i\colon i\in K)$. In particular,
\begin{equation*}
    G_{\{1\}}(x_1)=\frac{Q(x_1,0,\dots,0)-Q(0,0,\dots,0)}{x_1}
    \Big(=\frac{Q(x_1,0,\dots,0)}{x_1}\Big)
\end{equation*}
if $x_1\ne0$ and
\begin{equation*}
    G_{\{1,2\}}(x_1,x_2)=\frac{Q(x_1,x_2,0,\dots,0)-Q(x_1,0,0,\dots,0)
    -Q(0,x_2,0,\dots,0)+Q(0,0,\dots,0)}{x_1x_2}  
\end{equation*}
$x_1x_2\ne0$.
So, \eqref{1} and, equivalently, \eqref{3} may be considered difference analogues of the Maclaurin expansion for functions of $n$ variables.
Representations \eqref{1} and, equivalently, \eqref{3} are essentially unique: the values of $G_K(\x_K)$ are uniquely determined for all $\x$ and $K$ such that $p_K(\x)\ne0$. This can be verified by reasoning quite similar to \eqref{2}; this is left as an exercise for now.
