I was able to derive some closed form expressions using Mathematica for a binary loop Ising model, but they are unwieldy even for 3-loop and get complicated very fast with increasing loop size.
The basic idea is to use the fact that gradient of log Z gives marginal probabilities, so you invert that mapping algebraically, plug into the original equation of the joint and simplify.
If we have
For example, here's probability of 3-loop ising taking configuration {x1,x2,x3} where m1,m2,m3 are probabilities P(X1=X2),P(X2=X3),P(X3=X1) respectively, after simplification (Mathematica's FullSimplify) http://mathurl.com/324o4a3
For uniform potentials, it doesn't get a lot nicer. For instance take Ising model with uniform potentials on a loop of size n, then probability of all spins being +1 can be written in terms of marginal m = P(x1=x2) as
$exp(n j)/Z$ where Z=$\lambda_1^n + \lambda_2^n$, $\lambda_1=e^j+e^{-j}$, $\lambda_2=e^j-e^{-j}$ and j is the solution of
$\frac{\lambda_1^n(e^{2j}-1)+\lambda_2^n(e^{2j}+1)}{2 \lambda_1 \lambda_2}/Z=m$
Mathematica can solve the above equation for various values of n, but the expression becomes large very fast.
