Q:

Gibbs Sampling for Boltzmann Machines

David Mac Kay, in his book on machine learning talks about Boltzmann machines, and on pg. 3 here http://www.inference.phy.cam.ac.uk/itprnn/ps/521.526.pdf
He says "the second equation  $<x_ix_j>_{P(x|W)}$ is not so easy to evaluate, but it can be estimated by Monte Carlo methods, that is, by observing the average value $x_ix_j$ while the activity rule of the Boltzmann machine (43.3) is iterated". 
I found sample code that learns BMs, seen here https://stats.stackexchange.com/a/132555/9577, and in this code a np.random.rand(N)<sigmoid(b) call is used to generate samples from "the model" which represents the previous state of the probability distribution. This implementation makes sense, but I was wondering if this is what Mac Kay meant in the passage above. Maybe he was talking about iterating over values and multiplying by the probability values coming from sigmoid, not generating samples. 
Any ideas? 
Thanks! 

A:

Not an expert on this but I think $\langle x_i x_j\rangle_{P(x|W)}$ is not easy to evaluate exactly because we need to sum over all possible settings of $X:=(x_1,\ldots,x_n)$.
$\langle x_i x_j\rangle_{P(x|W)} = \sum_{x_1}\sum_{x_2}\cdots \sum_{x_n} x_i x_j p(X|W)$
which involves $2^N$ summands. Evaluating the normalizer is also difficult for the same reason.

